Bayesian and Classical Solutions for Binomial Cytogenetic Dosimetry Problem

Abstract: The main interest of the cytogenetic dosimetry is the prevision of an unknown radiation dosage based in cytogenetic analysis. In this paper the dosimetry problem is formulated as a linear calibration problem for binary response data. Two approaches are considered for inference on the quantity of interest, which is expressed as a calibration parameter in a discrete response variable situation. One is based on the maximum likelihood approach, which depends on large sample results and the second one is based on a Markov chain Monte Carlo (MCMC) simulation approach using BUGS. Application to a data set obtained from blood cultures exposed in vitro to Cobalt 60 (.60Co) at the Energetic Nuclear Research Center (IPEN Brasil) is considered.


Introduction
The calibration problem can be briefly described as follows.There are two related responses x and y, where x represents the true value of the characteristic of interest and y a variable related to it.The statistical controlled calibration problem considers that the variables x and y are related through a fuctional form as in the usual regression problem.However, the main interest center on estimation of an unknown value of x.At a first stage a collection of values x 1 , x 2 , . . ., x n are fixed and corresponding values y 1 , y 2 , . . ., y n are observed.At a second stage, k replications of the response are observed, y 01 , y 02 , . . ., y 0k , are observed, corresponding at an unknown value x 0 .The interest is on estimating x 0 given (x i , y i ), i = 1, 2, . . ., n and y 0j , j = 1, 2, . . ., k.The special case where y and x are linearly related and y is normally distributed has been extensively considered.A good exposition of this area is presented in Osborne (1991) and Brown (1993).Extensions for Student-t models and more generally, elliptical linear models are presented in Branco et al. (1998Branco et al. ( , 2000)).
In this paper, it is considered that the response variable y is binomial, with parameter p.In this case, a function of p is modeled as a linear function of the independent variable x.
The cytogenetic dosimetry is concerned with the relationship between dosage as some form of exposure to radiation and response, as some measure of genetic aberration.The exposure is typically very difficult to measure and the cytogenetic analysis can be helpful to make inference about the unknown exposure dosage.The cytogenetic dosimetry experiment in vitro considers samples of cell cultures exposure to a range of doses levels of a given agent.At each dose level some measure of cell disability is recorded.In the application developed at section 3, the agent is gamma radiation and cell disability is the presence of micronuclei (MN) at the cells.We recall that Cobalt 60 (. 60 Co) is a heavy radioactive isope of the mass number 60 produced in nuclear reactors and used as a source of gamma rays (as for radiotherapy).Therefore, k 1 , k 2 , . . ., k n cells are exposed to fixed dose levels d 1 , d 2 , . . ., d n , respectively.The response y 1 , y 2 , . . ., y n are the number of cells with MN, for each dose level.The interest centers on estimating the unknown dose, d 0 , to which k 0 new cells were exposed, based in the number of cells with MN, y 0 , and a dose-response model.See Madruga et al. (1996) for details on the data set.
The model used for describing the dose-response relationship considered in the paper is described in Finney (1971), where it is suggested the use of the logit and probit models to study the problem.Two components are considered: the stimuli (radiation, for example) and the response observed in a subject (blood cells, for example).The denomination dose is used to describe the intensity of the stimuli at which the subject is submit-ted.Tolerance, denoted by T , is the value used to specify the limit of the stimuli, after which a response is expected (cell deformation, for example).Moreover, tolerance is a population characteristic varying with the population units.Given a dose level d, a response is expected in subjects with T ≤ d.Thus, the expected proportion of subjects with positive response is p dt, where g(t) is the probability density function associated with T .Since T is a positive random variable, the transformation X = log T may be considered, taking values in and for which we consider p If f is the normal density then the probit model follows.To establish the calibration problem, let y be the positive response among n subjects submitted to a value x of the independent variable.Considering y|x ∼ Bin(k, p) in the probit model, x and y are related through the nonlinear model with β 1 = −µ/σ and β 2 = 1/σ, β 1 ∈ and β 2 ∈ + .Thus, a linear transformation is obtained relating x and a function of p, Φ −1 (p).In the logistic case, the transformation obtained is Estimates obtained by using the logistic or the probit model are similar, except for small (close to zero) or large (close to one) values of p, as considered, for example, in Lloyd (1999).Estimates for β 1 and β 2 can be obtained by using the maximum likelihood approach, which are computed by using numerical techniques.It is available in any statistical software.The Bayesian methodology for analyzing logistic regression models abound in the literature.See, for example, Zellner and Rossi (1984), Albert and Chib (1998) and Bedrick et al. (1997).The above references mainly address the issue of Bayesian inference on the regression coefficients.In this paper, the main interest is focused on the calibration problem which seems not to have been considered in the literature using either classical or Bayesian approaches.
As it happens in the case where interest centers on the regression coefficients, there is no analytical or closed form posteriors for the calibration problem.
Section 2 presents classical (based on the maximum likelihood approach) and Bayesian (based on the MCMC methodology) solutions to the calibration problem under the binomial model for logit and probit link functions.The problem of model comparisons is also investigated.An asymptotic approximation is considered for the posterior distribution for estimating x.
Finally in Section 3 we present an application to a data set reported in Madruga (1996) on the number of blood cells affected by . 60Co radiation.

The binomial calibration model
In this section, we consider the binomial calibration model, (2.1) i = 0, 1, . . ., n, where β 1 , β 2 and x 0 are unknown parameters and F is a (known) continue distribution function, which has a continue density function f .Note that if F is the distribution function of the standard normal distribution, then the probit model follows and if F is the distribution function of a logistic distribution, then the logit model follows.It follows from (2.1) that the likelihood function can be written as Thus, it is not simple to deal with the likelihood (2.2) in the sense of obtaining explicit expressions for the maximum likelihood estimator (MLE) and for the posterior distribution of x 0 .To overcome this difficulty, two different approximations are considered.One is based on the asymptotic distribution of the MLE and the other approximation is based on the Markov chain Monte Carlo approach to posterior approximation, by using BUGS (Spiegelhater et al. (1995)).

The maximum likelihood approach
It is well known that under certain regularity conditions the distribution of the MLE of (β 1 , β 2 , x 0 ) can be approximated (see Lehmann (1999)) by a normal distribution with mean (β 1 , β 2 , x 0 ) and the covariance matrix as the inverse of the Fisher information matrix evaluated at the MLE.In the following we discuss the derivation of the maximum likelihood estimators for the binomial calibration problem discussed above.As such, considering the reparametrization and taking the logarithm of the likelihood function (2.2), we obtain the log-likelihood given by Let p0 and ( β1 , β2 ) be the MLE of p 0 and (β 1 , β 2 ), respectively.Thus, from (2.3) it follows that p0 = y 0 /k 0 and ( β1 , β2 ) is a function of the calibration data (k i , x i , y i ), i = 1, . . ., n.Note that p0 and ( β1 , β2 ) are independent.To obtain the MLE x0 of x 0 , we note that so that by using the invariance property of the MLE, it follows that , for the probit model, and ii) x0 = (log{y 0 /(k 0 − y 0 )} − β1 )/ β2 , for the logit model.
As mentioned previously, it is known from the likelihood theory for generalized linear models (see Lloyd (1999)) that the MLE of (β 1 , β 2 ) can not be obtained explicitly, and numerical algorithms such as the Newton-Raphson must be used to compute them.Thus, from the MLE of (β 1 , β 2 ), the MLE of x 0 can be computed by using (2.5) and S-Plus subroutines for logit and probit link functions, for example.
The asymptotic variance of the MLE x0 is considered next.Let where θ θ θ = (β 1 , β 2 , p 0 ) and N = n i=0 k i , be the Fisher information matrix corresponding to the log-likelihood function (2.3).Thus, after some algebraic manipulations, it can be shown that (2.6) , where , i = 0, 1, . . ., n and f is the density function corresponding to the distribution function F .Note that w(t) = F (t)(1 − F (t)), for the logistic model.
Thus, after some algebraic manipulations we obtain that where Notice that the above asymptotic variance require f to be nonnull on .Further, for large N, λ i ≈ k i /N , i = 0, . . ., n, so that ∆(θ θ θ) can be estimated consistently by ∆( θ θ θ).

The Bayesian approach
As mentioned before, it is not possible to obtain explicit expressions for the posterior distribution of x 0 .In fact, from (2.1), it follows that The last integral is intractable even for logit and probit models or for the case where noninformative or reference priors are considered.So, to overcome such difficulties we consider the MCMC methodology for approximating to the posterior distribution.As is well known, the main idea behind MCMC is to build up a Markovian process whose stationary distribution (with density f ) is the one of interest.Among the MCMC methods, the most popular approach is the Gibbs sampler, introduced in Bayesian inference by Gemman and Gemman (1984) while studying problems related to image processing.The books by Robert and Casella (1999) and Chen et al. (2000) contain a comprehensive review of these methods with applications for logistic regression models.
In the case of the binomial calibration model with probit or logit links, the likelihoods are logconcave (Wedderburn (1976)).So the adaptive rejection algorithm (Gilks and Wild (1992)) can be used and implemented by using the software BUGS developed by Spiegelhater et al. (1995).It is a free software and can be obtained from the world wide web page http : ||www.mrcb su.com.ac.uk|bugs.We specify normal priors for β 1 and β 2 with large variances (flat prior) and x 0 ∼ N(m 0 , v 0 ).This will guarantee a proper posterior distribution.For a more recent discussion about this see, for example, Gelfand and Sahu (1999).
Remark 2.1.Another alternative to approximating the posterior distribution is to consider the normal approximation (see Section 2.1).Under general regularity conditions (Chen et al. (2000)), the posterior distribution of x 0 , can be approximated for large N by the normal distribution where θ θ θ = ( β1 , β2 , x0 ) is the MLE of θ θ θ = (β 1 , β 2 , x 0 ) and ∆(θ θ θ) is the asymptotic variance of √ N(x 0 −x 0 ) (see Section 2.1).Thus, the credibility interval for x 0 coincides with the classical interval that follows by using the normal approximation to the distribution of the MLE x0 .
Another aspect of interest is to decide which of the two link functions is more appropriate for a particular data set.The binomial calibration model with the logistic (probit) link function is denoted by M 1 (M 2 ).The Bayes factor can be computed with the aim of deciding for one of the two models.Let p i (y|θ θ θ i ) and π i (θ θ θ i ), respectively, the distribution of the data y = (y 1 , . . ., y n ) and the prior distribution for the parameter vector θ θ θ i under model M i , i = 1, 2. Thus, the Bayes factor for model M 2 against model M 1 is given by where m i (y) is the marginal (predictive) distribution of y under M i , i = 1, 2. The predictive distribution can be approximated by using Monte Carlo methods (see, for example, Bedrick et al. (1996) and Carlin and Chib (1995)).Because the Bayes factor can be extremely sensitive to the specified prior π(θ θ θ i ) (see, for example, O'Hagan (1995) and de Santi and Spezaferri (1997)), several authors have proposed the use of robust Bayes factors and Partial Bayes factors.One of them is the pseudo Bayes factor which is easy to compute and is implemented in the program BUGS.It was introduced in Geisser and Eddy (1979) (see also Gelfand et al. (1992) and Gelfand and Dey (1994)) and it is based on the conditional predictive densities p(y r |y (r) ), where y (r) = (y 1 , . . ., y r−1 , y r+1 , . . ., y n ).
The pseudo-Bayes factor for model M 1 against model M 2 is .
Using Monte Carlo methods and the fact that we can write (see Gelfand and Dey (1994)) where s is the size of the sample generated by using BUGS from the posterior of θ θ θ.
In our case, y r is independent of y (r) given θ θ θ, so that kr−yr , i = 1, . . ., s and r = 1, . . ., n.The estimates c r (l) = pl (y r |y (r) ) can be plotted against r for l = 1, 2, which together with c(l) = n r=1 c r (l) will give indication of which model to select.

Analysis of cytogenetic data
The data considered in the following is part of the data analyzed in Madruga et al. (1996).The experiment was conducted at the São Paulo Nuclear Institute.Presence of micronucley (MN) indicates cell aberration.We consider here only the presence (or absence) of the MN.Table 3.1 presents the frequency of cells with MN in blood samples from two healthy older subjects, which were exposed to gamma radiation (. 60 Co) with doses 20, 50, 100, 200, 300, 400 and 500.We consider the transformation x i = log(d i ), i = 1, . . ., 8, where d i represents the i-th dose value, which is previously fixed.For each one of the groups a model is specified by considering , where y i is the frequency of MN cells associated with the i-th dose value, k i is the number of cells exposed to dose d i , p i is the probability of a cell exposed to the i-th dose value will present micronucley and F is a distribution function.Using the maximum likelihood approach, three models are considered: the probit, the logit and the Student-t with ν = 8 degree of freedom.The three models are compared by using the mean squared error (MSE) computed using cross validation.The results are presented in the Table 3.2.Note, from the table, that the results for the logit and Student-t with 8 degrees of freedom links are very close which is not unexpected as the logistic distribution is well approximated by a Student-t distribution with 8 degrees of freedom (Mudholkar and George (1978)).We can see that the probit model performs worst according to the MSE criterion.
We recall that the main interest here is not fitting the model but the provision of a new dose associated with a new individual.Madruga et al. (1996) consider a new observed value, that is y 0 = 1117 cells with MN in a total of 2427 evaluated cells.In this case we do not know the value of the dose the cells are exposed to.Relative to Table 3.1, a dose value larger than 500 is suggested, yielding an extrapolation problem.3.4 presents classical and Bayesian point and interval estimates based on the probit and logit link functions.The Bayesian computation is based on normal prior specification for x 0 , with mean m 0 = x and variance v 0 = 10, where x is the average of the fixed value x i s.Also normal prior specifications for β 1 and β 2 are considered with a large variances (10 3 ) (flat prioris).The Gibbs samples were generated by using a program implemented in the software WinBUGS (see Spiegehalter et al. (1995)) with an average time of 46 seconds used to generate a sample of size 90,000 disregarding the 10,000 initial iterations.Convergence was verified by using the Geweke's statistics (Geweke (1992)) and also by looking at the graphics of the generated values.
The classical and Bayesian confidence intervals are somewhat close in length for probit and logit links.However, the estimates obtained with the probit link is larger than the estimates obtained with the logit link.The Student-t link estimates follow the logit link estimates.
Using BUGS, we also computed the conditional predictive densities.As we can see in the Figure 3.2, the logit link performs better than the probit link for the most part of the time.However there is not a uniform best model.

Conclusion
The present paper considers Bayesian and classical approaches for the calibration problem with binomial response under logit and probit links.
A new kind of link is proposed, the Student-t link.The three models are considered with classical, or asymptotic Bayesian solutions.The application in Section 3 suggests us that, results using Student-t link (ν = 8) are close to the results by using logit link.In the MCMC Bayesian solution, it is not straightforward to implement the Student-t model.However, that can be done by introducing latent variables as considered in Branco (1997).The Bayesian approach is very helpful for model comparison as we can see from Figure 3.2.As remarked before, here we consider the response variable as binomial, but in the original data set the response is multinomial.It is under current investigation using the multinomial calibration model and will be reported in future work.Some results in that direction are presented in Branco (1997) and Kottas, Branco and Gelfand (2001).The later one presents a Bayesian nonparametric approach for the multinomial calibration problem.However, in both cases easiness of computational implementation using BUGS is lost and more elaborated programs are required.

Table 3 .
1 Frequency of MN for binucleated cells

Table 3 .
2 Adjusted values for logit, probit and Student-t models

Table
Figure 3.1 Graphics (p i versus d i ) Table 3.4 MLE, posterior mean and 95% C.I. for d 0 , when y 0 = 1117