Radiat Environ Biophys (2014) 53:775–780 DOI 10.1007/s00411-014-0563-4
SHORT COMMUNICATION
A Bayesian hierarchical method to account for random effects in cytogenetic dosimetry based on calibration curves Shuhei Mano • Yumiko Suto
Received: 7 February 2014 / Accepted: 14 August 2014 / Published online: 26 August 2014 Springer-Verlag Berlin Heidelberg 2014
Abstract The dicentric chromosome assay (DCA) is one of the most sensitive and reliable methods of inferring doses of radiation exposure in patients. In DCA, one calibration curve is prepared in advance by in vitro irradiation to blood samples from one or sometimes multiple healthy donors in considering possible inter-individual variability. Although the standard method has been demonstrated to be quite accurate for actual dose estimates, it cannot account for random effects, which come from such as the blood donor used to prepare the calibration curve, the radiationexposed patient, and the examiners. To date, it is unknown how these random effects impact on the standard method of dose estimation. We propose a novel Bayesian hierarchical method that incorporates random effects into the dose estimation. To demonstrate dose estimation by the proposed method and to assess the impact of inter-individual variability in samples from multiple donors on the estimation, peripheral blood samples from 13 occupationally non-exposed, non-smoking, healthy individuals were collected and irradiated with gamma rays. The results clearly showed significant inter-individual variability and the standard method using a sample from a single donor gave anti-conservative confidence interval of the irradiated dose. In contrast, the Bayesian credible interval for irradiated dose calculated by the proposed method using samples from multiple donors properly covered the actual doses. S. Mano (&) Department of Mathematical Analysis and Statistical Inference, The Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa 190-8562, Japan e-mail: [email protected] Y. Suto Research Center for Radiation Emergency Medicine, National Institute of Radiological Sciences, Chiba, Japan
Although the classical confidence interval of calibration curve with accounting inter-individual variability in samples from multiple donors was roughly coincident with the Bayesian credible interval, the proposed method has better reasoning and potential for extensions. Keywords Dicentric chromosome assay Calibration Random effects Bayesian hierarchical model
Introduction Radiation exposure causes DNA strand breaks that lead to chromosomal aberrations. Among radiation-induced chromosomal aberrations, multicentric chromosomes represented by dicentric chromosomes are considered sensitive and specific bio-indicators for assessing radiation dose. For more than three decades, the dicentric chromosome assay (DCA) has been considered the ‘gold standard’ for biological dosimetry in accidental radiation exposure, in which a dicentric yield per cell from a radiation-exposed patient is referred to a calibration curve (dose–response curve) prepared in advance from dicentric scoring data at several dose points from blood samples of one healthy donor (Barquinerro et al. 1995) or by pooling scoring data from several donors (Thierens et al. 2005) in considering inter-individual variability. In the latter cases, Van Buul and Natarajan (1980) and Kakati et al. (1986) noted a remarkable inter-individual variability. On the other hand, some reports showed moderate or no significant interindividual variability so that data of all donors could be pooled (Vinnikov and Maznyk 2013; Al-Hadyan et al. 2014). In DCA, inter-individual variability seems to be a small confounding factor to variations in dose estimation (Prasanna and Coleman 2013), with the exception of the
123
776
habit of smoking versus non-smoking at the background level (Romm et al. 2009). In the standard method (IAEA 2011; Szluinska et al. 2005), the calibration curve is modeled by Poisson regression, where the number of dicentric chromosomes, yi ; i ¼ 1; . . .; m; are regressed onto prefixed doses, xi, in an in vitro irradiation In Poisson regression, we assume yi Poissonðni li Þ with a quadratic curve, li ¼ ax2i þ bxi þ c, where ni is the number of assessed cells for the prefixed doses xi. Hereafter, we will refer to the data (ni, yi) as training data. Maximum likelihood estimates of the parameters are obtained by the iterated reweighted least square (IRLS) method (McCullagh and Nelder 1989). The point estimates and their standard errors are used to construct confidence intervals of the calibration curve. Exposure dose is inferred by examining a blood sample from a radiation-exposed patient in which the number of assessed cells no and the number of dicentric chromosomes yo are recorded. Hereafter, we will refer to the data (no, yo) as test data. The inference is an inverse estimate, or a calibration problem (Brown 1994). For the linear model, the rigorous theory has been well developed due to its simplicity; however, for complex models such as generalized linear models, which include Poisson regression, the inference causes technical burden. In the standard DCA method, an ad hoc method is employed: the number of dicentric chromosomes is assumed to obey Poisson distribution, in which the mean is fixed to the observed count, and the confidence interval of the incidence of dicentric chromosomes is constructed. By using the calibration curve, the confidence interval of the incidence of dicentric chromosomes is mapped to the dose space (IAEA 2011; Szluinska et al. 2005). Although the standard method has been established and demonstrated to work very well for actual dose estimates for many years, several drawbacks have been recognized. Especially, (1) The inverse estimate in the standard method does not account for uncertainty of the observed count. In reality, the observed count is a single realization of a random variable. (2) The standard method in not flexible to account for uncertainty of the irradiated dose. (3) The Poisson model is too restrictive. In particular, the variance and the mean should be equal in the Poisson model, which could be violated by actual data. Such deviation may occur in a form of underdispersion or, much more frequency, an overdispersion. Bayesian methods have potential to overcome all of these drawbacks. The first serious attempts of Bayesian methods in cytogenetic dosimetry include Groer and Pereia (1987) and Bender et al. (1988) and a recent comprehensive review is in Ainsbury et al. (2013a). In the Bayesian methods, the inverse estimate is replaced by the posterior credible interval, the uncertainty of the irradiated dose is modeled by the prior distribution, and the mixture of the
123
Radiat Environ Biophys (2014) 53:775–780
Poisson distribution by a prior distribution can account for overdispersion. See Ainsbury et al. (2013a, b) for discussion and an implementation to computer software. In this paper, we propose a novel Bayesian method, which in addition to overcoming above drawbacks in standard method, accounts random effects that might be ubiquitous in modern cytogenetic dosimetry. The random effects include the blood donors used to prepare the calibration curve, the radiation-exposed patients, and the examiners. The random effects in the blood donors, or inter-donor variations, have been recognized by practitioner (Van Buul and Natarajan 1980; Kakati et al. 1986; Thierens et al. 1991, 2005; Vinnikov and Maznyk 2013; Prasanna and Coleman 2013; Al-Hadyan et al. 2014). For example, Tierens et al. (1991) proposed to construct 95 % confidence intervals of the calibration curve based on the standard deviation among samples in 10 donors with the normal approximation. Nevertheless, to date, it remains controversial how inter-donor variation impacts on dose estimate. To account for random effects in cytogenetic dosimetry, we propose a novel Bayesian method using data from multiple individuals based on a mixed model. The mixed model can be represented as a hierarchical model and the Markov Chain Monte Carlo (MCMC) method (Gilks et al. 1996) can be used to sample from the posterior distributions. In this paper, we investigate the presence of random effects in multiple blood donors, demonstrate dose estimation by the proposed method, and assess the impact of random effects on the dose estimation by using irradiated peripheral blood samples from 13 healthy individuals (Suto et al. 2010). Extensions to more general setting are also illustrated.
Materials and methods Data collection The present analysis was approved by the Institutional Review Board of the National Institute of Radiological Sciences (Registration No. 05-002). Heparinized peripheral blood samples were obtained from 13 occupationally nonexposed, non-smoking, healthy individuals (ages 22-59; 4 males and 9 females) with full informed consent, as previously described (Suto et al. 2013). The DCA dataset, version 1, dated 30 November 2010 (Suto et al. 2010), was used in this analysis. Briefly, blood samples were dispensed to 3 ml/tube and irradiated with one of seven different doses (0, 0.5, 1, 2, 3, 4, 5 Gy) using a Cobalt-60 source at a dose rate of 0.5 Gy/ min. Peripheral blood mononuclear cells were isolated and
Radiat Environ Biophys (2014) 53:775–780
cultured for 48 h with 0.05 lg/ml colcemid (GIBCO), and chromosome preparations were obtained and stained with Giemsa (Merck) as described previously (Suto et al. 2013). Metaphase images were selected in manual mode and captured in AutoCapt mode using an AXIO Imager, a Z2 (Carl Zeiss) equipped with a CCD camera and Metafer 4 (MetaSystems). According to (IAEA 2011), chromosome aberration recording and dicentric scoring were performed by three examiners until at least 100 dicentric or 1,000 metaphase chromosomes were counted per sample. Metaphases containing 45–47 centromeres with clear shapes were selected. A multicentric chromosome with n centromeres was scored as n-1 dicentrics. Each examiner analyzed metaphases that were assigned to each examiner. Then, the three examiners held meetings to make consensus on every metaphases among the three examiners. The dataset is the result of the consensus. Thus, the interexaminer variability, which is one of the main sources of the uncertainty in practical cytogenetic dosimetry, was deliberately eliminated by the quality assurance/quality control procedure. That allows us to illustrate correctly our new Bayesian hierarchical approach on the example, which accounts only for one source of random effects, i.e., blood donors, and ignore others. Nevertheless, we illustrate an extension of the proposed method to account for random effects in examiners.
777
Fig. 1 A graphical representation of the Bayesian hierarchical model used in the proposed method. Open circles are parameters, filled circles are irradiate doses, and boxes are observed counts. The arrows represent dependencies among them. The left hand side represents test data, while the right hand side represents training data. The training data consist of repeats for doses (i = 1,…,m) and individuals (j = 1,…,n)
Statistical analysis We introduce a mixed model extension of the well-established Poisson regression with a quadratic calibration curve: yij Poisson nij lij with ð1Þ lij ¼ ða þ aj Þx2i þ b þ bj xi þ c þ cj for i ¼ 1; . . .; m; j ¼ 1; . . .; n, where nij, yij are the number of assessed cells and dicentric chromosomes from the training data consist of n individuals, respectively, for the prefixed doses xi (i.e., 0, 0.5, 1, 2, 3, 4, 5 Gy) for the j-th individual. The parameters representing fixed effects (common effects) and random effects (effects assigned to each individual) are denoted by a, b, c and aj, bj, cj, respectively. A hierarchical Bayesian structure was assumed (Gilks et al. 1996), whose graphical representation is given in Fig. 1. The parameters for the fixed effects follow the priors a, b * Normal (0, 1002), c * Uniform (0.1) and the parameters for the random effects follow the priors aj * Normal (0, r2a ), bj * Normal (0, r2b ), cj * Normal (0, r2c ). Hyper-parameters were assumed to follow the priors 1=r2a ; 1=r2b ; 1=r2c * Normal (0,106). To simplify expressions, we have denoted the parameters for fixed and random effects and hyper-parameters as hf, hr and r,
respectively. Note that sometimes conjugate distributions are chosen for prior distributions in literatures because of the analytical tractability (for example, Groer and Pereia 1987), but we choose the above setting since easy to interpret and the analytical tractability is not our concern since all of the computations are done numerically. Under these setting, the joint posterior density of the parameters are given by. p hf ; hr ; rjyij / q yij jhf ; hr pðhr jrÞp hf pðrÞ: ð2Þ It is reasonable to assume that parameters for the random effect of the test data, ho , follow the same prior distribution as that of the training data. Hence, Eq. 2 is used as the prior distribution for the parameters and the joint posterior distribution of the parameters and observed counts becomes. p yo ; hf ; hr ; ho ; rjyij / q yo jhf ; ho ; x pðho jrÞp hf ; hr ; rjyij ; where q yo jhf ; hr ; x is the likelihood of the dose, which is the Poisson probability mass function and x is the dose of the irradiation-exposed patient. The marginal likelihood of the dose of the radiation-exposed patient is provided by an integration.
123
778
Radiat Environ Biophys (2014) 53:775–780
pðyo jx; yij Þ / r q yo jhf ; ho ; x pðho jrÞp hf ; hr ; rjyij dhf dhr dho dr;
where the integration was performed by using the MCMC sampling method (Gilks et al. 1996). By assuming the prior distribution of the dose (denoted by p(x)) follows Uni form(0,10), the posterior density is pðxjyo ; yij Þ / p yo jx; yij pðxÞ. Note that the posterior density obtained here is a natural extension of the posterior density given by Groer and Pereira (1987, see also Appendix of Ainsbury et al. 2013a).
Results Assessment of random effect We analyzed the dataset obtained from blood samples treated with gamma irradiation at 0, 0.5, 1, 2, 3, 4, and 5 Gy collected from 13 individuals. In the present study, we used the DCA dataset version 1 (dated 30 November 2010, Suto et al. 2010). We first assessed the applicability of the Poisson regression model to the dataset by fitting the model to the data from each individual separately using the IRLS method. The resulting deviance residuals are plotted in Fig. 2. No significant deviance residuals were observed (the limit for two-sided 5 % significance level is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v2ð4Þ;0:05=14 ¼ 3:95). However, this fact never rejects the existence of inter-individual variability, given that distinct models were fitted to the data from each individual. When a single model was applied to the whole dataset, we
2
1
0
-1
-2 0
1
2
3
4
5
Dose (Gy)
Fig. 2 Deviance residual when the Poisson regression model without accounting random effects is fitted to the data from each individual. Axes: the x-axis is the radiation dose, and the y-axis is in the units of deviance residuals. The limits for the two-sided 5 % significance level are 3.95
123
Table 1 Analysis of deviance table for the Poisson regression with Eq. 1 v2
DICa
Model
Degrees of freedom
deviance
1
89
16,834.2
15,470.3
Fixed
86
682.040
116.529
684.838
Fixed ? random Saturated
47 0
565.511 481.676
83.835
588.147
a
The deviance information criterion
observed significant deviance residuals. In fact, analysis of deviance showed the presence of highly significant interindividual variability (Table 1), where the model ‘‘1’’ is a Poisson regression with a single parameter (i.e., c in Eq. 1) and ‘‘saturated’’ model is a Poisson regression with parameters lij replaced by the observed proportions of the dicentric chromosomes, yij/nij, in Eq. 1. Although the degrees of freedom should be discounted in hierarchical models (Spiegelhalter et al. 2002), v2 = 116.529 (39 degrees of freedom, p \ 10-8). The deviance information criterion (DIC) (Spiegelhalter et al. 2002), which discounts degrees of freedom, also strongly supported the model with random effects. Performance of dose estimation To demonstrate dose estimation by the proposed method and investigate performance of the dose estimations, we used the dataset (Suto et al. 2013). To mimic actual dose estimation procedure, individual NIRS_12 was assumed to be one of the healthy blood donors for preparing the calibration curve and individual NIRS_11 was assumed to be the radiation-exposed patient. Although this choice of individuals is one of 156 possible choices, nevertheless, the tendency reported here roughly holds in other choice of pairs. By using the data obtained from NIRS_12 a calibration curve was constructed by the standard method implemented in CABAS software (Deperas et al. 2007). Then, assuming the number of dicentric chromosomes found in 1,000 cell counts to be 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, the 95 % confidence intervals for radiation dose were constructed by the standard method with the calibration curve constructed above. The results are shown as solid lines in Fig. 3. In the data obtained from NIRS_11, actual numbers of dicentric chromosomes found in 1,000 cell counts were 0, 22, 80 for 0, 0.5, and 1 Gy dose irradiations, respectively [shown by boxes (h)]. It can be seen that none of the three data points are covered by the confidence interval constructed by assuming the number of dicentric chromosome found in 1,000 cell counts to be 0, 20, and 80. This striking result demonstrates that 95 % confidence
Radiat Environ Biophys (2014) 53:775–780
779
Discussion
1.4 1.2
Dose (Gy)
1 0.8 0.6 0.4 0.2 0
0
10
20
30
40
50
60
70
80
90
Number of dicentrics
Fig. 3 Comparison of interval estimates of radiation dose. The dose estimates in assumed irradiated patient NIRS_11 were obtained using standard method for different numbers of dicentrics found in 1,000 cells. The results are presented as 95 % confidence interval provided by the standard method with training data from NIRS_12 (solid perpendicular lines). The 95 % Bayesian credible intervals provided by the proposed method with training data from other 12 individuals (dashed lines) are also shown. The actual observations are also shown by the boxes (h) for data points of NIRS_11 and the crosses (9) for data points of the other 12 individuals. Axes: the x-axis is the number of dicentrics found in 1,000 cells; the y-axis is the radiation dose. The classic 95 % confidence interval of the calibration curve based on multiple donors (according to Thierens et al. 1991) is presented as dotted curves
intervals determined by the standard method using single individual is anti-conservative. Moreover, actual numbers of dicentric chromosomes found in 1,000 cell counts in other 12 individuals are also shown by crosses (9). It can be seen that many data points are not covered by the 95 % confidence intervals which was constructed by the standard procedure with data obtained from NIRS_12. The 95 % Bayesian credible intervals of the dose were constructed based on the proposed method using training data from 12 individuals other than NIRS_11, which are shown as dashed lines in Fig. 3. Because of accounting for random effects, the credible intervals were significantly broader than the confidence intervals determined by the standard method. The credible intervals covered the actual observations for NIRS_11, as expected. The credible intervals also covered observations for the other 12 individuals, although it is not surprising because the 12 individuals contributed to construction of the credible interval. In addition, the classic 95 % confidence interval of the calibration curve based on multiple donors proposed by Thierens et al. (1991), which was described in Introduction, based on 12 individuals other than NIRS_11 is shown. The classic 95 % confidence interval agreed surprisingly well with the 95 % Bayesian credible interval constructed based on the proposed method.
To incorporate random effects in dose estimation, we propose a novel Bayesian method for dose estimation using training data from multiple individuals based on a hierarchical model. The model is a natural extension of the proposed Bayesian methods in literatures (Ainsbury et al. 2013a), and as far as the authors know, the first attempt of systematic treatment of random effects in terms of a hierarchical model. In addition to the three drawbacks introduced in Introduction (technical difficulty in inverse estimation, uncertainty of the irradiated dose, overdispersion), the proposed method naturally overcomes influence of random effects. The proposed method was applied to the dataset obtained from blood samples irradiated with gamma rays from 13 occupationally non-exposed healthy individuals especially to investigate the random effects in multiple blood donors. The result gives a striking caution: confidence intervals determined using standard methods using a sample from a single donor can be anti-conservative; the confidence intervals may not cover the actual dose with credited probability. However, the interesting finding is that the classic confidence interval based on multiple blood donors proposed by Thierens et al. (1991) showed good agreement with the Bayesian credible interval constructed by the proposed method. Although it is surprising if we consider the great methodological advancement achieved by the proposed method, the proposed method has better reasoning and richer potential for extensions than the classic method. Extensions of the proposed model to account for other random effects are straightforward. For example, if we would like to random effect in the examiners, we replace Eq. 1 by lijk ¼ ða þ aj þ a~k Þx2i þ b þ bj þ b~k xi þ c þ cj þ c~k ; where r~k * Normal(0,r~2a ), b~k * Normal(0,r~2b ), c~k * Normal(0,r~2c ) are the random effects come from the k-th examiner. Other possible direction for improving dose estimation is to collect covariates of the training and test data, including statuses of the blood donors, radiationexposed patients, and examiners. It is also straightforward to incorporate these covariates into the model. Acknowledgments We would like to thank Dr. Masanari Itokawa (Tokyo Metropolitan Institute of Medical Science, Tokyo, Japan) for arranging blood sample collection from volunteer donors and Dr. Momoki Hirai (National Institute of Radiological Sciences, Chiba, Japan) for useful discussions. We would also like to thank the reviewer for many useful comments and suggestions. This work was supported by the Research Organization of Information and Systems (Grant awarded by the president) and by the Institute of Statistical
123
780 Mathematics cooperative research program (Grant No. 2012-ISM CRP -4302). This project was partly supported by International Atomic Energy Agency (IAEA CRP E.3.50.08).
References Ainsbury EA, Vinnikov VA, Puig P, Higueras M, Maznyk NA, Lloyd DC, Rothkamm K (2013a) Review of Bayesian statistical analysis methods for cytogenetic radiation biodosimetry. Radiat Prot Dosimetry. doi:10.1093/rpd/nct301 Ainsbury EA, Vinnikov V, Puig P, Maznyk N, Rothkamm K, Lloyd DC (2013b) CytoBayesJ: Software tools for Bayesian analysis of cytogenetic radiation biodosimetry data. Mutat Res 756:184–191 Al-Hadyan K, Elewisy S, Moftah B, Shoukri M, Alzahrany A, Alsbeith G (2014) Establishing cytogenetic biodosimetry laboratory in Saudi Arabia and producing preliminary calibration curve of dicentric chromosomes as biomarker for medical dose estimation in response to radiation emergencies. 3 Biotech doi: 10.1007/s13205-014-0217-x Barquinero JF, Barrios L, Caballı´n MR, Miro´ R, Ribas M, Subias A, Egozcue J (1995) Establishment and validation of a dose-effect curve for gamma-rays by cytogenetic analysis. Mutat Res 326:65–69 Bender MA, Awa AA, Brooks AL, Evans HJ, Groer PG, Littlefield LG, Pereira C, Preston RJ, Wachholz BW (1988) Current status of cytogenetic procedures to detect and quantify previous exposures to radiation. Mutat Res 196:103–159 Brown PJ (1994) Measurement. Oxford University Press, Oxford, Regression and Calibration Deperas J, Szluinska M, Deperas-Kaminska M, Edwards A, Lloyd D, Lindholm C, Romm H, Roy L, Moss R, Morand J, Wojcik A (2007) CABAS: a freely available PC program for fitting calibration curves in chromosome aberration dosimetry. Radiat Prot Dosimetry 124:115–123 Gilks WR, Richardson S, Spiegelhalter DJ (1996) Markov Chain Monte Carlo in Practice. Chapman & Hall, Florida Groer PG and Pereira CAdeB (1987) Calibration of a radiation detrctor: chromosome dosimetry for neutrons. In: Viertl R (ed) Probability and Bayesian statistics. Plenum Press, New York
123
Radiat Environ Biophys (2014) 53:775–780 International Atomic Energy Agency (2011) Cytogenetic Dosimetry: Applications in Preparedness for and Response to Radiation Emergencies. International Atomic Energy Agency, Vienna Kakati S, Kowalczyk JR, Gibas Z, Sandberg AA (1986) Use of radiation induced chromosomal damage in human lymphocytes as a biological dosimeter is questionable. Cancer Genet Cytogenet 22:137–141 McCullagh P, Nelder JA (1989) Generalized Linear Models, 2nd edn. Chapman & Hall, Florida Prasanna PGS, Coleman N (2013) Cytogenetic Biodosimetry. Medical Consequences of Radiological and Nuclear Warfare, USA. Chapter 12, pp 272 Romm H, Oestreicher U, Kulka U (2009) Cytogenetic damage analyzed by the dicentric assay. Ann Inst Super Sanita 45:251–259 Spiegelhalter DJ, Best NG, Carlin BP, Van der Linde A (2002) Bayesian measures of model complexity and fit. J Roy Stat Soc B 64:583–616 Suto Y et al (2010) Radiation Cytogenetics: dicentric chromosome assay data collection, version 1 (Dated on 30 Nov 2010). Available at: http://www.nirs.go.jp/ENG/core/rem/rem02_1. shtml Suto Y, Hirai M, Akiyama M, Kobashi G, Itokawa M, Akashi M, Sugiura N (2013) Biodosimetry of restoration workers for Tokyo electric power company (TEPCO) Fukushima Daiichi nuclear power station accident. Health Phys 105:366–373 Szluinska M, Edwards AA, Lloyd DC (2005) Statistical Methods for Biological Dosimetry. Health Protection Agency, Didcot Thierens H, Vral A, de Ridder L (1991) Biological dosimetry with the micronucleus assay for lymphocytes: interindividual differences in doseresponse. Health Phys 61:623–630 Thierens H, De Ruyck K, Vral A, de Gelder V, Whitehouse CA, Tawn EJ, Boesman I (2005) Cytogenetic biodosimetry of an accidental exposure of a radiological worker using multiple assays. Radiat Prot Dosimetry 113:408–414 van Buul PP, Natarajan AT (1980) Chromosomal radiosensitivity of human leucocytes in relation to sampling time. Mutat Res 70:61–69 Vinnikov VA, Maznyk NA (2013) Cytogenetic dose–response in vitro for biological dosimetry after exposure to high doses of gammarays. Radiat Prot Dosimetry 154:186–197
SHORT COMMUNICATION
A Bayesian hierarchical method to account for random effects in cytogenetic dosimetry based on calibration curves Shuhei Mano • Yumiko Suto
Received: 7 February 2014 / Accepted: 14 August 2014 / Published online: 26 August 2014 Springer-Verlag Berlin Heidelberg 2014
Abstract The dicentric chromosome assay (DCA) is one of the most sensitive and reliable methods of inferring doses of radiation exposure in patients. In DCA, one calibration curve is prepared in advance by in vitro irradiation to blood samples from one or sometimes multiple healthy donors in considering possible inter-individual variability. Although the standard method has been demonstrated to be quite accurate for actual dose estimates, it cannot account for random effects, which come from such as the blood donor used to prepare the calibration curve, the radiationexposed patient, and the examiners. To date, it is unknown how these random effects impact on the standard method of dose estimation. We propose a novel Bayesian hierarchical method that incorporates random effects into the dose estimation. To demonstrate dose estimation by the proposed method and to assess the impact of inter-individual variability in samples from multiple donors on the estimation, peripheral blood samples from 13 occupationally non-exposed, non-smoking, healthy individuals were collected and irradiated with gamma rays. The results clearly showed significant inter-individual variability and the standard method using a sample from a single donor gave anti-conservative confidence interval of the irradiated dose. In contrast, the Bayesian credible interval for irradiated dose calculated by the proposed method using samples from multiple donors properly covered the actual doses. S. Mano (&) Department of Mathematical Analysis and Statistical Inference, The Institute of Statistical Mathematics, 10-3 Midori-cho, Tachikawa 190-8562, Japan e-mail: [email protected] Y. Suto Research Center for Radiation Emergency Medicine, National Institute of Radiological Sciences, Chiba, Japan
Although the classical confidence interval of calibration curve with accounting inter-individual variability in samples from multiple donors was roughly coincident with the Bayesian credible interval, the proposed method has better reasoning and potential for extensions. Keywords Dicentric chromosome assay Calibration Random effects Bayesian hierarchical model
Introduction Radiation exposure causes DNA strand breaks that lead to chromosomal aberrations. Among radiation-induced chromosomal aberrations, multicentric chromosomes represented by dicentric chromosomes are considered sensitive and specific bio-indicators for assessing radiation dose. For more than three decades, the dicentric chromosome assay (DCA) has been considered the ‘gold standard’ for biological dosimetry in accidental radiation exposure, in which a dicentric yield per cell from a radiation-exposed patient is referred to a calibration curve (dose–response curve) prepared in advance from dicentric scoring data at several dose points from blood samples of one healthy donor (Barquinerro et al. 1995) or by pooling scoring data from several donors (Thierens et al. 2005) in considering inter-individual variability. In the latter cases, Van Buul and Natarajan (1980) and Kakati et al. (1986) noted a remarkable inter-individual variability. On the other hand, some reports showed moderate or no significant interindividual variability so that data of all donors could be pooled (Vinnikov and Maznyk 2013; Al-Hadyan et al. 2014). In DCA, inter-individual variability seems to be a small confounding factor to variations in dose estimation (Prasanna and Coleman 2013), with the exception of the
123
776
habit of smoking versus non-smoking at the background level (Romm et al. 2009). In the standard method (IAEA 2011; Szluinska et al. 2005), the calibration curve is modeled by Poisson regression, where the number of dicentric chromosomes, yi ; i ¼ 1; . . .; m; are regressed onto prefixed doses, xi, in an in vitro irradiation In Poisson regression, we assume yi Poissonðni li Þ with a quadratic curve, li ¼ ax2i þ bxi þ c, where ni is the number of assessed cells for the prefixed doses xi. Hereafter, we will refer to the data (ni, yi) as training data. Maximum likelihood estimates of the parameters are obtained by the iterated reweighted least square (IRLS) method (McCullagh and Nelder 1989). The point estimates and their standard errors are used to construct confidence intervals of the calibration curve. Exposure dose is inferred by examining a blood sample from a radiation-exposed patient in which the number of assessed cells no and the number of dicentric chromosomes yo are recorded. Hereafter, we will refer to the data (no, yo) as test data. The inference is an inverse estimate, or a calibration problem (Brown 1994). For the linear model, the rigorous theory has been well developed due to its simplicity; however, for complex models such as generalized linear models, which include Poisson regression, the inference causes technical burden. In the standard DCA method, an ad hoc method is employed: the number of dicentric chromosomes is assumed to obey Poisson distribution, in which the mean is fixed to the observed count, and the confidence interval of the incidence of dicentric chromosomes is constructed. By using the calibration curve, the confidence interval of the incidence of dicentric chromosomes is mapped to the dose space (IAEA 2011; Szluinska et al. 2005). Although the standard method has been established and demonstrated to work very well for actual dose estimates for many years, several drawbacks have been recognized. Especially, (1) The inverse estimate in the standard method does not account for uncertainty of the observed count. In reality, the observed count is a single realization of a random variable. (2) The standard method in not flexible to account for uncertainty of the irradiated dose. (3) The Poisson model is too restrictive. In particular, the variance and the mean should be equal in the Poisson model, which could be violated by actual data. Such deviation may occur in a form of underdispersion or, much more frequency, an overdispersion. Bayesian methods have potential to overcome all of these drawbacks. The first serious attempts of Bayesian methods in cytogenetic dosimetry include Groer and Pereia (1987) and Bender et al. (1988) and a recent comprehensive review is in Ainsbury et al. (2013a). In the Bayesian methods, the inverse estimate is replaced by the posterior credible interval, the uncertainty of the irradiated dose is modeled by the prior distribution, and the mixture of the
123
Radiat Environ Biophys (2014) 53:775–780
Poisson distribution by a prior distribution can account for overdispersion. See Ainsbury et al. (2013a, b) for discussion and an implementation to computer software. In this paper, we propose a novel Bayesian method, which in addition to overcoming above drawbacks in standard method, accounts random effects that might be ubiquitous in modern cytogenetic dosimetry. The random effects include the blood donors used to prepare the calibration curve, the radiation-exposed patients, and the examiners. The random effects in the blood donors, or inter-donor variations, have been recognized by practitioner (Van Buul and Natarajan 1980; Kakati et al. 1986; Thierens et al. 1991, 2005; Vinnikov and Maznyk 2013; Prasanna and Coleman 2013; Al-Hadyan et al. 2014). For example, Tierens et al. (1991) proposed to construct 95 % confidence intervals of the calibration curve based on the standard deviation among samples in 10 donors with the normal approximation. Nevertheless, to date, it remains controversial how inter-donor variation impacts on dose estimate. To account for random effects in cytogenetic dosimetry, we propose a novel Bayesian method using data from multiple individuals based on a mixed model. The mixed model can be represented as a hierarchical model and the Markov Chain Monte Carlo (MCMC) method (Gilks et al. 1996) can be used to sample from the posterior distributions. In this paper, we investigate the presence of random effects in multiple blood donors, demonstrate dose estimation by the proposed method, and assess the impact of random effects on the dose estimation by using irradiated peripheral blood samples from 13 healthy individuals (Suto et al. 2010). Extensions to more general setting are also illustrated.
Materials and methods Data collection The present analysis was approved by the Institutional Review Board of the National Institute of Radiological Sciences (Registration No. 05-002). Heparinized peripheral blood samples were obtained from 13 occupationally nonexposed, non-smoking, healthy individuals (ages 22-59; 4 males and 9 females) with full informed consent, as previously described (Suto et al. 2013). The DCA dataset, version 1, dated 30 November 2010 (Suto et al. 2010), was used in this analysis. Briefly, blood samples were dispensed to 3 ml/tube and irradiated with one of seven different doses (0, 0.5, 1, 2, 3, 4, 5 Gy) using a Cobalt-60 source at a dose rate of 0.5 Gy/ min. Peripheral blood mononuclear cells were isolated and
Radiat Environ Biophys (2014) 53:775–780
cultured for 48 h with 0.05 lg/ml colcemid (GIBCO), and chromosome preparations were obtained and stained with Giemsa (Merck) as described previously (Suto et al. 2013). Metaphase images were selected in manual mode and captured in AutoCapt mode using an AXIO Imager, a Z2 (Carl Zeiss) equipped with a CCD camera and Metafer 4 (MetaSystems). According to (IAEA 2011), chromosome aberration recording and dicentric scoring were performed by three examiners until at least 100 dicentric or 1,000 metaphase chromosomes were counted per sample. Metaphases containing 45–47 centromeres with clear shapes were selected. A multicentric chromosome with n centromeres was scored as n-1 dicentrics. Each examiner analyzed metaphases that were assigned to each examiner. Then, the three examiners held meetings to make consensus on every metaphases among the three examiners. The dataset is the result of the consensus. Thus, the interexaminer variability, which is one of the main sources of the uncertainty in practical cytogenetic dosimetry, was deliberately eliminated by the quality assurance/quality control procedure. That allows us to illustrate correctly our new Bayesian hierarchical approach on the example, which accounts only for one source of random effects, i.e., blood donors, and ignore others. Nevertheless, we illustrate an extension of the proposed method to account for random effects in examiners.
777
Fig. 1 A graphical representation of the Bayesian hierarchical model used in the proposed method. Open circles are parameters, filled circles are irradiate doses, and boxes are observed counts. The arrows represent dependencies among them. The left hand side represents test data, while the right hand side represents training data. The training data consist of repeats for doses (i = 1,…,m) and individuals (j = 1,…,n)
Statistical analysis We introduce a mixed model extension of the well-established Poisson regression with a quadratic calibration curve: yij Poisson nij lij with ð1Þ lij ¼ ða þ aj Þx2i þ b þ bj xi þ c þ cj for i ¼ 1; . . .; m; j ¼ 1; . . .; n, where nij, yij are the number of assessed cells and dicentric chromosomes from the training data consist of n individuals, respectively, for the prefixed doses xi (i.e., 0, 0.5, 1, 2, 3, 4, 5 Gy) for the j-th individual. The parameters representing fixed effects (common effects) and random effects (effects assigned to each individual) are denoted by a, b, c and aj, bj, cj, respectively. A hierarchical Bayesian structure was assumed (Gilks et al. 1996), whose graphical representation is given in Fig. 1. The parameters for the fixed effects follow the priors a, b * Normal (0, 1002), c * Uniform (0.1) and the parameters for the random effects follow the priors aj * Normal (0, r2a ), bj * Normal (0, r2b ), cj * Normal (0, r2c ). Hyper-parameters were assumed to follow the priors 1=r2a ; 1=r2b ; 1=r2c * Normal (0,106). To simplify expressions, we have denoted the parameters for fixed and random effects and hyper-parameters as hf, hr and r,
respectively. Note that sometimes conjugate distributions are chosen for prior distributions in literatures because of the analytical tractability (for example, Groer and Pereia 1987), but we choose the above setting since easy to interpret and the analytical tractability is not our concern since all of the computations are done numerically. Under these setting, the joint posterior density of the parameters are given by. p hf ; hr ; rjyij / q yij jhf ; hr pðhr jrÞp hf pðrÞ: ð2Þ It is reasonable to assume that parameters for the random effect of the test data, ho , follow the same prior distribution as that of the training data. Hence, Eq. 2 is used as the prior distribution for the parameters and the joint posterior distribution of the parameters and observed counts becomes. p yo ; hf ; hr ; ho ; rjyij / q yo jhf ; ho ; x pðho jrÞp hf ; hr ; rjyij ; where q yo jhf ; hr ; x is the likelihood of the dose, which is the Poisson probability mass function and x is the dose of the irradiation-exposed patient. The marginal likelihood of the dose of the radiation-exposed patient is provided by an integration.
123
778
Radiat Environ Biophys (2014) 53:775–780
pðyo jx; yij Þ / r q yo jhf ; ho ; x pðho jrÞp hf ; hr ; rjyij dhf dhr dho dr;
where the integration was performed by using the MCMC sampling method (Gilks et al. 1996). By assuming the prior distribution of the dose (denoted by p(x)) follows Uni form(0,10), the posterior density is pðxjyo ; yij Þ / p yo jx; yij pðxÞ. Note that the posterior density obtained here is a natural extension of the posterior density given by Groer and Pereira (1987, see also Appendix of Ainsbury et al. 2013a).
Results Assessment of random effect We analyzed the dataset obtained from blood samples treated with gamma irradiation at 0, 0.5, 1, 2, 3, 4, and 5 Gy collected from 13 individuals. In the present study, we used the DCA dataset version 1 (dated 30 November 2010, Suto et al. 2010). We first assessed the applicability of the Poisson regression model to the dataset by fitting the model to the data from each individual separately using the IRLS method. The resulting deviance residuals are plotted in Fig. 2. No significant deviance residuals were observed (the limit for two-sided 5 % significance level is qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v2ð4Þ;0:05=14 ¼ 3:95). However, this fact never rejects the existence of inter-individual variability, given that distinct models were fitted to the data from each individual. When a single model was applied to the whole dataset, we
2
1
0
-1
-2 0
1
2
3
4
5
Dose (Gy)
Fig. 2 Deviance residual when the Poisson regression model without accounting random effects is fitted to the data from each individual. Axes: the x-axis is the radiation dose, and the y-axis is in the units of deviance residuals. The limits for the two-sided 5 % significance level are 3.95
123
Table 1 Analysis of deviance table for the Poisson regression with Eq. 1 v2
DICa
Model
Degrees of freedom
deviance
1
89
16,834.2
15,470.3
Fixed
86
682.040
116.529
684.838
Fixed ? random Saturated
47 0
565.511 481.676
83.835
588.147
a
The deviance information criterion
observed significant deviance residuals. In fact, analysis of deviance showed the presence of highly significant interindividual variability (Table 1), where the model ‘‘1’’ is a Poisson regression with a single parameter (i.e., c in Eq. 1) and ‘‘saturated’’ model is a Poisson regression with parameters lij replaced by the observed proportions of the dicentric chromosomes, yij/nij, in Eq. 1. Although the degrees of freedom should be discounted in hierarchical models (Spiegelhalter et al. 2002), v2 = 116.529 (39 degrees of freedom, p \ 10-8). The deviance information criterion (DIC) (Spiegelhalter et al. 2002), which discounts degrees of freedom, also strongly supported the model with random effects. Performance of dose estimation To demonstrate dose estimation by the proposed method and investigate performance of the dose estimations, we used the dataset (Suto et al. 2013). To mimic actual dose estimation procedure, individual NIRS_12 was assumed to be one of the healthy blood donors for preparing the calibration curve and individual NIRS_11 was assumed to be the radiation-exposed patient. Although this choice of individuals is one of 156 possible choices, nevertheless, the tendency reported here roughly holds in other choice of pairs. By using the data obtained from NIRS_12 a calibration curve was constructed by the standard method implemented in CABAS software (Deperas et al. 2007). Then, assuming the number of dicentric chromosomes found in 1,000 cell counts to be 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, the 95 % confidence intervals for radiation dose were constructed by the standard method with the calibration curve constructed above. The results are shown as solid lines in Fig. 3. In the data obtained from NIRS_11, actual numbers of dicentric chromosomes found in 1,000 cell counts were 0, 22, 80 for 0, 0.5, and 1 Gy dose irradiations, respectively [shown by boxes (h)]. It can be seen that none of the three data points are covered by the confidence interval constructed by assuming the number of dicentric chromosome found in 1,000 cell counts to be 0, 20, and 80. This striking result demonstrates that 95 % confidence
Radiat Environ Biophys (2014) 53:775–780
779
Discussion
1.4 1.2
Dose (Gy)
1 0.8 0.6 0.4 0.2 0
0
10
20
30
40
50
60
70
80
90
Number of dicentrics
Fig. 3 Comparison of interval estimates of radiation dose. The dose estimates in assumed irradiated patient NIRS_11 were obtained using standard method for different numbers of dicentrics found in 1,000 cells. The results are presented as 95 % confidence interval provided by the standard method with training data from NIRS_12 (solid perpendicular lines). The 95 % Bayesian credible intervals provided by the proposed method with training data from other 12 individuals (dashed lines) are also shown. The actual observations are also shown by the boxes (h) for data points of NIRS_11 and the crosses (9) for data points of the other 12 individuals. Axes: the x-axis is the number of dicentrics found in 1,000 cells; the y-axis is the radiation dose. The classic 95 % confidence interval of the calibration curve based on multiple donors (according to Thierens et al. 1991) is presented as dotted curves
intervals determined by the standard method using single individual is anti-conservative. Moreover, actual numbers of dicentric chromosomes found in 1,000 cell counts in other 12 individuals are also shown by crosses (9). It can be seen that many data points are not covered by the 95 % confidence intervals which was constructed by the standard procedure with data obtained from NIRS_12. The 95 % Bayesian credible intervals of the dose were constructed based on the proposed method using training data from 12 individuals other than NIRS_11, which are shown as dashed lines in Fig. 3. Because of accounting for random effects, the credible intervals were significantly broader than the confidence intervals determined by the standard method. The credible intervals covered the actual observations for NIRS_11, as expected. The credible intervals also covered observations for the other 12 individuals, although it is not surprising because the 12 individuals contributed to construction of the credible interval. In addition, the classic 95 % confidence interval of the calibration curve based on multiple donors proposed by Thierens et al. (1991), which was described in Introduction, based on 12 individuals other than NIRS_11 is shown. The classic 95 % confidence interval agreed surprisingly well with the 95 % Bayesian credible interval constructed based on the proposed method.
To incorporate random effects in dose estimation, we propose a novel Bayesian method for dose estimation using training data from multiple individuals based on a hierarchical model. The model is a natural extension of the proposed Bayesian methods in literatures (Ainsbury et al. 2013a), and as far as the authors know, the first attempt of systematic treatment of random effects in terms of a hierarchical model. In addition to the three drawbacks introduced in Introduction (technical difficulty in inverse estimation, uncertainty of the irradiated dose, overdispersion), the proposed method naturally overcomes influence of random effects. The proposed method was applied to the dataset obtained from blood samples irradiated with gamma rays from 13 occupationally non-exposed healthy individuals especially to investigate the random effects in multiple blood donors. The result gives a striking caution: confidence intervals determined using standard methods using a sample from a single donor can be anti-conservative; the confidence intervals may not cover the actual dose with credited probability. However, the interesting finding is that the classic confidence interval based on multiple blood donors proposed by Thierens et al. (1991) showed good agreement with the Bayesian credible interval constructed by the proposed method. Although it is surprising if we consider the great methodological advancement achieved by the proposed method, the proposed method has better reasoning and richer potential for extensions than the classic method. Extensions of the proposed model to account for other random effects are straightforward. For example, if we would like to random effect in the examiners, we replace Eq. 1 by lijk ¼ ða þ aj þ a~k Þx2i þ b þ bj þ b~k xi þ c þ cj þ c~k ; where r~k * Normal(0,r~2a ), b~k * Normal(0,r~2b ), c~k * Normal(0,r~2c ) are the random effects come from the k-th examiner. Other possible direction for improving dose estimation is to collect covariates of the training and test data, including statuses of the blood donors, radiationexposed patients, and examiners. It is also straightforward to incorporate these covariates into the model. Acknowledgments We would like to thank Dr. Masanari Itokawa (Tokyo Metropolitan Institute of Medical Science, Tokyo, Japan) for arranging blood sample collection from volunteer donors and Dr. Momoki Hirai (National Institute of Radiological Sciences, Chiba, Japan) for useful discussions. We would also like to thank the reviewer for many useful comments and suggestions. This work was supported by the Research Organization of Information and Systems (Grant awarded by the president) and by the Institute of Statistical
123
780 Mathematics cooperative research program (Grant No. 2012-ISM CRP -4302). This project was partly supported by International Atomic Energy Agency (IAEA CRP E.3.50.08).
References Ainsbury EA, Vinnikov VA, Puig P, Higueras M, Maznyk NA, Lloyd DC, Rothkamm K (2013a) Review of Bayesian statistical analysis methods for cytogenetic radiation biodosimetry. Radiat Prot Dosimetry. doi:10.1093/rpd/nct301 Ainsbury EA, Vinnikov V, Puig P, Maznyk N, Rothkamm K, Lloyd DC (2013b) CytoBayesJ: Software tools for Bayesian analysis of cytogenetic radiation biodosimetry data. Mutat Res 756:184–191 Al-Hadyan K, Elewisy S, Moftah B, Shoukri M, Alzahrany A, Alsbeith G (2014) Establishing cytogenetic biodosimetry laboratory in Saudi Arabia and producing preliminary calibration curve of dicentric chromosomes as biomarker for medical dose estimation in response to radiation emergencies. 3 Biotech doi: 10.1007/s13205-014-0217-x Barquinero JF, Barrios L, Caballı´n MR, Miro´ R, Ribas M, Subias A, Egozcue J (1995) Establishment and validation of a dose-effect curve for gamma-rays by cytogenetic analysis. Mutat Res 326:65–69 Bender MA, Awa AA, Brooks AL, Evans HJ, Groer PG, Littlefield LG, Pereira C, Preston RJ, Wachholz BW (1988) Current status of cytogenetic procedures to detect and quantify previous exposures to radiation. Mutat Res 196:103–159 Brown PJ (1994) Measurement. Oxford University Press, Oxford, Regression and Calibration Deperas J, Szluinska M, Deperas-Kaminska M, Edwards A, Lloyd D, Lindholm C, Romm H, Roy L, Moss R, Morand J, Wojcik A (2007) CABAS: a freely available PC program for fitting calibration curves in chromosome aberration dosimetry. Radiat Prot Dosimetry 124:115–123 Gilks WR, Richardson S, Spiegelhalter DJ (1996) Markov Chain Monte Carlo in Practice. Chapman & Hall, Florida Groer PG and Pereira CAdeB (1987) Calibration of a radiation detrctor: chromosome dosimetry for neutrons. In: Viertl R (ed) Probability and Bayesian statistics. Plenum Press, New York
123
Radiat Environ Biophys (2014) 53:775–780 International Atomic Energy Agency (2011) Cytogenetic Dosimetry: Applications in Preparedness for and Response to Radiation Emergencies. International Atomic Energy Agency, Vienna Kakati S, Kowalczyk JR, Gibas Z, Sandberg AA (1986) Use of radiation induced chromosomal damage in human lymphocytes as a biological dosimeter is questionable. Cancer Genet Cytogenet 22:137–141 McCullagh P, Nelder JA (1989) Generalized Linear Models, 2nd edn. Chapman & Hall, Florida Prasanna PGS, Coleman N (2013) Cytogenetic Biodosimetry. Medical Consequences of Radiological and Nuclear Warfare, USA. Chapter 12, pp 272 Romm H, Oestreicher U, Kulka U (2009) Cytogenetic damage analyzed by the dicentric assay. Ann Inst Super Sanita 45:251–259 Spiegelhalter DJ, Best NG, Carlin BP, Van der Linde A (2002) Bayesian measures of model complexity and fit. J Roy Stat Soc B 64:583–616 Suto Y et al (2010) Radiation Cytogenetics: dicentric chromosome assay data collection, version 1 (Dated on 30 Nov 2010). Available at: http://www.nirs.go.jp/ENG/core/rem/rem02_1. shtml Suto Y, Hirai M, Akiyama M, Kobashi G, Itokawa M, Akashi M, Sugiura N (2013) Biodosimetry of restoration workers for Tokyo electric power company (TEPCO) Fukushima Daiichi nuclear power station accident. Health Phys 105:366–373 Szluinska M, Edwards AA, Lloyd DC (2005) Statistical Methods for Biological Dosimetry. Health Protection Agency, Didcot Thierens H, Vral A, de Ridder L (1991) Biological dosimetry with the micronucleus assay for lymphocytes: interindividual differences in doseresponse. Health Phys 61:623–630 Thierens H, De Ruyck K, Vral A, de Gelder V, Whitehouse CA, Tawn EJ, Boesman I (2005) Cytogenetic biodosimetry of an accidental exposure of a radiological worker using multiple assays. Radiat Prot Dosimetry 113:408–414 van Buul PP, Natarajan AT (1980) Chromosomal radiosensitivity of human leucocytes in relation to sampling time. Mutat Res 70:61–69 Vinnikov VA, Maznyk NA (2013) Cytogenetic dose–response in vitro for biological dosimetry after exposure to high doses of gammarays. Radiat Prot Dosimetry 154:186–197
