Reanalysis of cancer mortality in Japanese A-bomb survivors exposed to low doses of radiation: bootstrap and simulation methods

Analyses of cancer mortality in survivors of the atomic bombing of Japan in 1945 have played a central role in establishing risk estimates and radiation protection limits for workers and the wider public. While many people exposed at Hiroshima and Nagasaki received fairly high doses, the recommended occupational and environmental dose limits aim to control the risks from much lower doses. Typically, analyses of the full cohort of Japanese survivors, whose doses ranged from 0 - 4000 mSv or more, are extrapolated downwards to predict the risks from doses below 20 mSv [1]. However, even if a particular model gives a good overall account of the full range of data, fitting that model and extrapolating the results may misrepresent the risks at low doses where the dose response may be quite distinct.

My previous paper [2] considered the cancer mortality data for those survivors whose external (flash) dose from the A-bombs was below 20 mSv, the annual occupational dose limit recommended by the International Commission on Radiological Protection (ICRP). I applied several different models, each of which used as their exposure variable the radiation dose lagged by an unknown latency. One model adopted the standard assumption that Excess Relative Risk (ERR) is proportional to dose, in this case the lagged dose. Other models allowed that ERR might be a non-linear function of dose, as suggested by radiation cell biology studies of the bystander effect. All the models detected ERR values for a dose of 10 mSv which were far higher than predicted by extrapolating the standard analyses of the full cohort. The linear model found large, significant results for the liver, while the preferred non-linear model found large, significant results for the stomach, liver, and lung. These results occurred over a range of latencies.

The statistical method I used to assess the significance of the findings with non-linear models was criticised by Mark Little [3]. In this paper the data is re-analysed with other methods of statistical inference for the preferred non-linear model.

The models use Poisson regression to analyse the grouped data in the 0 - 20 mSv subcohort of the 1950-1990 Life Span Study LSS12 [4]. The models take the form λ = λ_o(1+ERR) where λ is the cancer risk, λ_o depends only on control variables, and ERR depends only on the lagged radiation dose. Profile Likelihood Confidence Intervals for ERR were based in [2] on the Likelihood Ratio Test (LRT) for comparisons between nested models. Once latency is fixed, the models depend smoothly on all other parameters. If Wilks Theorem [5] holds, the asymptotic null distribution of LRT is χ² on d degrees of freedom when the null hypothesis is specified by fixing the values of d parameters in the wider model. For example, the LRT comparison of the linear model ERR = βD with the control model ERR = 0 is asymptotically χ² on 1 d.f. if the control model holds. Here the null hypothesis is β = 0 and Wilks Theorem does apply to this nested pair. As Little [3] pointed out, the non-linear models are less well behaved. The transient model ERR = σDexp(-τD) and the two-phase model ERR = βD + σDexp(-τD) are indeterminate in τ when σ = 0, so one of the regularity conditions for Wilks Theorem fails. In any case, even if the asymptotic null distribution of LRT were known, the distribution for a given finite set of records may differ from the asymptotic case.

However, the actual distribution of LRT can be estimated by simulation before comparing two competing models or constructing confidence intervals. Bootstrap methods [6, 7] can also be used to construct confidence intervals more directly.

In response to questions raised during peer review I also consider whether the data justifies fitting the two-phase model, investigate uncertainties in the "optimal latency" as estimated from the data, extend the analysis to include the covariate "city", and compare results from the 0 - 20 mSv and 5 - 500 mSv dose ranges.

An initial analysis of latency is carried out for all cancer sites with at least 100 cases (deaths) in the 0 - 20 mSv subcohort. Sites which appear to have a raised dose response are then analysed by the bootstrap and the distribution of LRT is estimated by simulation. A fuller simulation is carried out for the stomach, the site with the largest number of cases in this subcohort. The analysis focuses on the behaviour of the two-phase model (defined below) and its relation to the linear model. It aims to compare confidence intervals obtained by bootstrap and LRT simulation methods and to discover without relying on Wilks Theorem, whether the apparent elevation and non-linearity of the dose response in the 0 - 20 mSv subcohort are statistically significant. Further analysis considers variation in latency, the impact of "city", male and female subcohorts, and the 5 - 500 mSv dose range. Computations use the statistical freeware R[8], giving a further check on previous results. Possible sources of error in the Japanese data itself and the potential for confounding by other covariates as discussed in [2], are not reconsidered here.

As previously source data for the 1950-1990 mortality cohort LSS12 was obtained from RERF [9] via the CEDR [10], and the 0 - 20 mSv subcohort defined by restricting the weighted adjusted colon dose. The subcohort comprises 3011 data cells with 1690391.75 p-y. The control model, linear model, and two-phase model are defined as in [2]. Briefly, the control model is log-linear in 14 indicator variables for 5 year age-at-exposure categories, an indicator for gender, and the numerical variable log mean attained age, each of which is well-defined for the data cells. The control model is specified as λ = λ₀ = exp(α+Σβ_jx_j) where α and the β_j are unknown parameters, x_j are the covariates, exp denotes exponential, and the Poisson parameter in cell i with T_i person-years is λ_iT_i where λ_i is evaluated using the covariates in cell i. At this stage, as in [2], "city" (Hiroshima or Nagasaki) is not included as a covariate.

Other models are defined by λ = λ₀(1+ERR) where ERR depends only on D_ϕ, the radiation dose lagged by a latency parameter ϕ. D_ϕ is defined as the weighted adjusted colon dose in 10 mSv units when time-since-exposure ≥ ϕ, and 0 otherwise.

For the linear model, ERR = βD_ϕ

For the two-phase model, ERR = βD_ϕ + σD_ϕexp(-τD_ϕ) where exp denotes exponential.

The models require 1+ERR > 0 in all cells so that the Poisson distributions are defined. For the two-phase model τ > 0, slightly modifying the previous approach where τ ≥ 0. However, τ = 0 reduces to the linear model as does σ = 0.

In fitting a model at latency ϕ by Poisson regression, the unknown parameters α,β_j and β, σ, τ (as appropriate to the model) are chosen to maximise the likelihood of the observed data assuming independent Poisson distributions in each cell. Equivalently, fitting minimises the sum over the 3011 data cells of E_i - O_iln(E_i), where O_i and E_i = λ_iT_i are the observed and expected number of cases in cell i.

Latency ϕ is fixed during the minimisation. The models depend smoothly on all other parameters. ERR_D,ϕ denotes the Excess Relative Risk evaluated at D for the fitted model with latency ϕ. Note ERR_1,ϕ is the estimated risk at 10 mSv lagged by ϕ years. For model I nested within model J, LRT_J-I,ϕ is the Likelihood Ratio Test computed at latency ϕ.

Code written in R emulates the previous model fitting with Excel-Solver, using the Newton-Raphson minimisation routine "optim" [11]. Results for the stomach at latency ϕ = 5, 6,... 44 are checked against the values obtained previously.

At a given cancer site, the control, linear, and two-phase models are fitted at latency ϕ = 5, 6,... 44 to give LRT_J-I,ϕ for each pair of nested models: linear against control, two-phase against control, and two-phase against linear. Sites are selected for further analysis if LRT_lin-con,ϕ > 4 for some ϕ, or if LRT_2p-lin,ϕ > 6 for some ϕ (not assuming that this criterion establishes statistical significance). For sites meeting the criterion for the linear model, Profile Likelihood Confidence Intervals for β = ERR_1,ϕ are computed using Wilks Theorem at integer values of ϕ. For sites meeting the criterion for the two-phase model, ϕ = ϕ_mabs is chosen to maximise LRT_2p-con,ϕ (not LRT_2p-lin,ϕ) so that fitting at ϕ_mabs gives the absolute Maximum Likelihood Estimate (MLE) for the model. The "optimal latency" ϕ_m is chosen to maximise LRT_2p-con,ϕ subject to the constraint ERR_1,ϕ ≥ 0. Bootstrap confidence intervals and LRT simulations are then computed at ϕ_m . Details are given in the next two more technical subsections, and examples with commentary are supplied with Additional File 1. The method for taking account of uncertainty in the "optimal latency" is outlined in a final subsection.

Bootstrap

Fitting the two-phase model at ϕ = ϕ_m gives a maximum likelihood estimate for the vector of all 20 model parameters and determines fitted values of λ_i in each data cell. The parametric bootstrap assumes the observed data arose from sampling this fitted model, and considers the variation arising if other data had been sampled from that same model.

Sampling the independent Poisson distributions with parameters λ_iT_i gives simulated observations in each data cell. Fitting the two-phase model to this simulated data gives simulated parameter values including , and , and therefore a simulated value of ERR = ERR_D,ϕ = ). The entire process is repeated B times to give the bootstrap replications ERR*.

There are many well-established methods to obtain confidence intervals from bootstrap replications [7, 12]. The two options chosen here are the original "percentile" method, which is easy to describe, and the "Bias-corrected accelerated" method which is much faster to compute, and in general more accurate.

The 1-2α "percentile" confidence interval has endpoints

where denotes the Bα^th ordered value of ERR*. With B = 1000 replications and α = 0.025, the endpoints of the 0.95 percentile interval are the 25^th and 975^th order statistics of ERR*.

For the Bias-corrected accelerated (BC_a) confidence interval (see [7] Chapter 5 pp. 203-207 and p.249) simulation is again carried out using λ_iT_i. However, the simulated observations are now used to fit the model to the "Least Favourable family" of distributions, obtained by restricting the parameter vector ψ to a line parametrised by ζ, passing through the MLE and defined by where

Here i^-1( ) is the inverse of the Fisher Information evaluated at and is the gradient of the function h which defines the variable of interest (in this case ERR) in terms of the model parameters ψ.

For each bootstrap replication, fitting the simulated data along the LF family determines ζ and thus the simulated value ERR = . If ERR* _LF denotes the resulting set of B replications, the BC_a confidence interval has endpoints

where

Here Φ is the standard normal cumulative distribution and z^(α) is the 100αth percentile point of a standard normal distribution. The parameters w and a are estimated from the bootstrap replicates. The bias- correction w is computed as

i.e. the fraction of bootstrap replications along the LF family which are below the original fitted value of ERR, in normal units.

The acceleration a is computed as

where m₃ and m₂ are the 3^rd and 2^nd moments of , the partial derivative of the log likelihood of the simulated data with respect to ζ, evaluated at ζ = 0.

Code for this calculation is shown in Additional File 1, Code File AC1.txt, which calls Additional File 1, Data Files stomdat1.txt, lambda2.txt and theta2.txt. Additional File 1, Commentary File ACom1.doc discusses this code.

Despite the intricate formulas, computation is much faster for BC_a than for the simpler percentile method, because optimisation is now restricted to a line within the full 20-dimensional parameter space. The percentile CIs are computed at B = 1000, the recommended minimum number, while BC_a results are obtained at B = 5000.

Coverage for the percentile and BC_a methods is tested using the stomach data. For each of 400 sets of simulated observations generated by sampling the Poisson distributions determined by the original fitted parameters, the model is refitted to the simulated data. Confidence intervals for the simulated ERR are obtained by the bootstrap methods (B = 1000) using the refitted parameters and counted as 1 or 0 if they contain or exclude the "true" value of ERR from the original fitted parameters. To streamline these computations of coverage for the percentile method, data is restricted to the Hiroshima subcohort (N = 1536 cells) and the model simplified to use 10-year rather than 5-year age-at-exposure categories. Estimating coverage for the BC_a method uses the full 0 - 20 mSv dataset and model. The application of both methods to compute CIs from the observed data uses the full dataset and model.

To check R against Excel, the two-phase model is fitted to the stomach data for latencies ϕ = 5, 6,... 44 and LRT_2p-con,ϕ computed in both programmes. The difference (LRT_2p-con,ϕ,R) - (LRT_{2p-con, ϕ, Excel}) has a range of (-0.079, 0.000), mean = -0.004, s.d. = 0.014. At the optimal latency ϕ_m = ϕ_mabs = 11.89, LRT_{2p-con, ϕ} = 21.1903 by either method. The parameter estimates in R are = 0.4586, = 120.0076, = 82.6729, confirming the point estimates reported in [2]. On that basis, the methods are taken as comparable. Results are also comparable when computing LRT-based confidence intervals.

For the linear model ERR = βD_ϕ with ϕ = 5, 6,... 44 only the liver and urinary tract showed significant results. Figure 1 shows LRT_lin-con,ϕ against ϕ for these two cancers, while Figure 2 shows the fitted values ERR₁ = when ϕ ≥ 30 with 95%CIs (Profile Likelihood, assuming LRT ~ χ² on 1 d.f. for testing β = β₀ against β ≠ β₀). For the liver with this model, ϕ_m = ϕ_mabs = 38.58, at which = 0.69 (0.25, 1.26). Of the 540 cases in the subcohort, 171 have ϕ > 38.58. For the urinary tract, ϕ is locally optimal at 41.32 at which = 2.14 (0.51, 5.12). Of the 103 cases in the subcohort, 23 have ϕ > 41.32.

The stomach, liver, lung, pancreas, and leukaemia were selected for further analysis as LRT_2p-lin,ϕ > 6 for some ϕ, an indication that the dose response may be non-linear. For the pancreas and leukaemia, Figure 3 shows LRT_2p-lin,ϕ and LRT_2p-con,ϕ while Figure 4 shows ERR_1,ϕ = as estimated from fitting the two-phase model at latency ϕ. Note that may differ greatly from the value obtained by fitting the linear model at ϕ. Related graphs for the stomach, liver, and lung were shown in [2], though their statistical significance is reconsidered below.

Coverage for the BC_a method is estimated for the stomach at 94% for ERR₁ (376 of 400) and 93% for ERR_0.025 (372 of 400). Coverage for the percentile method is estimated for the stomach using the 0 - 20 mSv Hiroshima subcohort with the simplified two-phase model (10 year age-at-exposure categories) giving 94.3% for ERR₁ (377 of 400) and 93% for ERR_0.025 (372 of 400). All CIs reported in Table 2 use the full model (5 year categories, both cities).

For the stomach with latency ϕ = 5, 6,... 21, Figure 5 shows the BC_a 95% confidence intervals (B = 5000) for ERR₁ and ERR_0.025. Each of these CIs is strictly positive. Since the minimum time-since-exposure in the subcohort is 6.08 years, for ϕ = 5 or 6 the model involves no correction for latency and ERR is a function of the unmodified colon dose.

More generally, for each site (stomach, liver, lung, pancreas, leukaemia) there is a neighbourhood surrounding ϕ_m over which BC_a 95% CIs for ERR₁ and ERR_0.025 are strictly positive, which was evaluated restricting to integer latencies and using B = 5000 bootstrap replications. For the stomach, both CIs are strictly positive when 5 ≤ ϕ ≤ 21. For the liver, both CIs are strictly positive when 32 ≤ ϕ ≤ 43. For the lung, both CIs are strictly positive when 7 ≤ ϕ ≤ 21. For the pancreas, both CIs are strictly positive when 5 ≤ ϕ ≤ 11. For leukaemia, the 95% CI for ERR₁ is strictly positive when 24 ≤ ϕ ≤ 26, but the 95% CI for ERR_0.025 is strictly positive when 19 ≤ ϕ ≤ 36.

For each site, simulating data (200 sets) from the fitted model at ϕ_m and re-determining the optimal latency (restricting to integer values) from the simulated data shows the uncertainty in ϕ_m itself. For the stomach 5 ≤ ≤ 21 for 196 out of 200 simulated datasets, while 5 ≤ ≤ 21 for 187 out of 200 simulated datasets. For the liver, 32 ≤ ≤ 43 for 200 out of 200 simulated datasets, while 32 ≤ ≤ 43 for 197 out of 200 simulated datasets. For the lung, 7 ≤ ≤ 21 for 184 of 200, while 7 ≤ ≤ 21 for 173 of 200 simulated datasets (nb for the lung ϕ_m ≠ ϕ_mabs). For the pancreas the results are weaker, as 5 ≤ ≤ 11 for 157 of 200, and 5 ≤ ≤ 11 for 151 of 200 simulated datasets. For leukaemia 24 ≤ ≤ 26 for 132 of 200 and 24 ≤ ≤ 26 for 130 of 200, but 19 ≤ ≤ 36 for 188 of 200 and 19 ≤ ≤ 36 for 184 of 200 simulated datasets.

Additional File 1, Tables S1, S2, and S3 give the simulated null distributions of LRT for the two-phase model with the null hypothesis H₀: (β, σ, τ) = (β₀,σ₀,τ₀) vs the alternative (β, σ, τ) ≠ (β₀,σ₀,τ₀), for the stomach, liver, lung, pancreas, and leukaemia. Additional File 1, Table S1 shows simulation for each cancer at the fitted values (β₀,σ₀,τ₀) = ( ) and for randomly selected points in the neighbourhood N₁₀ of the fitted values (see Methods). Additional File 1, Table S2 shows simulation for each cancer with σ₀ = 0, where τ₀ is arbitrary. Additional File 1, Table S3 shows additional simulation for the stomach when σ₀ ≠ 0, but (β₀,σ₀,τ₀) is outside the neighbourhood N₁₀. For the simulations in S2, Additional File 1, Table S4 shows the results for the linear model and the comparison with the two-phase model. These simulations are carried out at the optimal latency for the specific cancer.

As the control parameters are optimised subject to H₀ before simulating LRT, perhaps other choices of control parameters while retaining H₀ might yield different estimates of the null distribution of LRT subject to H₀. However, four sets of simulations (×500) for the stomach with (β₀,σ₀,τ₀) = (0, 0, arbitrary) using very different choices of control parameters and total cases expected, yielded fairly similar LRT* and fitted gamma distributions as shown in Table 3. Optimised control parameters for the observed data are shown in the last row. This is evidence that the null distribution of is approximately independent of the parameters used to simulate the data, provided that H₀ holds. In the terminology of [7] (Chapter 4 p. 139) is an "approximate pivot" even though Wilks Theorem fails. In turn, this gives the basis for rejecting H₀ given the observed data.

α	β1	β2	β3	β4	β5	β6	β7	β8	β9	β10	β11	β12	β13	β14	β15	β16	N	qlrt	s	r	qgam	expect
-5	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	0	500	8.23	2.38	0.64	8.32	11389.77
-15	1	2	0	0	0	0	0	0	0	0	0	0	0	0	0	0	500	8.61	2.36	0.63	8.50	7528.62
-8	-1	0	2	0	0	0	0	5	0	0	0	0	0	0	0	0	500	8.39	2.35	0.63	8.42	1418.11
-14	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	1	500	8.45	2.24	0.60	8.59	1013.12
-18.13	-0.84	3.11	-1.05	-0.77	-1.09	-0.69	-0.27	-0.31	-0.29	0.20	0.15	0.29	0.40	0.37	0.45	0.16	1000	8.88	2.42	0.60	8.96	1482.00

α,β_j : specified values of control parameters

N: number of simulations

qlrt: 95^th percentile of

s, r: shape and rate parameters of the fitted gamma distribution for

qgam: 95^th percentile of the fitted gamma distribution for

expect: expected number of cases

The simulated null distributions of from Additional File 1, Tables S1, S2, and S3 are then applied to estimate t_.95, and as defined in Methods. All simulations showed LRT and its fitted gamma distribution had 95^th percentiles qlrt(0.95) < 10 and qgamma(0.95) < 10. Thus t_.95 = 10 is a conservative estimate throughout the parameter space, for all cancer sites considered here, and the choice of N₁₀ is appropriate (see Methods).

Likewise from Additional File 1, Table S4 and using the notation in that file, all simulations of the null distribution of LRT_2p-lin showed qlrt2(0.95) < 8 and qgamma2(0.95) < 8 so T_.95 = 8 is a conservative estimate throughout the parameter space, for all cancer sites considered here. For any given site is the supremum of qlrt2(0.95) and qgamma2(0.95) over simulations for that site. is the supremum of qgamma2(0.95) over simulations for that site. Likewise and are the suprema of qgamma2(0.99) and qgamma2(0.999) over simulations for that site. Table 5 shows these critical values and the observed value of LRT_2p-lin (see Table 2) for each site.

Site	ϕm	T.95					LRT2p-lin
Stom	11.89	8	7.542	7.276	10.394	14.708	19.64
Liver	36.90	8	6.779	6.603	9.329	13.082	27.64
Lung	13.60	8	7.418	6.930	9.795	13.740	14.74
Panc	11.86	8	7.372	7.007	10.049	14.265	14.80
Leuk	23.66	8	6.269	5.832	8.507	12.245	25.29

T_.95, , : conservative, global and refined estimates of the 95% critical value of LRT for comparison of two-phase and linear models. See Methods.

, : refined estimates of the 99% and 99.9% critical value of LRT for comparison of two-phase and linear models. See Methods.

LRT_2p-lin : observed value of LRT for comparison of two-phase and linear models. See Table 2.

Null distributions of and LRT_2p-lin for the stomach with ϕ = 5, 10,... 40 show little variation with latency and are similar to those obtained at ϕ_m as shown in Table 6.

ϕ	N	qlrt2pcon	s	r	qgam2pcon	qlrt2plin	sl	rl	qgam2plin
5	500	9.53	2.39 (0.14)	0.59 (0.04)	9.26	7.72	2.07 (0.12)	0.66 (0.04)	7.39
10	500	8.79	2.56 (0.15)	0.65 (0.04)	8.67	6.87	2.09 (0.12)	0.70 (0.05)	7.02
15	500	8.90	2.46 (0.15)	0.64 (0.04)	8.59	6.89	2.02 (0.12)	0.72 (0.05)	6.63
20	500	8.21	2.38 (0.14)	0.63 (0.04)	8.48	6.87	2.00 (0.12)	0.72 (0.05)	6.60
25	500	8.70	2.62 (0.16)	0.65 (0.04)	8.84	7.58	2.01 (0.12)	0.66 (0.04)	7.20
30	500	8.51	2.31 (0.14)	0.63 (0.04)	8.28	6.71	1.84 (0.11)	0.68 (0.05)	6.62
35	500	8.84	2.31 (0.14)	0.61 (0.04)	8.54	7.09	1.79 (0.10)	0.65 (0.04)	6.79
40	500	8.46	2.45 (0.15)	0.67 (0.04)	8.11	6.95	1.86 (0.11)	0.70 (0.05)	6.46
11.89	1000	8.88	2.42 (0.10)	0.60 (0.03)	8.96	7.18	1.94 (0.08)	0.65 (0.03)	7.10

ϕ: latency

N: number of simulations

qlrt2pcon: 95^th percentile of

s, r: shape and rate parameters with their estimated standard errors for the fitted gamma distribution for

qgam2pcon: 95^th percentile of the fitted gamma distribution for

qlrt2plin: 95^th percentile of

sl, rl: shape and rate parameters with their estimated standard errors for the fitted gamma distribution for

qgam2plin: 95^th percentile of the fitted gamma distribution for

Simulation of the optimal latency as described earlier also has implications for non-linearity. Again restricting to integral latencies, for the stomach LRT_2p-lin > 9.4 when 5 ≤ ϕ ≤ 21 and 196 of 200 simulated datasets gave in this range. For the liver LRT_2p-lin > 16.7 when 34 ≤ ϕ ≤ 38 (198 of 200 simulations). For the lung LRT_2p-lin > 7.7 when 7 ≤ ϕ ≤ 21 (184 of 200 simulations). For leukaemia LRT_2p-lin > 9.4 when 17 ≤ ϕ ≤ 36 (191 of 200 simulations). Each of these results shows non-linearity (see Table 5) and for the stomach and liver the evidence is very strong. For the pancreas, LRT_2p-lin > 9.5 when 7 ≤ ϕ ≤ 11 but only 156 of 200 simulations give in this range.

The covariate "city" has little impact on the estimates of ERR when taken as an additional control. For the liver including "city", BC_a 95%CIs for ERR₁ and ERR_0.025 are strictly positive when 29 ≤ ϕ ≤ 43, again restricting to integral latencies. With "city" included in both models LRT_2p-con attains its maximum subject to ERR₁ ≥ 0 when ϕ_m = 37 and LRT_2p-con = 28.41. At ϕ_m = 37, ERR₁ = 1.26 with BC_a 95%CI (0.61, 2.14) and ERR_0.025 = 0.68 with BC_a 95%CI (0.18, 1.25). For comparison, without "city" the optimal integral latency is ϕ_m = 37, at which ERR₁ = 1.22 with BC_a 95%CI (0.59, 2.11) and ERR_0.025 = 0.74 with BC_a 95%CI (0.23, 1.43). As a test of non-linearity (with "city" included in both models) LRT_2p-lin > 6.8 when 29 ≤ ϕ ≤ 41, and LRT_2p-lin > 18.15 when 34 ≤ ϕ ≤ 38, whereas without "city" LRT_2p-lin > 6.86 when 32 ≤ ϕ ≤ 41 and LRT_2p-lin > 16.73 when 34 ≤ ϕ ≤ 38.

A limited analysis by simulation of LRT for the model with "city" was carried out over the 0 - 20 mSv and 5 - 500 mSv dose ranges, focused at the optimal (integral) latencies on the null distributions of LRT_2p-con and LRT_2p-lin at (β₀,σ₀,τ₀) = (0, 0, arbitrary) and (β₀,σ₀,τ₀) = ( , 0, arbitrary) where is obtained by fitting the linear model. Results are shown in Table 8 and the fitted gamma distributions are used to assign p-values to the observed LRT in Table 7.

Dose range	Site	Φm	LRT2p-con Shape	LRT2p-con Rate	t 0.95	LRT2p-lin shape	LRT2p-lin rate	T 0.95
0-20	Stom	11	2.38	0.59	9.06	1.98	0.65	7.25
5-500	Stom	8	2.03	0.62	7.72	1.30	0.55	6.42
0-20	Liver	37	2.50	0.73	7.63	1.69	0.64	6.58
5-500	Liver	37	1.62	0.49	8.39	1.14	0.53	6.16
0-20	Lung	16	2.47	0.58	9.43	2.09	0.68	7.20
5-500	Lung	16	1.81	0.61	7.26	1.38	0.66	5.58
0-20	Panc	10	2.35	0.61	8.65	2.13	0.71	6.93
5-500	Panc	42	1.72	0.52	8.18	1.18	0.60	5.61
0-20	Leuk	25	2.11	0.62	7.94	1.59	0.67	6.01
5-500	Leuk	28	1.87	0.62	7.34	1.25	0.56	6.14

gamma distributions are fitted to 500 simulations of LRT_2p-con for the null hypothesis (β₀,σ₀,τ₀) = (0,0, arbitrary) and of LRT_2p-lin for the null hypothesis (β₀,σ₀,τ₀) = ( ,0, arbitrary) where is the fitted parameter for the linear model. Models include "city" as a control covariate and are fitted at the optimal (integral) latency Φ_m for the site and dose range.

shape: shape parameter of the fitted gamma distribution

rate: rate parameter of the fitted gamma distribution

t _0.95 : critical value of the fitted gamma distribution for

T _0.95 : critical value of the fitted gamma distribution for

From Table 8 (with "city") and Additional File 1, Table S4 (without "city"), the criterion LRT_2p-lin > 6 is a test of non-linearity at 90% level. Accordingly, Table 7 also shows the latency ranges surrounding Φ_m for which LRT_2p-lin > 6.

Other aspects of the analysis (simulated null distributions of LRT_2p-con and LRT_2p-lin across the parameter space and simulated variation in optimal latency ϕ_m) were not repeated with "city".

This paper develops and partially corrects the approach in [2] and the two should be read in conjunction. I focus initially on the 0 - 20 mSv dose range, which provides over 60% of the person-years of observation in the Life Span Study 12 cohort. Over the full data range 0 - 8000 mSv and the truncated range 0 - 4000 mSv the dose response for cancer mortality is known to be approximately linear and risk estimates for 1000 mSv have been obtained through the ongoing Life Span Study project. The ICRP recommendations for radiation protection are then based on linear extrapolation, as the Excess Relative Risk from 10 mSv is presumed to be 0.01 times the ERR from 1000 mSv.

If we were confident that the approximate linearity of the dose response extends to the low dose region, then linear extrapolation would be justified. But this is not known, and the response may be approximately linear over a wide dose range but highly non-linear at lower doses. There are many examples of non-linear dose response in radiation cell biology [16, 17].

Thus, if the Japanese data is used to derive risk estimates for doses such as 10 mSv or 1mSv which are directly relevant to occupational and public exposure, it makes sense to consider the 0 - 20 mSv data in its own right and in relation to wider dose ranges, rather than to pool the data and apply linear extrapolation.

Secondly, a latency parameter ϕ is used to assess the delay between exposure and cancer mortality. Although it would be better to allow latency to modify dose by some smooth function rather than an abrupt switch from no effect to full effect when Time-Since-Exposure passes ϕ, this is much harder to compute and was not attempted here.

As pointed out during peer review, we should consider at the outset whether the data is sufficient for fitting a non-linear model based on mean dose. Figure 6 shows scatterplots of p-y observation against mean dose, for the ranges 0 - 500 mSv, 5 - 20 mSv, 0 - 5 mSv, and 0 - 0.5 mSv. Clearly the 0 - 500 mSv and 5 - 500 mSv data contains a wide variety of doses. The 5 - 20 mSv data shows a spread of mean doses. Whilst the 0 - 5 mSv data is concentrated below 1 mSv, there is a spread of mean doses in that region. Each dose range determines a very large dataset.

The fact that for the stomach, liver, lung, and leukaemia rather similar results for ERR₁ and latency are obtained from fitting the two-phase model across the 5 - 500 mSv dose range as from fitting across the 0 - 20 mSv dose range suggests that there is also sufficient information in the 0 - 20 mSv data to justify fitting this model. The stomach, with the largest number of cases, shows fairly tight confidence intervals for Excess Relative Risk over a broad range of latencies, which supports the validity of the model.

If the dose response were linear, similar estimates of its slope (β) should be obtained over different dose ranges. As the liver and lung data illustrate (Table 10), fitting the linear model separately to subcohorts of the 0 - 20 mSv dose range yields very different estimates of the dose response. For the liver, there are significant results at comparable latencies in the 0 - 20 mSv dose range and all the subranges, but the estimates of β from doses below 5 mSv (or below 0.5 mSv) are completely different from those obtained from doses above 5 mSv (or 0.5 mSv), in fact around 20 times higher when latency = 36 years. For the lung, the linear model detects some significant activity in restricted dose ranges, but the results are very different, and combining the two dose ranges yields insignificant results. By contrast, see Table 7, the two-phase model shows significant results for the lung with the 0 - 20 mSv data for latencies from 10 to 21 years.

It is not unusual for a model to achieve significance on portions of the data but fail to do so when they are combined. Even univariate Ordinary Least Squares regression on a bimodal independent variable may illustrate this, as in Figure 7 with synthetic data.

As discussed in [2] the graphs of ERR against latency show a clear separation into positive and negative regions, as is also apparent with the linear model over restricted dose ranges for the liver and lung. For any given model, the "optimal latency" ϕ_m identifies the time lag which maximises LRT within the region of non-negative ERR₁. The existence of a significant raised risk is of interest whether or not other latencies show significant decrease in risk, with possibly higher LRT. Frequently ϕ_m coincides with ϕ_mabs which maximises LRT without constraint, but the distinction is relevant for the lung (Table 2) and more generally when simulating the variation in latency.

The point estimates of optimal latency ϕ_m and ERR obtained previously are confirmed by Table 2. The next stage of analysis concerns confidence intervals for ERR at the optimal latency. Statistical inference in [2] was based on the Likelihood Ratio Test but it was wrongly assumed, as Little [3] pointed out, that Wilks theorem held for the non-linear models so that LRT would be asymptotically χ² distributed, if the null hypothesis (specifying parameters in the model) were true.

One of the regularity conditions required for Wilks Theorem is that distinct values of the model parameters produce distinct probability distributions. In the two-phase model when σ = 0, ERR = βD_ϕ irrespective of the value of τ. The asymptotic properties of models with an indeterminate parameter (which cannot be identified from the distribution) have been analysed [18–21] but I was unable to use these methods and therefore took an intensive computational approach, which applies to the given set of records rather than asymptotically.

The parametric bootstrap is an established technique for finding confidence intervals in the absence of any information on the distribution of LRT. The percentile method refits the model to simulated data generated by samples from the model whose parameter values were obtained by fitting to the actual data. The BC_a method adjusts for the fact that simulation may produce biased estimates of the statistic of interest (in this case ERR) and of its variance. Both methods gave good coverage with almost identical rates when tested on the stomach data, and comparable confidence intervals for all the sites analysed. The BC_a method is much faster and generally preferable.

However, the underlying assumption of the parametric bootstrap is open to question. Perhaps the observed data arose from sampling at parameter values other than those obtained by fitting the model to that data. Confidence Intervals based on LRT avoid this assumption but depend on the critical value t_.95 used to reject the null hypothesis if LRT > t_.95.

If Wilks Theorem had applied, then LRT for testing H₀: (β, σ, τ) = (β₀,σ₀,τ₀) vs. the alternative (β, σ, τ) ≠ (β₀,σ₀,τ₀) would have the asymptotic distribution LRT ~ χ² on 3 d.f. if H₀ holds, for any H₀. Thus t_.95 = F^-1(0.95) = 7.8147 where F is the cumulative χ² on 3 d.f. would be the appropriate critical value, rejecting H₀ if LRT > t_.95.

In the absence of any theoretical prediction, the null distribution of LRT can still be estimated by Monte Carlo simulation but may depend on the choice of H₀ as seen in Additional File 1, Tables S1, S2, S3 and S4. In fact Additional File 1, Tables S2 and S3 show the null distribution of departs from χ² on 3 d.f. as σ₀ → 0, while Additional File 1, Table S4 shows that the null distribution of LRT_2p-lin departs from χ² on 2 d.f. everywhere. Table 3 indicates that the simulated null distribution of is fairly independent of the parameters for the control covariates, and depends only on H₀: (β, σ, τ) = (β₀,σ₀,τ₀). The null distributions when σ₀ = 0 are quite different from those when σ₀ >> 0, as expected since it is σ₀ = 0 which causes the indeterminacy of the model and the failure of Wilks Theorem on this hyperplane. Gamma distributions give a reasonable fit to the simulated null distribution of in general (see Additional File 1, Tables S1, S2 and S3). When H₀ has σ₀ = 0, typical shape and rate parameters are s ~ 2.24 and r ~ 0.63, while for σ₀ >> 0 s ~ 1.5 and r ~ 0.5 as predicted by Wilks Theorem. Figure 8, the probablity plot for the stomach with (β₀,σ₀,τ₀) = (0,0, arbitrary) shows the failure of Wilks Theorem and the close fit of simulated = LRT_2p-con with a gamma distribution (s = 2.418, r = 0.604).

The Kolmogorov-Smirnov test against the fitted distribution is extremely sensitive and even when it distinguishes LRT from the fitted gamma, a probability plot may show close agreement. For example for the liver, with (β₀,σ₀,τ₀) = (7.297, 0, arbitrary), the K-S test has D = 0.058 with simulated p = 0.001 but the plot (Figure 9) is fairly close.

Likewise, fitted gamma distributions generally give a reasonable fit to the null distribution of LRT_2p-lin with typical parameters s ~ 1.79, r ~ 0.68.

For the linear model, Additional File 1, Table S4 shows simulation of the null distribution of generally conforms to a χ² distribution on 1 d.f. as expected from Wilks theorem. However as β₀ → -0.5, the null distribution approaches a γ distribution with shape = 0.5, rate = 1. Since β₀ = -0.5 specifies a boundary of the parameter space, the theory of mixture models [22] predicts, at the boundary, the null distribution of ~ 1/2 χ² (1 d.f.) = 1/2 γ (0.5,0.5) = γ (0.5,1). The confidence intervals shown in Figure 2 for liver and urinary cancers involve β values far from the boundary and Wilks theorem is valid in this region.

For all simulations and all cancer sites (Additional File 1, Tables S1, S2 and S3) the 95^th percentile of the fitted gamma distributions and the 95^th percentile of the simulated itself are all below 10. Additional File 1, Tables S1, S2, and S3 summarise nearly 50,000 simulations of LRT at a wide variety of randomly chosen values of H₀ . On that basis, t_.95 = 10 is chosen as a conservative estimate applying throughout the parameter space for all cancer sites considered here. H₀ can be rejected whenever > 10, i.e. whenever (β₀,σ₀,τ₀) is outside the neighbourhood N₁₀ of the fitted parameters ( ). N₁₀ is thus an appropriate neighbourhood, as defined in Methods.

Inside this neighbourhood Additional File 1, Table S1 shows gamma distributions generally fit the simulated very well, justifying the definition of from fitted gamma distributions within N₁₀. The intermediate choice is based on simulations inside and outside N₁₀ and uses both the simulations and the fitted distributions to adjust the critical value upwards as either may be underestimates.

For the stomach, liver, and lung, confidence intervals for ERR obtained by the bootstrap (Table 2) or from simulation of (Table 4) are broadly similar, whichever critical value or bootstrap method is chosen. LRT intervals using show closer agreement with the bootstrap methods. For these 3 cancers, all the LRT-based intervals for ERR₁ and ERR_0.025 exclude 0 as do the bootstrap intervals, which are tighter. This validates the claim in [2] that ERR₁ is significantly elevated for the stomach, liver, and lung.

All five cancers show significant non-linearity of the dose response by LRT comparison of the two-phase and linear models (Table 5). For the stomach, liver, and leukaemia the two-phase model is preferable with p << 0.001, while for the lung and pancreas p < 0.01.

The previous discussion concerns behaviour at the optimal latency ϕ_m, an unknown parameter estimated from the data. BC_a intervals, comparable at ϕ_m to CIs obtained by other methods, were used to investigate ERR across a range of latencies, restricted to integers 5, 6,... 44. Each cancer site has latency ranges in which BC_a CIs for ERR₁ and ERR_0.025 are strictly positive. Figure 5 illustrates this for the stomach and also shows that ERR is significantly raised when the dose is not modified by latency (ϕ = 6 is below the minimum time-since-exposure in the 0 - 20 mSv cohort). With the linear model, liver and urinary cancers are significantly elevated over a range of long latencies, showing similar variation (Figures 1 and 2).

Fitting the two-phase model at ϕ_m, simulating new data and finding the simulated optimal latency gives an estimate of the probability that the true value of ϕ_m lies in any given region. For the liver and stomach, this approach gives very strong evidence that ϕ_m lies in the regions where ERR₁ and ERR_0.025 are strictly positive, and where the two-phase model is clearly superior (by LRT) to the linear model. For the lung, the evidence on these points is still valid, though weaker. For leukaemia, there is strong evidence that ϕ_m lies in the region where the dose response is non-linear, and likewise in the region where the 95% CI for ERR_0.025 is strictly positive. However, with the variation in ϕ_m we cannot be confident that the 95% CI for ERR₁ is strictly positive.

Many results for the pancreas are considerably weaker than for other sites. Even at the optimal latency, the confidence intervals are extremely wide. With the BC_a method, the lower limit is attained as the minimum bootstrap value, and is assigned a probability > 0.025. The "bias" and "acceleration" of the BC_a bootstrap sample are both large, unlike those for other sites. When the uncertainty in the optimal latency is determined by simulation at ϕ_m we cannot be confident that the true optimal latency for the pancreas lies within the range where ERR > 0 or where non-linearity is detected by LRT_2p-lin.

Up to this point, the primary aim was to test the validity of the conclusions in [2], and accordingly "city" was not included as a control covariate and the dose range remained at 0 - 20 mSv. A limited reanalysis including "city" in the baseline model and comparing the 0 - 20 mSv and 5 - 500 mSv dose ranges is shown in Tables 7 and 8. The higher dose range was chosen to include 10 mSv, but to exclude the lowest dose category in the LSS12 dataset. Above 5 mSv the revised dosimetry DS02 generally agrees with the DS86 dosimetry used here (and in [2]). Comparing these ranges tests whether the results from 0 - 20 mSv simply reflect some anomaly in the lowest dose category. With the exception of the pancreas, there is a striking similarity between results from the 0 - 20 mSv and 5 - 500 mSv dose ranges, which show comparable estimates of the optimal latency, latency ranges, and the Excess Relative Risk from 10 mSv.

I do not know if any other approach to baseline risk may explain these effects. However, significant non-linearity and positive ERR are also found for the stomach, liver and lung when the data is restricted to male or female subcohorts, with "city" included in the model (Table 9). Any possible interaction effects involving gender would disappear with the all-male or all-female data. As mentioned in [2] the use of log mean attained age was an adequate alternative to the full set of attained age categories, and alternative controls affected the estimates of ERR by no more than a factor of 2.

The results at 0.25 mSv (ERR_0.025) show that the risk detected here arises at extremely low doses. If there are doubts about the accuracy of the DS86 dosimetry, the question remains as to why the DS86 dose is a significant predictor of risk.

Overall, the results have several possible interpretations. They may reflect confounding by other covariates not shown in the public dataset, which may be correlated with the external dose used here (and in LSS12) and may also interact with radiation and may vary over time. Differences between the rural and urban populations may contribute to the non-linear effects. Or, these may be caused by gross underestimates of the external dose (unlikely) or by internal doses arising, for example, from "black rain". Such interpretations would be specific to the Japanese cohort, and would therefore cast doubt on using this data to predict radiation response elsewhere.

However, the results may also indicate that low doses of radiation do impact on cancer mortality, provided that latency is included in the model and the dose response is not constrained to be linear. In radiation cell biology, non-linear dose response is known for the bystander effect using a variety of endpoints and the two-phase model was proposed as a simplified form of a model for the bystander effect [16].

This paper does not reconsider any of these possible interpretations [2]. However, I hope it establishes the statistical validity of this analysis of the available public data for the 1950-1990 mortality cohort.

This reanalysis validates the main conclusions of [2]. Bootstrap methods and Monte Carlo simulation of LRT show that Excess Relative Risk for cancer mortality in Japanese A-bomb survivors exposed to external radiation doses below 20 mSv is positive, large, and significant for various cancers, as detected by models incorporating latency. The dose response is highly non-linear for the stomach, liver, lung, pancreas, and leukaemia. In each case the two-phase model shows large Excess Relative Risk at 10 mSv external dose lagged by the optimal latency for the cancer. LRT-based 95% Confidence Intervals are strictly positive, except for leukaemia. Bootstrap BC_a 95% CIs are strictly positive for all five cancers over a range of latencies. Large, positive Excess Relative Risk is also found when the male and female data is analysed separately for the stomach, liver and lung.

When the optimal latency varies by simulation, the stomach, liver, lung and leukaemia still show non-linear dose response, and likewise ERR > 0 at 95% level for the stomach, liver, lung.

With "city" included as a control covariate in the two-phase model, similar estimates of latency and ERR at 10 mSv are obtained for the stomach, liver, lung, and leukaemia whether the dose range is 0 - 20 mSv or 5 - 500 mSv. Dose response for the liver, lung, and leukaemia is significantly non-linear in the 5 - 500 mSv range, particularly for the lung. Such results cannot be explained by any anomaly of the 0 - 5 mSv data.

The linear model finds significant results for the liver and urinary tract over a range of long latencies.

This analysis of cancer mortality in Japanese A-bomb survivors exposed to low doses of external radiation in 1945 shows significant non-linearity of dose response and significant large Excess Relative Risk over a range of latencies. These findings do not support the current ICRP recommended annual occupational dose limit of 20 mSv.