Dataset Sources and Treatments

Datasets. Six publicly available cohort datasets were analyzed: two Techa River cohorts (incidence and mortality) (9, 21), the INWORKS worker study (22), and three A-bomb survivor studies (7, 8, 23). Analysis began with pairs of counts and person-years of observation obtained for a sequence of dose categories. For the three A-bomb survivor studies, the needed pairs were extracted from their publicly available, stratified datasets (24–26). Data for the Life Span Study (LSS) cohorts were chosen for analysis. Data were grouped either directly or by adjusting person-years for covariates (standardization), using mean corrections determined in the respective studies, as discussed later in this section. For the default analyses of the A-bomb datasets, the maximum dose included was ∼0.6 Gy to keep the focus on the lower dose range and to better match dose ranges in the other studies analyzed. Even in studies with large regression sample sizes, the dose ranges can be restricted to the lower regions when looking for nonlinear effects at low doses (8).

Pairs of counts and person-years for the INWORKS study (22) and both the Techa River incidence cohort study (9) and the mortality cohort study (21) were obtained from data found in the published literature. For ease of reproducibility and to show the differences in reference dose levels, the 6 sets of paired data are collected in the Supplementary Material ( tables S1 (rare-199-04-03_s02.docx)– S3 (rare-199-04-03_s02.docx);  https://doi.org/10.1667/RADE-21-00217.1.S2). The number of dose categories (the regression sample size) was, 11, 7 and 7, respectively, for the INWORKS and the two Techa River studies. In contrast, the A-bomb datasets have 14–16 dose categories over the same dose range and around 25 over the total dose range.

The Techa River datasets had an unusual pattern of counts per category, where the counts in the lowest category (reference category) were much lower than the counts in the next two higher dose categories. To explore the impact of such a reversal on the frequency of false positives, modifications were made to the datasets for sensitivity analysis. The count in the reference category was scaled upward to match the count in the second dose category, while simultaneously scaling person years to keep risk the same.

For the 2017 A-bomb dataset, we chose the Not-In-City group (NIC) as the reference level to avoid questions about the reference level being above a possible dose region of hormesis (27, 28). Only the 2017 dataset separates data for the NIC group. Because such a reference level is not standard, sensitivity analyses were performed. For example, the Not-In-City group was combined with the group of survivors at closer distances, namely those located between 3 and 10 km from the epicenter at the time of the bombings.

It was pointed out to us by a reviewer that, in an earlier set of simulations for the two-stage clonal model, which can produce highly discontinuous dose response functions, reliable parameter estimates for the dose response function simulated were only obtained by individualized regression. Parameter estimates for this model were not reliable when individuals were grouped into Poisson dose categories (29). Although the dose response models considered in this article are all continuous in dose, the likelihoods for the threshold and 2-slope spline models can have discontinuous derivatives. We therefore looked to see if individualizing risks would make a difference in the results for these two dose-response models.

To this end, two new datasets were created with 17,500 simulated Techa River individuals to be used in dose threshold analysis and 1,750 to be used in the 2-slope spline analysis. Each simulated individual was assigned a dose and a survival history, which incorporated radiation-induced mortality. The radiation contribution to survival was generated by a Cox proportional hazards model, with hazard function linear in dose ( Supplementary text S-1 (rare-199-04-03_s01.pdf);  https://doi.org/10.1667/RADE-21-00219.1.S2).

Linearization of Risks

The data were linearized by adjusting person-years of observation at each dose category, which allows the risk at any dose category to be arbitrarily set, while preserving relative variance of risk. The unadjusted risk at each dose category in a cohort study is determined by the number of cancer cases per person-year at risk (counts/pyr). The relative error (standard deviation over the mean) in these risk values depends only on the count variable, if any experimental error in person-years is neglected, as is standard. To obtain a synthetic set of pure linear data with Poisson errors and therefore relative errors matching the underlying study, the counts in each dose category were left unchanged. Only the number of person years in a dose category was modified, which allowed changes in risk without change in count. Person-years were scaled at each dose category, k= 1 to n, as in Eq. (1), so that the ratio of count to person-years increased linearly with the category's dose at index k, designated as d_k.

The variable s = slope. Transformation stages for each dose category in the 2017 A-bomb dataset are shown in  Supplementary tables S-2 (rare-199-04-03_s02.docx) and  S-3 (rare-199-04-03_s02.docx) ( https://doi.org/10.1667/RADE-21-00219.1.S2). After linearization, the numbers and positions of the dose categories and the corresponding counts and count variances remain. The process is equivalent to assigning every individual study subject in a dose category the same value of person-years multiplied by a dose dependent factor whose inverse is linear in dose.

Once the risk was set to pure linearity of dose response, Monte Carlo simulation techniques were used (30) to quantify prevalence of nonlinearity. Poisson variations were added to the count data using the open source R-function, rpois (31). Each of 5,000 unique dose-response curves was generated for each of the 6 linearized datasets by drawing from the Poisson distribution centered on the average count data in each dose category. Random samples of the modified data were then drawn and used to compute test statistics whose distributions provided estimates of the test statistic under the null (30). Sensitivity runs were made with the same starting random number seed, so that the same replications were analyzed for each different set of assumptions.

Person-years were kept unchanged from replication to replication. The choice of linear slope is a free parameter that can be chosen to match published values or varied to see if changes in slope affect the results of analysis. A published linear slope may have been affected by some underlying and uncaptured nonlinearity, and, of course, its value is subject to statistical and other errors so it may deviate from other studies.

Relative Risk Data

In addition to the absolute risks discussed above, we also explored a nonstandard relative risk model, where the data were made relative by dividing all risk values by the reference risk (risk at the lowest dose). This model is not based on a reparameterization that leaves the data alone but makes the parameters relative. Such a reparameterization was not considered here, because it does not produce results that differ from the absolute risk model in our simple case without covariates (results not shown). To distinguish this type of relative risk model from the reparametrized versions, it is labeled here as a “relative risk data model.” An analyst might consider such a relative risk model for exploring nonlinearity, which is easy to analyze with well-established and accessible software routines, should proprietary software, such as AMFIT, not be accessible for analyzing the reparametrized version. Alternatively, an analyst might simply be exploring alternatives to standard approaches.

Consideration of Study Covariates and Standardization

Only univariate analysis with dose as the independent variable was carried out in this article. In the case of the A-bomb datasets, it was possible to adjust risks by important sex-averaged covariates, such as age at exposure and age of attainment. This was possible because the necessary sex-averaged adjustment functions for relative risk models (not absolute risk models) and their fitted coefficients were included in the published A-bomb articles ( Supplementary table S-4 (rare-199-04-03_s02.docx);  https://doi.org/10.1667/RADE-21-00219.1.S2). Sex-averaged coefficients for absolute risk models were not available for all three A-bomb datasets, so the coefficients for relative risk models were used.

As a sensitivity test, these adjustment factors were used to produce standardized person-years for each stratum prior to linearization. Use of these covariate-adjusted risks in linearization of the datasets produced a change in absolute risk scale factor only. There was no change in relative risk at any dose category, a result that was initially surprising to us but was confirmed by algebraic analysis ( Supplementary text S-2 (rare-199-04-03_s01.pdf);  https://doi.org/10.1667/RADE-21-00219.1.S2), as well as confirmed numerically by running simulations and obtaining the same false finding frequencies.

The fact that standardization does not change the linearization by more than a scale factor is a fortunate result, because no covariate adjustment functions are publicly available for either the Techa River Cohort or the INWORKS worker study. Neither could we adjust background by city and location for the A-bomb data, as has routinely been done in studies carried out at the Radiation Effects Research Foundation (RERF). We also were unable to adjust for smoking in the 2017 dataset as was done in the corresponding LSS study (8). Although the relative variation between our crude and standardized risk estimates for the A-bomb datasets amounted to less than 10% in the 0 to 0.6 Gy range (results not shown), such a modest change would not be expected for the INWORKS study, where worker age would be correlated with cumulative dose.

Methods for Question 1 (Are there High False Positive Rates for a Single Test?)

Indicators of Nonlinearity. All of the five indicators of nonlinearity used or inferable from smoothed graphs in the six published studies were tested in each replication of linearized data. Specifically, we catalogued dose thresholds, curvature and quadratic terms in the replicates, as well as the dose response shapes, supralinearity and sublinearity (including hormesis). The tests are described further in Table 1. Individual test results were tabulated and aggregated for each replicate. For null-hypothesis testing, a (false) positive was counted when 95%-confidence intervals for the relevant parameter excluded the linear null.

TABLE 1

Overview and Description of Indicators of Nonlinearity Assessed in Replicates^a

A separate count was made for false-positive curvatures having absolute magnitudes greater than 1 and thus a greater than twofold difference in risk at 1 Gy compared to what would be expected from the slope at low doses. These we deem to be of regulatory interest.

Perhaps, the most important finding of a nonlinearity for regulatory policy would be evidence of a dose threshold. The threshold model used here is a simple, 1-parameter extension of the linear default:

Where d is dose, D is the threshold dose value, B is the slope above threshold, and Y is the risk of cancer.

The case of most interest is when D is found to have a narrow confidence range that excludes the null of 0. However, in addition to counting false positives for D, tabulations were made for all D values above zero in the simulations, no matter how weak the statistical evidence. This was done for several reasons. First, results about dose thresholds that do not rise to the level of a false positive are routinely published with P values, forming a collection that could be used for meta-analysis. Second, dose threshold findings are sometimes presented or discussed in the literature as indicators of nonlinearity without considering the strength of the evidence (32, 33), making it of interest to quantify the likelihood of such indications being statistical fluctuations. It can be argued that such unsupported claims are so obviously flawed that no quantitative response is needed. However, dismissal out of hand may not be sufficient to carry the day with audiences who influence regulation; quantitative illustrations of the frequency of false indications may be helpful.

Average values for the threshold dose, which might serve as a measure of bias away from the null, were also calculated. In an attempt to gain insight into above-zero threshold values, separate counts were made for those threshold dose values that satisfied the AIC criterion for an improved fit, which penalizes the statistical likelihood for each extra parameter introduced (34).

Evaluation of Nonlinear Shapes

Nonlinear shapes beyond the 3-parametric functions, threshold, quadraticity and curvature, were addressed using four methods. The methods were loess smoothing, 2-slope spline, the occurrence of five points lying on one side of the linear line, and multi-model averaging. A 2-slope spline dose response function was of particular interest, because it was used in the Techa River mortality study. For AIC selection and multi-model averaging in our analyses, the 2-slope spline function provided the only shape analysis. The 2-slope spline fits were made to each replicate using the R-functions, ‘lSpline’ and optimizer,‘mle.’

The 5-consecutive-point method was of interest because such dose-response curves tend to look nonlinear. They are an example of what has been called the “clustering illusion” (35). To qualify as a nonlinearity, the required five points had to start with the second data point and cluster on one side of the linear fit to the data, all by at least 0.1 standard deviation.

Locally estimated scatterplot smoothing (loess) was taken to be the default shape analysis. To qualify as an apparent sublinear or supralinear dose response, the 95%-loess confidence band had to fall below or rise above, respectively, the linear fit to the data, somewhere in the 0–0.25 Gy interval. The loess method was chosen as the default for shape analysis because it is easy to understand and interpret. It also has been used in numerous radiation epidemiologic cohort studies (6–9), including three of those analyzed here. The loess technique is of particular interest in an investigation of false findings, because the method is not recommended for datasets with a small number of data points (36) and could reasonably be expected to have non-standard error rates. Further details are given in  Supplementary text S-3 (rare-199-04-03_s01.pdf) ( https://doi.org/10.1667/RADE-21-00219.1.S2).

Additional details for indications of nonlinearity:

Fitting Cancer Count Data

Each replicate's data were analyzed for the five nonlinearities by software either by T-test, profile likelihood or likelihood ratio test. Software was set to call a result a positive when the relevant test statistic reached 95%-confidence (alpha set to 0.05). Deviations from a 5% false-positive rate, which are known to occur (38), could then be assessed by counting the number of software-assigned positives in simulations. Only if the rate in simulations exceeded 10% did we call attention to it, for two reasons. First, some radiation studies use 10% to delineate positives (22, 39, 40). Second, different tests are available to test essentially the same null hypothesis but can give different results for finite regression sample sizes (41). Finally, test values have their own uncertainty (42). Furthermore, an incorrect assumption about the true value of alpha and the width of the distribution for a test statistic would likely have little practical impact were the deviation modest.

For scoping purposes, brute force step searches involving many thousands of small increments in dose were made to find the maxima and confidence intervals for profile likelihoods predicted by software. This profile approach allowed graphical visualization of the complexities of profile likelihoods as functions of parameter values, which made it clear that in some cases software algorithms were inappropriate for finding maxima and confidence intervals. As in Neriishi et al. (43), 95% confidence limits were taken at 3.84/2 = 1.92 units above and below the parameter value at the log of the maximum profile likelihood (44).

Most of the nonlinear models have only one more parameter than the linear default, so the profile likelihood and likelihood ratio tests were generally the same. For the 2-slope spline model, there were two extra parameters, so lower confidence bounds for likelihood ratio tests were found at 5.99/2 units below the maximum loglikelihood.

When it was not necessary to perform step searches, e.g., for a linear-quadratic dose response, the maximum likelihoods were found using standard optimization software, specifically maximum likelihood estimator, “mle” in the R-statistical language, which sped up calculations. For the linear-quadratic fits, false-positive designations were based on the standard errors produced by the software, which are based on assumptions of smooth behavior at maximum likelihood. However special treatment was necessary to use standard software optimizers to speed up calculations for threshold breakpoint and 2-slope spline, which can have multiple likelihood peaks and discontinuous derivatives. For these two functions, threshold and 2-slope spline, segmented searches between dose category boundaries were used to obtain maximum likelihoods and critical points for profile likelihood, as well as for use in likelihood ratio tests. Between datapoints, the likelihood met regularity conditions allowing searches with the mle function to be used to find local maxima within the segments. These calculations were also checked with brute force step searches.

Additional details of the fitting process are listed below.

Methods for Question 2 (What is the Expected Rate of at Least One False Positive in a Single Study?)

The number of replicates with at least one false positive among the five tests for nonlinearity were simply counted for the 5,000 replicates of study data. This accounted for any correlations between the five tests.

Binomial Calculations for a Single Study

With the standard choice of 95% confidence bands, false positives for a single test of nonlinearity will happen in ∼5% of the simulations, assuming required assumptions about the data and modeling are met. The probability, then, that at least one false positive will show up when k independent tests for nonlinearity are made in a single study can be obtained by a binomial calculation and compared to simulation rates. With k set to 5, the number of tests of nonlinearity available in the A-bomb studies that we considered, the result is 23% assuming independence of the tests.

Methods for Question 3 (What is the Expected Rate of at Least one False Positive in Multiple Studies?)

The generalization of Eq. (3) to multiple studies is given in Eq. (4), where f_j is the aggregate false positive rate in the jth of n studies.

For the case of five independent statistical tests of nonlinearity in 6 studies, where the rate in each study is 23%, Eq. (4) predicts a 79% chance that at least one false positive would show up in 6 independent studies. This triplet of numbers, 5%, 23%, and 79% is to be compared with modified values discussed later, when the assumption of independent, perfect tests is relaxed. There are two cases to consider. In the first case it is assumed that all five tests of nonlinearity are available for all 6 studies, which is the assumption that is made to generate the numbers in most of the tables and graphs to follow. The total number of tests considered across all studies, then, is 30. In the second case, which is used to compare simulation results with actual study results, the number of test results available in each specific study was used to determine the reduced f values to insert in Eq. (4). In both cases, the f values were determined from simulations, which automatically accounts for correlation between the different statistical tests of nonlinearity.

To improve the match with actual results in the second case, the dose range considered was increased for the A-bomb datasets to 0–2 Gy, and weighted linear regression of absolute data was used to preserve linearity in this dose range. Poisson regression of absolute data was used for the other datasets.

A study-by-study listing of the primary test results available in the 6 studies, which total to 20, is given in  Supplementary table S-5 (rare-199-04-03_s02.docx) ( https://doi.org/10.1667/RADE-21-00219.1.S2). A primary test means a test performed on all data, without any subgrouping. To fully compare simulation results to A-bomb findings, which included a positive finding for nonlinearity in dose response for males in the 2017 A-bomb study, it was necessary to go beyond primary tests. To this end, linearization was carried out for the 2017 A-bomb dataset separately by sex.

To put the curvature result for males in perspective, a Bonferroni correction of 40 = 2 × 20 tests was used for comparing with simulation results. The factor of 2 comes from the assumption that all of the six study groups would have looked at their test results by sex and would only have reported a result by sex, if it had a very low P value. Thus, we inferred that there might have been tests conditional on sex that did not reach 95% confidence, which led us to double the Bonferroni adjustment to account for this conservative possibility.

AIC Selection of Best Model

One way to avoid choosing a privileged null, as null hypothesis testing does, is to rank models by the highest statistical likelihood found in fits to the data. However, with this method, adding extra parameters to a model will always lead to smaller model residuals. To compensate, a model's fitted likelihood can be penalized by a term related to the number of parameters. A common choice for the likelihood penalty in radiation epidemiology (23, 46–49), which we adopt here, is based on the AIC criterion (34, 50), which in turn has roots in information theory (34). Researchers have other choices, including “corrected” AIC, “consistent” AIC, and Bayesian Information criterion (50), which may change results depending on sample size (50).

To apply the standard AIC method, unity was subtracted from the statistical loglikelihood determined for each fitted parameter. The model generating the highest penalized likelihood (lowest AIC score) was then chosen in each replication of the 5,000 datasets. If AIC selection picked any of the nonlinear dose response functions over the linear default model, a false positive was declared and counted. Although a privileged null is avoided with the AIC penalty method, there is an implicit linear null for nested models like threshold, quadraticity, curvilinearity, and 2-slope spline. In such cases, the AIC criteria produces the same results as would selecting a model using a null-hypothesis test, but with a higher cutoff P value than 0.05, namely 0.157 (50).

Multi-Model Averaging

Multi-model averaging (51) is an example of multi-model inference. The averages were obtained by weighting the five dose-response functions with normalized AIC weights (52). The weights were based on the negative exponentiation of half the AIC values (51, 52), which is algebraically equivalent to taking weights proportional to the likelihood determined by software, reduced by a factor of e^–p, where p is the number of fitted model parameters. The standard error of the multi-model fit at each point on the fitted curve was based on the square root of the weighed sum of the squares of the standard errors produced by fitting software for each model. Multiplication by 1.96 produced the confidence bands.

The final dose-response curve is analogous to a smoothed curve and conceivably might be less likely to produce false positives. However, the choice of a method to compare a multi-model curve to the linear curve used in regulation is not obvious. The approach taken here for illustrative purposes is to require at least a doubling or halving of the crude slope (ERR/Gy) at low doses, while at the same time requiring the 95% confidence band for the fitted curve to exclude the linear fit. Without both of those conditions being met, we presume that regulators would have minimal interest in considering any modifications to regulatory policy.

Some authors narrow down the number of nested models before including them in the mix of models analyzed with AIC weights (52). In our case, we only have nested models to begin with, so we include them all in the model average. Note that some authors pick replicate risk medians of multi-models rather than averages (52). There are, thus, different versions of multi-model inference, which means that our results will not necessarily generalize.

Methods for Question 4 (What is the Expected Frequency of Above-Zero Dose Thresholds?)

Toy Threshold-Dose-Response Model. To aid in understanding frequency results for threshold doses above zero, a toy threshold-dose-response function, with a threshold dose between the first two datapoints, was introduced for use with absolute risk data. The model can be analyzed graphically, without the need for algebraic fitting of data ( Supplementary text S-4 (rare-199-04-03_s01.pdf) and the figure embedded therein;  https://doi.org/10.1667/RADE-21-00219.1.S2).

Fitting Simulated Individualized Survival Data

The Cox proportional hazard model was used to first create individualized linearized datasets using a hazard function linear in dose to predict cancer mortality history. In each replicate, randomization was introduced into the timing of death for each simulated individual. Next, the Cox model was used again, this time to fit the previously created data with nonlinear hazard functions. If a parameter for nonlinearity in a hazard function were declared a positive finding by the software, it would have to be a false positive, because the original dataset was built with a linear function. An embedded flowchart found in  Supplementary text S-1 (rare-199-04-03_s01.pdf) ( https://doi.org/10.1667/RADE-21-00219.1.S2) shows the process.

Two nonlinear hazard functions were used, dose threshold and 2-slope-spline. Hypothetical dose breakpoints in the modeling were chosen at each of the 17,500 individually assigned doses, not between them. The “Surv” function from the “survival” package in the R-statistical language (31) was used to provide a response variable for use in fitting with the R-function, “coxph.” For loess fitting, survival time of cancer deaths was the outcome fitted to dose. To directly compare individualized results with the standard cohort approach, the individualized results were grouped into dose categories and analyzed with Poisson regression.

Except for loess fitting of survival times, it was not feasible to run 5,000 replications of the Cox hazard analyses, as was done in the main analyses in this paper. We found 200 replications were sufficient to compare the individualized results to those obtained for the grouped data. A similar approach was taken with the 2-slope spline model, but the number of persons analyzed was reduced to 1,750.

Results Related to Question 1 (Are there High False Positive Rates for Single Tests?)

False positives rates over 10% occurred for tests of individual dose response functions such as curvature in some nonstandard situations. Thus, there were analysis methods in radiation epidemiology that could overstate evidence for nonlinearity in simulations, but this happened only in two situations. For instance, false positive rates from 13 to 40% occurred across datasets for dose threshold if some standard optimization methods were used to find maximum likelihoods. The excess could be eliminated by using software capable of analyzing data with multiple peaks and discontinuous derivatives in parameter likelihood.

High rates of false positives were also found for some dose-response functions when datasets were analyzed that had a small reference count, even as high as 150. For example, false positives for curvature of relative risk data analyzed with Poisson regression were 66% for Techa River incidence and 28% for Techa River mortality.

Other dose-response functions that gave elevated false positive rates for Techa River were the “5 point test” (5 consecutive points on one side of the linear line). It showed 13% for sublinear and 11% for supralinear when analysis was carried out on relative risk incidence data for Techa River. Even absolute risk data could show elevated rates for Techa River datasets. For example, weighted linear regression of absolute risk data produced false error rates in excess of 10% for the Techa River when it came to threshold (11%), 2-slope spline (16% sublinear, 20% supralinear), and hormesis with the reference being the first data point (38%). The best performing case was Poisson regression of absolute risk data, which had no rates above 10% except for 1 of 3 definitions of hormesis.

A sensitivity analysis demonstrated that it was indeed the low counts, 17 for Techa River incidence data and 150 for the mortality data, that caused the high false positive rates. When the Techa River reference count numbers were artificially increased to match counts in the next highest dose category, 630 and 955 for the two datasets, respectively, all the false positive rates were below 10–11%. Note that in this sensitivity analysis, the reference level person-years were also scaled to keep reference risk values unchanged. Similar problems did not appear in the other four studies analyzed, because the lowest count number in their reference categories was 4,600.

For completeness, results above 10% for all false positives for nonlinearity are listed in  Supplementary table S-6 (rare-199-04-03_s02.docx) ( https://doi.org/10.1667/RADE-21-00219.1.S2). All results, including those below 10%, can be found in  Supplementary tables S-7a (rare-199-04-03_s02.docx)–  S-7d (rare-199-04-03_s02.docx). Count numbers for each dose category in each of the six studies can be found in  Supplementary tables S-1 (rare-199-04-03_s02.docx) and  S-2 (rare-199-04-03_s02.docx). Four examples of loess fitting to the linearized 2017 A-bomb data are presented in  Supplementary fig. S-5 (rare-199-04-03_s01.pdf). Additional examples of false dose response shapes found with loess smoothing, 50 in total, are available in Supplementary file 1 ( https://doi.org/10.1667/RADE- 21-00219.1.S1).

Although almost all of the results for the A-bomb datasets in this article were calculated for a restricted mid-dose region of 0 to 0.6 Gy, it was of interest to assess a higher dose range for comparison. Table 2 gives such a comparison, showing rates when the dose range is extended from 0.6 Gy to 2 Gy. In the extended dose range case, simple Poisson regression will introduce its own nonlinearities, so weighted linear regression was used in the simulations. The resulting false positive rates for individual tests of nonlinearity as shown in the first five rows of Table 2 were similar for both dose ranges. All values were less than 10%.

TABLE 2

False Positive Rates for Nonlinearity, Alone and in Aggregate, for the LSS 2017 A-Bomb Linearized Dataset for Two Dose Ranges, 0 to 0.6 Gy and 0 to 2 Gy (Absolute Risk Data for 5,000 Replications with 95% Confidence Intervals)^a,b

Because the true dose response is not expected to be perfectly linear, without any higher order terms (2), not all nonlinearities would be of regulatory interest. As stated earlier, a curvature of regulatory interest was defined for this paper to be one with magnitude greater than 1, which means that risks at very low doses will be around half the value predicted by the LNT or less. Although this restriction had modest impacts on the frequency of false positives for curvature in the INWORKS, Techa River, and A-bomb datasets restricted to ∼0.6 Gy, there was a large reduction in frequency in simulations when the A-bomb upper dose limit for regression was raised to 2 Gy. The frequencies found for curvature magnitude exceeding unity were 0.16%, 0.3%, and 0.1% for the 2007, 2012, and 2017 datasets, respectively. The rate was also 0.1% for the 2017 dataset, when data for males were analyzed separately. Results for the 0–2 Gy range turn out to be important later, when considering the likelihood of positives for curvature reported in the 2017 A-bomb study.

Results for Individualizing Risk

As was largely the case for Poisson regression of absolute risk data, individualizing risks for the Techa River mortality dataset did not produce any false positive findings above 10% for the two indicators of nonlinearity tested, dose threshold and 2-slope spline.

Linearization Graph

Figure 1 shows a comparison plot of the (standardized) excess relative risks before and after linearizing the 2017 A-bomb incidence study. It can be seen that there is no nonlinearity left in the linearized curve. The “before” data (dashes) in Fig. 1 were standardized to age 70 and age at exposure of 30 years by adjusting person-years to better match published data. Doing so did not change the linearized data beyond a scaling. The “after” data (solid triangles) were obtained by scaling person-time for each data point so that the revised risk lay exactly on a linearized line. No change was made in cancer count data. What is left of the original study are the number and positions of the dose categories and the corresponding counts and count variances, which are used in regressions. The slope of the linearized curve was taken to match the published version computed over the entire dose range. It is slightly lower than the slope obtained by fitting the dose range used (0 to 0.6 Gy). For plotting purposes, the risks were normalized to the lowest fitted risk value.

Results Related to Question 2 (What is the Expected Rate of at Least One False Positive in a Single Study?)

False positives for nonlinearity from null-hypothesis testing were aggregated to obtain a net rate per study replicate, given five tests per study. Figure 2 shows, for each of the 6 datasets considered, the aggregate false-positive rate per study obtained in simulations for the different methods of regression and data treatment. The graph dramatizes the poor performance of relativized risk data when the counts in the reference level are small (Techa River). When the count rates were artificially increased in the sensitivity exercise mentioned earlier, the high aggregate rates for relativized Techa River data dropped to levels comparable to the results for the other datasets shown in Fig. 2 (results not shown).

FIG. 2

The percentage of times that at least 1 false positive showed up in a single dataset, given five tests for nonlinearity per study. Regression sample sizes are highest at the left, lowest on the right. The dashed line at 15% is the average over the six datasets for Poisson regression of absolute risk data. “Weighted” stands for weighted linear regression. The high values for Techa River relative risk data are caused by low count rates in the reference category, which was used as the divisor.

Figure 2. also shows the better performance obtained with Poisson regression of absolute risk data, even for the Techa River datasets with their low count rate in the reference level. In fact, the Techa River datasets performed better than the A-bomb datasets for Poisson regression of absolute data. Why are they better than the A-bomb datasets in this case? The false positive rates for shapes after loess fitting to data made the difference. The Techa River rates were only 30% of the corresponding rates for the A-bomb studies.

The percentage of studies with at least one false positive averaged over all cases in Fig. 2 was ∼25%, which is higher than the 23% predicted by Eq. (3) for independent statistical tests. How can non-independent tests combine to reach and exceed 23%? This happens because in the higher cases, the individual false positive rates are above 5%, compensating for test overlap. Results for single study, aggregate false positive rates that incorporated shape methods other than the default loess fitting, specifically 2-slope spline fits or the five data-point criterion, can be found in  Supplementary tables S-8a (rare-199-04-03_s02.docx) through  S-8d (rare-199-04-03_s02.docx) ( https://doi.org/10.1667/RADE-21-00219.1.S2).

Even without a low reference count, relativized risk data can be a particularly unreliable choice for Poisson regression, because it is a ratio of Poisson counts. The variance in the relative risks is not Poisson and is underestimated in Poisson regression of this type of data, thereby increasing the frequency of false positives.

The average over the six studies for Poisson regression of absolute risk data was 15%. It is shown as the dashed line in Fig. 2. The value of 15% per study is lower than the 23% frequency predicted in Eq. (3) for independent tests with an individual rate of 5%, implying that the tests for nonlinearity are not independent.

AIC selection did not eliminate the multiple comparisons problem associated with null-hypothesis testing of multiple types of dose response function. This was evident with all six datasets, including those whose results are presented in Table 3, specifically, the 2007 and 2017 A-bomb, INWORKS, and Techa River mortality datasets. As shown in Table 3, a nonlinear model was falsely selected as superior to the underlying linear model 24 to 29% of the time, using Poisson regression of absolute data, compared to 15% of the time that null-hypothesis testing falsely selected a positive for nonlinearity averaged over the 6 datasets for the same underlying regression method and data treatment (Fig. 2). Using weighted linear regression of absolute risk data, the corresponding numbers for AIC selection were 35% to 61%, compared to 21% for null-hypothesis testing averaged over the 6 datasets. High frequencies of aggregated false positives were not a surprise given the high effective false-positive cutoff rate (alpha = 0.157) implicit in AIC model selection of nested models. These high frequencies could be reduced using versions of AIC with stronger likelihood penalties (50), such as “corrected” or “consistent” AIC (results not shown). AIC selection avoids the need to explicitly pick a privileged null hypothesis, but it is implicitly a multiple comparison and may make false selections.

TABLE 3

Fraction of Time in 5,000 Replications that False Deviations from Linearity were Found in Linearized Datasets using Aikake Information Criterion (AIC)

Multi-model averaging of nested models with AIC weights did not decrease the frequency of false indications of nonlinearity compared with null-hypothesis testing in most cases, using our two-fold requirement for a false positive. When only the first requirement was met, namely that the confidence band for the average shape excluded the linear line at some dose below 0.25 Gy, the frequency of false positives found for multi-model averaging were similar to the results for AIC selection (results not shown). When the requirement of a major deviation at low doses was added, the fraction of time that a multi-model fit qualified in simulations as a false positive under our definition ranged from 10 to 25% for Poisson regression of absolute data (Table 3). The average per-study result of 17.5% was slightly higher than the corresponding average result of 15% for null-hypothesis testing. The two numbers are not directly comparable, because the nonlinear shape functions incorporated were different, specifically loess fitting for null-hypothesis testing and 2-slope spline for multi-model fitting. Loess fitting tends to pick up more false positives. Regardless, given the uncertainties in definitions and calculations, the conclusion is that both techniques, null-hypothesis testing and multi-model averaging can produce isolated false positives for nonlinearity at comparable rates.

Dataset by dataset numbers for both types of shape analysis, loess and 2-slope spline, are archived in  Supplementary tables S8a (rare-199-04-03_s02.docx)– S8d (rare-199-04-03_s02.docx) ( https://doi.org/10.1667/RADE-21-00219.1.S2) for null hypothesis testing that can be compared to the results for AIC selection and multi-model fitting given in Table 3.

Results Related to Question 3 (What is the Expected Rate of at Least One False Positive in Multiple Studies?)

The expected rate of at least one false positive occurring in multiple studies was obtained by substituting into Eq. (4) the appropriate single-study, aggregate false positive rates (f_j). In the first case considered, it was assumed that all five tests for nonlinearity were available for all six studies, leading to a predicted rate of 61% for Poisson regression of absolute data with dose restricted to a maximum of ∼0.6 Gy. For comparison with actual results in published studies, it was necessary to use f_j based on the actual test results, which were less than five in some studies. In this second case, the chances of at least one false positive occurring in the studies was 50%. Based on this 50% figure, there would be a good chance that the single positive for nonlinearity found for primary tests in the actual studies (8) could be a false positive, at least if results for subgroup data are ignored. The more general case of subgroup findings is considered in the discussion section.

The simulation-derived, per-study rates (f values) that were inserted into Eq. (4) for the first case equaled 0.18, 0.19, 0.19, 0.13, 0.088, and 0.084. The listed sequence begins with the rates for the three A-bomb studies in chronological order, followed by INWORKS and the two Techa River studies, incidence and mortality. The f-values in the second case were, 0.19, 0.11, 0.19, 0.0, 0.085, and 0.06.

Results Related to Question 4 (What is the Expected Frequency of Above-Zero Dose Thresholds?)

An above-zero threshold falsely occurred more than 50% of the time in the simulations for most of the datasets, due to a bias in the threshold dose, which to our knowledge has not previously been discussed. Figure 3a shows the distribution of dose threshold values found in 5,000 simulations of the linearized 2017-A-bomb dataset determined using Poisson regression methods of absolute risk data below 0.6 Gy.

FIG. 3

Histograms of parameter values fitted to 5000 dataset replicates. Threshold dose values are shown in upper panel (3a). Values for the quadratic term in the LQ model are shown in lower panel (3b). All non-zero values are false, failing to capture the underlying linearized dose response. Values were obtained from Poisson regression of absolute risk data for the linearized 2017 A-bomb dataset below 0.6 Gy. The peaks in the upper panel for dose threshold match the dose values assigned to each dose category. About 2/3 of the dose values above zero are contained in the peaks, including the peak very close to zero dose.

Experimental error or statistical variation at a true threshold dose of zero must be zero or positive, but the high frequency of skewing was a surprise before the data were plotted. Percentages ranged from 46% to 69%, depending on the regression method and data treatment, when averaged over the 6 datasets ( Supplementary table S-9 (rare-199-04-03_s02.docx);  https://doi.org/10.1667/RADE-21-00219.1.S2). The percentages were also high (73%), when data were individualized.

Because the datasets had been linearized, all of the above-zero threshold dose values were false, caused by statistical noise. No threshold dose peaking occurred, as expected, in the 200 replications of individualized data simulated for the Techa River mortality study, because no dose categorization was involved.

The peaks for dose threshold in Fig. 3a match the dose values assigned to each dose category. About 2/3 of the threshold dose values above zero are contained in the peaks, including the peak very close to zero dose. Data points in Fig. 3a that qualified as false positives did not occur until doses exceeded the 4th peak at 0.029 Gy.

There were four threshold dose values reported among the 6 published studies (7, 8, 21, 23), whose average value agrees with the simulation results, 0.043 Gy vs. 0.045 Gy, respectively. These findings are consistent with the published above-zero threshold doses being caused by statistical bias in the threshold dose parameter.

The quadratic parameter whose distribution is shown in Fig. 3b also has a high percentage of false non-zero values, in fact nearly 100%, but the values are equally likely to be negative as positive. Therefore, there is no potential problem for meta-analysis across studies.

Only 1.7% of the above-zero thresholds in Fig. 3a rose to the level of a false positive, so there is no increase in false positives above the expected value of 5%. However, without knowledge of the high expectation rate for dose thresholds with high P values, a reader might see an apparent inconsistency between the reported high P values and the appearance of above-zero dose thresholds in 3 out of 4 cohort studies that tested for threshold nonlinearity (7, 8, 21, 23).

Requiring an AIC test to label a threshold model fit as an improvement was another way to discount most of the above-zero dose thresholds found in the simulations. Only 6.3% of dose threshold values in the upper panel were associated with improvements in fit sufficient to satisfy the Aikake information criterion. Across regression methods and data treatment, the percentage of dose threshold values meeting an AIC criterion varied from 6.3 to 23%. There were three above-zero dose thresholds with P values reported in the 6 studies considered, all of which failed to qualify under AIC selection as an improvement in fit. Their P values were too high, lying above the 0.157 boundary that, for nested models, determines whether a fit is an improvement under the AIC criterion (50).

Results Related to Statistical Modeling

2-Slope Spline and Threshold. The 2-slope spline function does not have the standard asymptotic properties (53, 54) on which standard software may rely to quantify error rates of type I, and convergence of a likelihood test with regression sample size is slow (55). The profile likelihood curve was usually discontinuous at 0 breakpoint and always constant thereafter within the first and last regions defined by dose category boundaries. Discontinuous jumps in profile likelihood derivatives, as well as multiple peaks, could appear in other profiled regions.

Although the profile likelihood for the threshold model was continuous at zero breakpoint, and the likelihood was only constant for the last dose region, it too could have discontinuous profile likelihood derivatives and multiple peaks within its boundaries. Note that a two-peak likelihood surface for threshold dose, described as unusual, was reported in the 2017 A-bomb study (8). Example profile likelihood curves for 2-slope spline and dose threshold are shown in  Supplementary figs. S-6 (rare-199-04-03_s01.pdf) and  S-7 (rare-199-04-03_s01.pdf) ( https://doi.org/10.1667/RADE-21-00219.1.S2).

This nonstandard profile likelihood behavior generally prevented off-the-shelf software from finding the proper parameter values and critical points for threshold dose response and 2-slope splines. The software would not always find the maximum of the maxima within segments, which is why segmented searches were made between dose category values.

Although the 2-slope spline breakpoint showed the same peaking structure above zero breakpoint as found in the threshold model, the lower slope was not constrained in any way. It split in simulations evenly between positive and negative values (results not shown), showing no preference for either sublinear or supralinear dose response. Thus, false findings across studies would also split between sublinear and supralinear shapes and would not present a problem for meta-analysis.

The frequency of false positives for dose threshold, sublinearity and supralinearity with 2-slope spline were <10% for Poisson regression of absolute risk data, but accurate results required careful optimization of the parameter searches. Calculated rates were high if a naïve search for optimum dose threshold value were made using derivative methods for obtaining confidence intervals ( Supplementary table S-7b (rare-199-04-03_s02.docx);  https://doi.org/10.1667/RADE21-00219.1.S2). As expected, differences in false positive rates for 2-slope spline fits between regression methods decreased as regression sample sizes were expanded ( Supplementary table S-10 (rare-199-04-03_s02.docx) and  Supplementary figs. S-8 (rare-199-04-03_s01.pdf) and  S-9 (rare-199-04-03_s01.pdf);  https://doi.org/10.1667/RADE-21-00219.1.S2).

The false positive rates did not converge to 5%. This was initially surprising to us but was consistent with the literature in cases when the null hypothesis (zero breakpoint) lies on the boundary of a constrained parameter space (53, 54). Quite a bit is known in such situations, “... but knowledge is scattered across the literature and considerably less well known among practitioners” (56). A mixture distribution of chi-square statistics is expected that can depend on the configuration of the datapoints (53). As a result, the standard numerical values for setting critical points for false positives are not valid in these special cases. This fine point made little practical difference in the simulations: false positive rates for threshold did not increase above 10% for standard Poisson regression and barely increased over 10% for the case of weighted linear regression of absolute data (11% for threshold dose in the Techa River incidence case).

Curvature

Dose response functions with breakpoints were not the only examples of unusual profile likelihood curves. Curvature had a profile likelihood shape that almost always contained both a valley and a peak, an example of which is shown in  Supplementary fig. S-10 (rare-199-04-03_s01.pdf) ( https://doi.org/10.1667/RADE-21-00219.1.S2).

Individualizing Data

As for loess shape fitting to individualized survival data, it was possible to carry out the full 5,000 replications. Sublinear and supralinear false positive rates were 10% and 9.3%, respectively. This is 4 to 5 times higher than the results for loess fits to grouped cancer risks at Techa River. However, such a comparison is not between equivalents. One endpoint is survival time, and the other endpoint is accumulated mortality counts.

What was most notable about the loess fitting to individualized survival times for the Techa River mortality dataset was the low frequency of hormesis, with hormesis reference level equal to the average survival time at zero dose. The value of 0.6 per thousand was 300 times smaller than the corresponding hormesis frequency of 18% accumulated using loess fitting to grouped risk data. This finding is not surprising in light of the recommendation to use a loess algorithm only when the number of datapoints is large (36), which was not the case with the analysis of grouped counts.

Graphs corresponding to fits to a dose threshold model for a single replication of individualized Techa River mortality data are shown in  Supplementary fig. S-11 (rare-199-04-03_s01.pdf) ( https://doi.org/10.1667/RADE-21-00219.1.S2). In the individualized 2-slope spline case (not shown), bumps and oscillations in the profile likelihood for breakpoint were more common than in the corresponding dose threshold plots.  Supplementary table S-11 (rare-199-04-03_s02.docx) shows that false positive rates found when simulated individuals were analyzed with a Cox model were similar to the false positive rates obtained when groups of these individuals were analyzed as Poisson counts.

Sensitivity Analyses

Sensitivity analyses included varying the dose response slope used in linearization and the smoothing parameter in loess fitting. None of these variations led to any changes in conclusions. Nor did variations in A-bomb parameters, specifically, the choice of reference level, the upper dose range, and the number of persons considered. Details of sensitivity analyses, including those mentioned earlier in the text, can be found in  Supplementary text S-5 (rare-199-04-03_s01.pdf) ( https://doi.org/10.1667/RADE-21-00219.1.S2).

Regulatory Default

The choice of a linear dose response function for cancer causation as the null model is a policy choice by regulators, which in the past has coincided with long-standing recommendations by almost all national and international regulatory and scientific organizations (2–4). We have assumed in this paper that when considering changes to the status quo, regulators of radiation exposure would continue to give the LNT privileged status as the default dose response model. Were a threshold model to be the regulatory norm for a substance, it would have the privileged status and it would be the linear model that would be expected to meet a “clearly superior’ ranking before adoption (57). False positives would still be an issue, but it would be false positives for linearity or other nonconforming responses that would need to be considered.

From the scientific perspective, the linear model is simple and transparent, satisfying Occam's razor in adding only one free parameter to the constant term. It and the pure quadratic function have the least number of parameters of standard models, which makes them experimentally attractive as the models to use first when count numbers are limited.

The pure quadratic model is rarely considered in current radiation cancer epidemiology. Pure linearity, however, is the usual starting point, as indicated by its use as the default in all the epidemiologic studies considered for this paper. Yet, there is potential merit in analyzing dose response in ways that do not privilege the linear LNT model. As we have shown in this paper (Table 3), these methods can also exhibit substantial false indications of nonlinearity within a study, which may be of interest to regulators.

Multi-model fitting, like other data smoothing techniques, has an interesting potential advantage in regulatory debates about dose response. Every multi-model fit is likely to be nonlinear to some degree, which could change the conversation from, “Is the fit nonlinear?” to “How much nonlinearity is there, and does it matter?” Although not explored here, one complexity of multi-model averaging made evident by the results of this paper is that the threshold dose bias will be averaged into the multi-model shape to some extent, if a threshold dose response is one of the models included in the multi-model average.

Threshold Dose Bias

The simulation threshold results dramatize the fact that a bias exists for threshold dose. Intuitively, after seeing the simulation results in Fig. 3a, one might guess that half a bell curve would fall above zero dose and the other half would be clamped at zero, but an analytic proof is not obvious. Will more data and more dose categories reduce the numbers? The results from analysis of the toy threshold-dose-response function suggest not. Graphical analysis of the toy function indicates that, regardless of the magnitude of the spacing between dose categories, the appearance of above-zero dose thresholds in fits to data will be common. There exists a location for a threshold dose between the first two datapoints that will always reduce the residuals of the fit to the data, if two common conditions are met (Materials and Methods section and  Supplementary text S-4 (rare-199-04-03_s01.pdf);  https://doi.org/10.1667/RADE-21-00219.1.S2). This predicted insensitivity to the distance between dose categories was supported by the high frequency of dose thresholds found in the simulation of 17,500 individuals in the Techa River mortality dataset, where the spacing between individual dose categories was very small.

Hormesis

Hormesis is a biphasic dose response in which effects at low doses are opposite to those at high doses. In the present context, this suggests that low-dose radiation exposures are beneficial. Any model that predicts protective effects in a low-dose range would be difficult to fit into a regulatory regime. One reason is that there are concerns about the generalizability of hormetic responses (58). A major obstacle in regulatory application would be that the actual breakpoint dose for a hypothetical protective effect would be uncertain, and if one were identified, it could well vary among subgroups owing to genetic variation in susceptibility (59).

Multiple Comparisons

The results for aggregated false positives per study (Tables 2 and 3) were restricted to primary tests of nonlinearity and thus underestimate the possible instances of multiple comparisons because they do not account for secondary tests. The 2017 A-bomb study by Grant et al. (8) assessed a large number of dose ranges, historical dosimetry schemes, and other sensitivity-test subgroupings. We counted 8 confidence intervals for quadraticity, 10 P values for curvature in the main text, and almost 100 P values for curvature in the appendices. These indicators are not independent. While it is scientifically necessary to analyze all such fits to data to avoid missing any nonlinearities (60), the presentation of so many results without any guidance about multiple comparisons might lead some readers astray.

As several thoughtful reviews and commentaries suggest, even the research community is not immune to misinterpretations of statistical significance (35, 61, 62). This may be an appropriate moment to pay attention in regulatory analysis to the complications of multiple comparisons and the large numbers of statistical inferences that are made, within and across studies.

Should adjustments be made for multiple comparisons? Adjusting significance to account for multiple comparisons may lead to an increase in the number of false negatives, which could lead to important research findings being overlooked (60). From a decision-maker's perspective, however, adjusting for multiple comparisons makes statistical sense, but not necessarily with simple approaches such as a Bonferroni correction. It would be better to consider the costs of the false positives and false negatives (63, 64), which might be measured in social, health, and/or monetary units.

A Bonferroni correction may still be useful as a sensitivity calculation, because it gives an extreme, where the cost of a false negative is assumed negligible compared to the cost of a false positive. When a great number of tests are carried out, a false discovery rate calculation will be less conservative than a Bonferroni correction (65), but may still be incomplete from a decision-support perspective without some idea of cost functions for false positives and negatives. Overall, given the debate over adjustment and its complexity, there is no simple recommendation that can be made.

Simulation Frequency Results Compared to Study Results

As for false positives, there was only one study among the six that reported a positive for its primary analysis using all of the data without subgrouping. That study result was a curvature test in the 2017 A-bomb study with a P value of 0.03. To compare the result to simulation results is straightforward for this primary analysis, as discussed in the Results section, which predicted a 50% chance of finding a false positive among the 20 tests for nonlinearity in the six studies. Thus, both a true and false positive are compatible with the A-bomb primary finding, making the finding uninformative.

However, there is additional evidence to consider beyond the crossing of the bright-line P value of 5%. In the case of the curvature finding for males in the 2017 A-bomb study, the magnitude of the curvature was 1.3 and the P value was 0.002. Even multiplying 0.002 by a conservative Bonferroni correction of 40 (Materials and Methods) would leave the chances of a false finding to be low at 0.08, which argues against the curvature finding for males being explainable by multiple comparisons, although still possible. In simulations, a false positive for curvature for males equal to or greater than 1.3 occurred only with a frequency of 0.001, consistent with the low P value of 0.002.

There are findings in other A-bomb studies that also might be considered by a review committee. For instance, the 2012 A-bomb study reported a positive for upward curvature (P = 0.02), not as a primary test, but as a test within a 6-category subgroup consisting of a 2-category dose range for three time periods (their table 7). Whether or not quadraticity and dose threshold were also tested in these 6 subgroupings was not specified, and we did not try to simulate this situation. The 2000 A-bomb study (6), whose data we did not linearize, presented a smoothed loess fit that trends in the opposite direction for nonlinearity, appearing graphically to be a positive finding for supralinearity, not sublinearity. Uncertainty in the neutron component of dose may also be relevant to interpretations of the A-bomb curvature finding. Increasing the neutron RBE to the higher end of its uncertainty range made the curvature in the A-bomb data negative, as has been reported (66).

Unknown Error Rates of Type I

It can be difficult to argue that one fitting method or data treatment is superior to the others. Poisson regression seems natural for radiation cohort studies, but errors other than count variations are assumed to be negligible. If an analyst felt that error types other than count fluctuations were important, such as errors in dosimetry and/or cancer mortality classification (67), it might be reasonable to prefer weighted linear regression, which accounts for residual errors using the data directly. On the other hand, it does seem from the results in this paper (Fig. 2) that a relative risk data model can produce very high error rates of Type I when both the number of data categories and the number of counts in the reference level are small. An alternative way of getting relative risks in such cases might be considered, for instance by reparametrizing the regression equations before making maximum likelihood estimates or by relativizing absolute risks after fitting, not before.

Linearization in Case/Control Studies

The linearization method presented here for cohort studies is straightforward. For case-control studies, which do not have person-years to adjust, count data would have to be adjusted. A reasonable approach to linearization in these kinds of studies would be to adjust both case and control counts, so that at each dose category the standard deviation of the log odds ratio is kept constant upon linearization ( Supplementary text S-6 (rare-199-04-03_s01.pdf);  https://doi.org/10.1667/RADE-2100219.1.S2).

Methodological Limitations

Because covariates such as age, sex, and city of bombing were not included in the regressions, the results strictly apply only to univariate regressions and could change somewhat with multivariate analysis. Thus, comparison with univariate study results would be the most direct use of the method described here. It has been assumed that the datasets are independent. Yet, the 2007 and 2017 A-bomb studies have overlapping cancer incidence data, although the methodology, including dosimetry and choice of covariates, has changed between publication dates, introducing some additional effective randomness.

In some studies, uncertainty in dose values has been taken into account in novel ways (68). Although variance in dose values could have been introduced when adding Poisson variation to linearized count data, we did not analyze this possibility. The methods considered here included Poisson regression, weighted linear regression, regression of relative risk data, AIC selection, and multi-model fitting. Other analysis methods that were not covered may be of interest to some readers (32, 46, 47, 68, 69), including one study that used simulation of data replicates to check confidence limits (47).

False Positive Rates for a Single Test of Nonlinearity

According to simulations of the six linearized cohort studies, there were two situations where false positive rates for null-hypothesis tests of nonlinear cancer dose response exceeded 10%. The first situation was software related. Some common regression software could not correctly analyze models like dose threshold that had discontinuous derivatives and could have multiple peaks for parameter values on the likelihood surface. The excess rates were limited to dose threshold and 2-slope spline models, but the excesses appeared in all datasets.

The second problematic situation was caused by unusually low counts in the reference categories of two datasets, specifically the Techa River incidence and mortality datasets. The excess rates only occurred if nonstandard regression methods or data treatment were used. High false positive rates did not appear for standard Poisson regression of absolute risk data, even with the low counts in the reference category, except in a single instance of the three definitions of hormesis assessed.

Linearization combined with simulation proved to be an informative diagnostic tool for identifying problematic methods for detecting nonlinearities in dose response. The analysis carried out in this paper has largely assumed a regulatory default that is linear without threshold, but the methods introduced here can be adapted to use other types of dose response as the default, should a different one be selected in the future by regulators.

Multiple Tests in a Single Study

Even when individual tests of nonlinearity have low false positive rates, conducting multiple tests for nonlinearity in a study can present a multiple comparisons problem. As many as five tests per study for nonlinear dose-response functions can be identified in epidemiological studies of radiation and cancer (7–9, 21, 23). Consequently, the expected frequency of false positives in a study can be substantial, about 19% per study for the three A-bomb datasets according to simulations of linearized data, using Poisson regression of absolute risk data below 0.6 Gy. The value was 15%, when averaged over all six studies (Fig. 2).

The false positive rates quoted above were obtained for null-hypothesis testing. Considering approaches other than null-hypothesis testing, such as AIC model selection or multi-model fitting, did not solve the problem of false nonlinearities turning up in simulations (Table 3).

Multiple Tests in Multiple Studies

When datasets were considered as a group, the chance of at least 1 false positive occurring in the 20 primary test results available in the 6 published studies was estimated to be 50%. Thus, the 1 positive found in the actual six studies using all the data (primary finding) with a P of 0.03 could easily be a false positive. In subgroup analysis by sex, there was strength of evidence to evaluate in the case of the finding of curvature for males in the 2017 A-bomb study. The magnitude of the curvature parameter was greater than 1, and the P value was 0.002. Binomial calculations combined with simulations and a Bonferroni correction suggested that it was unlikely that such a result was explainable by multiple comparisons.

Above-Zero Dose Threshold Findings

In the six studies considered, four of them fitted their data to a dose threshold model. Of these, three reported above-zero dose threshold values averaging to about 0.05 Gy. Although none rose to the level of a positive, their high frequency was unexpected. Simulation provided a simple explanation. Analysis of linearized data showed that a high frequency was to be expected simply on the basis of random fluctuations in counts coupled with a statistical bias in the threshold dose parameter.

Without accounting for threshold dose bias or correcting for the extra model parameter, meta-analyses of above-zero dose thresholds, formal or informal, will be invalid. AIC methods provided a correction method by penalizing threshold fits for the extra parameter introduced. In simulations, this dramatically reduced the frequency of above-zero threshold fits that could be considered an improvement over the linear null. Furthermore, AIC methods could be applied to published P values and were sufficient to declare the published cases poor fits to the data, with no improvement over a linear fit.

Implications for Regulatory Analysis

A single isolated result from a study is sometimes presented to policy makers and the public as “statistically significant,” and may therefore be judged by some as convincing evidence of nonlinearity at low dose levels and falsification of the linear regulatory default model. These two words, statistically significant, have developed a magic in the wider culture that can overwhelm nuanced conversation and affect public policy. However, without screening for false positivity, isolated findings should not be considered valid evidence against the regulatory default. Binomial calculations and simulation as part of regulatory analysis may help in judging if a result might be a false positive.

Implications for Radiation Research

When using nonstandard and novel methods of analysis on cohort datasets or subsets of them, or when there is a small number of count numbers or dose categories, it could be helpful for authors to quantify false positive rates. The linearization and simulation techniques presented here may be useful for this purpose.

Regulatory policy is not set by researchers, but their language and the visuals they present can affect it. For instance, epidemiologic studies can be cited and debated by those trying to influence regulatory policy, legislation and public opinion. In discussing their finding of nonlinearity for male cancer dose response, Grant et al. (8) were sensitive to the authors' role in shaping regulation, warning against trying at this time to “confidently guide the development of modified radiation protection policies.” Their reasoning went beyond study uncertainties to include the contrasting results found in other studies. This is recognition of the multiple comparisons problem and an indirect warning of the possibility of false positives.

The simulation results in our paper strengthen the idea that explicit discussion of the potential for false positives (at least qualitatively), which is common practice in many fields, could be a useful addition to a study in radiation cancer epidemiology if the reporting or discussion of nonlinearity is deemed relevant to regulatory analysis. The recognition in papers of the possibility of false positives could reinforce the need to look at the family of studies before reaching decisions, which is what expert review committees do.

Another option that could be considered is to avoid the shorthand words, “statistically significant.” Authors and commenters might instead consider language that is less dependent on the drawing of bright lines, for instance, by discussing degrees of compatibility (70).

Attempting to correct for multiple comparisons is not necessarily a better solution for researchers, given the controversy over its desirability in the research context (60, 64). Furthermore, adjusting P values or confidence intervals in one study does not account for the comparisons across studies that can raise to high levels the chances of finding a false positive. When results are of regulatory interest, it can be constructive to point out the importance of considering the multiple comparisons problem and the possibility of associated false positives