Methods to Derive Uncertainty Intervals for Lifetime Risks for Lung Cancer Related to Occupational Radon Exposure.

INTRODUCTION
lifetime risks describe the probability of developing or dying from a specific disease (here: lung cancer death related to radon exposure) in the course of one’s lifetime and play an important role in the epidemiological approach for radon dose conversion of exposure in Working Level Months (WLM) to effective dose in millisieverts (mSv) (ICRP 1993, 2007, 2010) and radiation detriment, a concept used by the ICRP to quantify the overall harm caused by exposure to ionizing radiation (ICRP 2022). A lifetime risk estimate depends on several components, each imposing possible errors or uncertainties on the result. Lifetime lung cancer risk estimates related to radon exposure depend on (a) the exposure scenario, (b) baseline mortality rates for all causes of death and lung cancer, and (c) complex risk models describing the shape of the exposure-response relationship between radon exposure and lung cancer mortality. For occupational radon exposure, these risk models are derived from uranium miners cohorts. We focus on the latter two components (b) and (c) in the upcoming analysis, as previous work (Sommer et al. 2025) demonstrated their major influence on lifetime risk estimates.
Globally, millions of workers are potentially occupationally exposed to radon, particularly in uranium mining regions such as Kazakhstan, Canada, and Australia, as well as in industries such as non-uranium mining, caves, spas, and waterworks in countries like Germany (Beck and Ettenhuber 2006; Fan et al. 2016; Daniels and Schubauer-Berigan 2017; NEA 2023). Exposure to radon and radon progeny is a recognized leading cause of lung cancer. This association has been confirmed in uranium miners and residential radon studies (NRC 1999; Kreuzer et al. 2018; UNSCEAR 2021). Both uranium miners and residential radon studies support a linear relationship between radon exposure and lung cancer risk within the considered exposure ranges (Darby et al. 2005; UNSCEAR 2021), which is in case of uranium miners studies additionally influenced by factors like age, time since exposure, and exposure rate (Richardson et al. 2022; Kelly-Reif et al. 2023). The intricate nature of these models can make comparisons across different cohorts challenging. Lifetime risk estimates offer a valuable tool for comparison and interpretation of risk models and enable clearer public risk communication. However, current literature often lacks uncertainty intervals for lifetime risk estimates. Lifetime Excess Absolute Risk (LEAR) estimates derived from various international miners studies range from 2.5 ×10−4 to 9.2 ×10−4 per WLM as reported by the large pooled uranium miners study PUMA (Kreuzer et al. 2024). The present study employs advanced statistical methods to quantify lifetime risk uncertainty with 95% uncertainty intervals.
Uncertainty intervals together with point estimates provide a more complete and more nuanced picture, allowing for informed decision-making, meaningful comparisons, and transparent communication. Uncertainty assessment for lifetime risks is a complex endeavor because final lifetime risk estimates are composite quantities depending on multiple, independently derived results.
The literature addressing lifetime risk estimation uncertainty for lung cancer related to radon exposure based on uranium miners cohorts is limited and generally does not prioritize uncertainty quantification (Chiu and Brattin 2003; Tomasek 2020). Existing studies typically employ Monte Carlo simulation techniques, incorporating various distributional assumptions for calculation components. For example, Tomasek (2020) investigated the lifetime excess absolute risk (LEAR) and assumed a multivariate normal distribution for risk model parameters, while in the report of the United States Environmental Protection Agency (US EPA), an additional complex distribution for average residential radon exposure is assumed (Chiu and Brattin 2003). A comparison of uncertainties across different lifetime risk measures could not be found in the literature. Uncertainties in risk measurement estimates similar to lifetime lung cancer risks related to radon exposure, like “attributable cases,” are better understood. In the Biological Effects of Ionizing Radiation (BEIR) VI report (NRC 1999), Monte Carlo simulations were used to comprehensively quantify uncertainties for “attributable cases,” similar to methods for lifetime lung cancer risks in Tomasek (2020).
Several software tools have been developed for calculating lifetime cancer risks and associated uncertainties. These tools, however, primarily rely on risk models derived from the Atomic Bomb Survivors Life-Span Study (LSS) (Preston et al. 2007), which involves acute external exposure to other radiation types, fundamentally different from the chronic internal exposure to radon progeny in occupational uranium miners studies. Notable among these tools are “CONFIDENCE” (Walsh et al. 2020), “RadRAT” (Berrington de Gonzalez et al. 2012), and “LARisk” (Lee et al. 2022). “CONFIDENCE” is a European software designed for cancer risk assessment post-radiation exposure from nuclear accidents. “RadRAT”, a free tool based on the BEIR VII report (National Research Council 2006), estimates lifetime cancer risks for the US population based on a user-specific exposure scenario. “LARisk” extends “RadRAT” by adding flexibility, such as modified baseline incidence rates, to tailor risk calculations for specific populations. While “CONFIDENCE” and “RadRAT” rely on Monte Carlo simulations and sampling from probability distributions for uncertainty assessment, “LARisk” employs parametric bootstrap methods.
Despite these advanced tools, there is a notable gap in the literature regarding the uncertainty quantification of lifetime risks for lung cancer specifically related to occupational radon exposure. Addressing this, our study places special emphasis on refining the understanding of these uncertainties, using risk models derived from uranium miners cohorts. We specifically focus on the LEAR measure while also examining alternative excess lifetime risk measures: Risk of Exposure Induced Death (REID), Excess Lifetime Risk (ELR) (Thomas et al. 1992; Kellerer et al. 2001; Ulanowski et al. 2019), and Radiation-Attributed Decrease of Survival (RADS) (Ulanowski et al. 2019). We elaborate on sources of uncertainty for lifetime risks, discuss existing techniques, and introduce advanced methods to quantify these uncertainties by calculating 95% uncertainty intervals.
To derive such intervals, both frequentist methods derived from likelihood theory and Bayesian techniques, chosen for their flexibility, are employed to quantify uncertainty in risk model parameter estimates. These methods are applied to risk models derived from the German Wismut uranium miners cohort study as a practical example. Mortality rate uncertainty is assessed based on information from the WHO mortality database (WHO 2022). All uncertainty assessment methods are realized through Monte Carlo simulations.
This study builds upon previous research on lifetime risk sensitivities (Sommer et al. 2025) and introduces a novel and comprehensive statistical framework for deriving uncertainty intervals for lifetime risk estimates. In addition to advancing these methodologies, the study also provides a thorough analysis and discussion of the derived intervals, thereby contributing to a more stable basis for radon protection strategies.

MATERIALS AND METHODS

Lifetime risk definition and calculation
For a lifetime risk measure of choice, we employ the Lifetime Excess Absolute Risk (LEAR). This measure is often used (Vaeth and Pierce 1990; Kellerer et al. 2001; Tomasek et al. 2008a; Tomasek 2020, Kreuzer et al. 2023) and is defined as:
The lifetime risk of dying from a specific disease (here: lung cancer) under exposure (here: radon and radon progeny) is and is the corresponding baseline lifetime risk without exposure. r0(t) is the baseline lung cancer mortality rate, and rE(t) is the lung cancer mortality rate under exposure at age t. S(t) = ℙ(T ≥ t) is the survival function and describes the probability to attain age t with T ≥ 0 the unknown random retention time until death. The survival function is modeled as with baseline mortality rates q0(u) at age u for overall death (all-cause mortality rates). ERR(t; Θ) specifies the excess relative risk at age t with unknown parameter set Θ. Overall, we assume the following risk projection model:
The exact structure and complexity of the ERR(t) = ERR(t; Θ) term depends, in addition to age t, on the chosen risk model with parameters Θ and its included effect-modifying variables. The LEAR from eqn (1) is estimated and finally calculated with the approximation
For lifetime risk calculations, parameter estimates are plugged into the corresponding risk model structure. approximates the true survival function S(t) and is based on the Nelson-Aalen estimator of the cumulative hazard rate (Nelson 1972). The maximum age tmax was set to tmax = 94 for comparability with previous studies on lifetime risks (Tomasek et al. 2008a; Kreuzer et al. 2023).
A lifetime risk estimate further depends on the chosen risk model shaping ERR (t; Θ), the mortality rates r0(t), q0(t) for all ages 0 ≤ t ≤ tmax and the assumed radon exposure in WLM for all ages 0 ≤ t ≤ tmax (exposure scenario). In line with all considered risk models, the latency time L in risk models, i.e., the minimal time between age at exposure and age at lung cancer risk amplification, was chosen as L = 5 y.
Here, a lifetime risk estimate for a given risk model without uncertainty quantification is always calculated using the risk model specific parameter estimates, with an exposure scenario of 2 WLM y−1 from age 18–64 y (94 WLM total cumulative exposure over lifetime) to represent an occupational scenario, in line with the literature (Tomasek et al. 2008a; Tomasek 2020; Kreuzer et al. 2023) and fixed mortality rates , for all ages 0 ≤ t ≤ 94 from the sex-averaged Euro-American-Asian mixed reference population as presented in ICRP Publication 103 (ICRP 2007). Such a risk model specific lifetime risk estimate does not reflect uncertainties and is referred to as “reference estimate (ICRP 103)” or simply “ref. estimate” in this paper. Further, for enhanced readability, lifetime excess risk estimates such as the LEAR are presented as percentages (LEAR in %, i.e., LEAR ×100).

Data and risk models
To estimate the risk models associated with the German uranium miners cohort, this analysis employs cohort data, including cumulative radon exposure, age, and other variables of individual miners. These cohort data provide the basis for modeling and assessing the developed parametric (continuous) and categorical risk models, particularly for the quantification of uncertainties. All considered risk models include unknown parameters (indicated by Greek letters), which are estimated using maximum-likelihood (ML) methods. Access to the German uranium miners cohort data allows for refitting the subsequent described risk models and thereby for an uncertainty assessment of the resulting lifetime risks.
The following parametric (continuous) risk models are considered:
with cumulative radon exposure W(t) at age t in WLM and continuous effect-modifying variables age at median exposure AME(t) at age t in years, time since median exposure TME(t) at age t in years, and binary variables erj for j = 1, …, 6 for six categories of exposure rate at age t in units of working level (WL). Further, categorical BEIR VI exposure-age-concentration models (National Research Council 2006) are considered:
with θ(1) = (θ5 − 14, θ15 − 24, θ25 − 34, θ35+), θ(2) = (θ5 − 14, θ15 − 24, θ25+) and where W5 − 14, W15 − 24, W25 − 34, W25+, W35+ is the cumulative radon exposure in WLM in windows 5–14, 15–24, 25–34, 25+ or 35+ years ago and ϕage and γrate are factors for attained age in years and exposure rate in WL, respectively.
All presented risk models are derived from the German Wismut cohort of uranium miners with follow-up up to 2018 and have been initially introduced in Kreuzer et al. (2023). These models are used in this study as an application example of the subsequent introduced methods to quantify risk model parameter uncertainty. The models eqn (5) “Parametric 1960+ sub-cohort,” eqn (6) “Simple linear 1960+ sub-cohort,” eqn (8) “BEIR VI 1960+ sub-cohort” are derived from the Wismut sub-cohort with miners hired in 1960 or later (1960+ sub-cohort), and models eqn (4) “Parametric full cohort” and eqn (7) “BEIR VI full cohort” are derived from the corresponding full Wismut cohort. The derived parameter estimates are shown partly in the Supplemental Digital Content (Supplemental Digital Content; http://links.lww.com/HP/A317) or can be inspected in the source paper (Kreuzer et al. 2023). The simple linear model eqn (6) is included in this analysis for comparability. All other models were chosen here because in earlier work they showed the best fits for the corresponding cohort (Kreuzer et al. 2023). The 1960+ sub-cohort models allow investigation of lung cancer risk under contemporary conditions of low to moderate radon exposure, are characterized by uniform exposure rates, provide higher quality exposure data and more detailed information on confounding variables, and thereby better reflect current occupational scenarios (Laurier et al. 2020; Kreuzer et al. 2023). The chosen risk models (categorical/BEIR VI structure, parametric/continuous with effect modifying variables, and simple linear risk model) offer a diverse range of risk models. Here, the terms “categorical” and “parametric/continuous” refer to the categorical or continuous nature of the effect-modifying variables.
All considered risk model parameters are estimated with maximum likelihood (ML) methods based on internal Poisson regression methods applied to grouped cohort data, with estimation of baseline rates based on cohort data (internal) and not on external data. The corresponding likelihood function is based on the assumption that the number of lung cancer deaths Ci in cell i with i = 1, …, n are Poisson-distributed via
with offset for person-years at risk PYi, excess relative risk ERRi(Θ) and unknown parameter vector Θ. The specific shape of ERRi(Θ) depends on the prescribed risk model structure. The baseline risk predictor describes the stratified baseline with K levels, where is a categorical variable with K levels. Each level represents a unique combination of conditions or classifications that are relevant variables in assessing the baseline risk. Here, the levels correspond to different groups categorized by age, calendar year, and duration of employment. Setting Δ = (δ1, δ2, …, δK), the full likelihood model consists of the parameters Ω = (Δ, Θ).

Uncertainty assessment

Uncertainty intervals for lifetime risks
Uncertainties for a parameter of interest are quantified by deriving 95% uncertainty intervals, which is a range of values that is calculated to cover the true unknown parameter value with 95% certainty. It consists of a lower and an upper bound, enclosed in parentheses. The precise interpretation of the uncertainty interval depends on the underlying statistical inference system (frequentist or Bayesian). Lifetime risks are composite quantities depending on multiple, independently derived results. Hence, we rely on sampling techniques to derive lifetime risk uncertainty intervals. The derived uncertainty intervals for lifetime risks like the LEAR reflect the expectable range of potential values that arise from the inherent variability in each calculation component and cannot be interpreted as classical confidence intervals.
Here, uncertainties are determined by quantifying the variability in statistical estimates that results from drawing a sample from the entire population (sampling uncertainty). Risk model parameter uncertainties and mortality rate uncertainties are considered. Depending on the investigated lifetime risk calculation component, previously fixed values for calculation components are replaced by random variables. The uncertainty quantification is carried out with Monte Carlo simulations. Here, we focus on the excess lifetime risk measure LEAR. N = 100,000 samples from the underlying assumed probability distribution are drawn independently and a LEAR is calculated for each sample resulting in N independent LEAR estimate samples. The two-sided 95% uncertainty interval is the span of observed LEAR samples by disregarding the 2.5% lowest and 2.5% highest samples. All upcoming methods use this approach to derive uncertainty intervals unless explicitly stated otherwise but differ in the calculation component analyzed and the assumed probability distribution. Corresponding probability density functions for the LEAR distribution are derived from the histogram of LEAR samples with a smoothing kernel density estimate. Note that for lifetime risks, we additionally present the relative span of the uncertainty interval in brackets, calculated as the interval span divided by the reference estimate (“relative uncertainty span”). This enables easier comparison across various estimates, making it clearer when an uncertainty interval is relatively large or small, especially when compared to another interval with different absolute length.

Risk models‐ANA approach
The idea behind the Approximate normality assumption (ANA) approach is the following: Risk model parameter estimates are obtained by fitting the ERR(t; Θ) model to miners cohort data with ML and Poisson regression methods assuming Poisson-distributed numbers of lung cancer cases (eqn 9). The ML method provides statistically efficient estimates, meaning they tend to be close to the true values on average given a sufficiently large sample size. The estimated parameters are subject to sampling uncertainty. By statistical theory, under mild regularity conditions met in our case (Amemiya 1985), the parameter estimator in a given risk model is asymptotically (for an infinitely large cohort size) normally distributed with expectation equal to the unknown true parameter value Θ and covariance matrix equal to the inverse Fisher information (Held and Bove 2021; Pawitan 2013). The true parameter Θ is approximated by the cohort-specific maximum likelihood estimate (MLE) . The inverse Fisher information can be approximated by an estimate of the parameter covariance matrix.
The ANA approach follows the frequentist inference where the true risk model parameter Θ is treated as a fixed but unknown value. The parameter estimator Θˆ is considered as a random variable subject to variability depending on the specific sample used in the estimation process. The estimate is a realization of by applying specific sample data. Following this framework, risk model parameter uncertainty is quantified by assuming a multivariate-normal distribution on the parameter estimator ,
with cohort-specific ML estimate and covariance matrix estimate . By generating a large number of samples from this approximate distribution and calculating the corresponding lifetime risks for each sample, the distribution of lifetime risk estimates reflects the inherited parameter sampling uncertainty. This approach is called the “Approximate normality assumption” (ANA) approach. In this study, we employ Epicure software (Preston et al. 1993) to obtain parameter estimates and their associated covariance matrix estimate . Note that this method implicitly incorporates knowledge from the estimation of baseline parameters Δ, as the covariance estimate is adjusted accordingly (Preston et al. 1993). Importantly, after this theoretical introduction, the manuscript simplifies notation by referring to estimates using without further distinguishing the estimator notation. Note that the ANA approach requires access to the ML estimate and covariance matrix estimate , but no access to underlying cohort data since a re-estimation of parameters is not necessary.

Risk models‐Bayesian approach
In contrast to frequentist inference underlying the above ANA approach, in the Bayesian approach (Bayesian inference), the unknown parameter Θ is interpreted as a random variable itself. To apply Bayesian statistics to account for prior information about the risk model parameters we assume the generic Bayesian framework
where P(Ω ∣ X) is the posterior probability density function of observing the parameter Ω given cohort data X. Here Ω = (Δ, Θ), compare eqn (9). P(Ω) is the prior probability density for Ω and L(X) = L(X ∣ Ω) is the likelihood function. Using concepts from Higueras and Howes (2018) allows to derive the marginal posterior distribution P(Θ∣ X) for the risk model parameters of interest Θ analytically. This requires assuming independence between prior distributions for Δ and Θ and a non-informative prior for Δ. The resulting marginal posterior given cohort data X reads
with lung cancer cases in strata k, Sk = ∑i ∣ xi = k Ci for k = 1, …, K and normalizing constant M. Parameter estimates are derived as the values that maximize the posterior distribution (mode), denoted as , analogous to the ANA approach. Note that the Bayesian approach requires a re-estimation of parameters and therefore access to the original cohort data in contrast to the ANA approach.
This approach is applied to the 1960+ sub-cohort models with eqns (5) and (6). Identical risk model structures were used in Tomasek et al. (2008b), which we employ as prior information. Due to increasing computational complexity, it was not feasible to apply this approach to full cohort risk models eqns (4), (7), and sub-cohort model eqn (8), which involve larger cohort data size and/or more parameters. Non-informative, uniform prior distributions for Θ are applied to obtain the true marginal likelihood. Otherwise, the prior information about β in the simple linear risk model eqn (6) is modeled with a gamma distribution with from Tomasek et al. (2008b). The prior information about β, α and ε for model (5) is modeled as , and normal distributions with from Tomasek et al. (2008b). Prior certainty is increased by either raising the gamma shape parameter a for β or reducing the standard deviation σ for α and ɛ. By construction, the modes of the marginal prior distributions align with the corresponding parameter estimate from Tomasek et al. (2008b). All three components are assumed to be independent in the prior; i.e., P(Θ) = P(β) · P(α) · P(ε). Importantly, estimates from Tomasek et al. are not assumed as true “prior” knowledge but illustrate integrating diverse cohort information using Bayesian methods (see Discussion).
For the simple linear model eqn (6), Rejection Sampling (Held and Bove 2021) was applied to obtain N = 100,000 samples of the posterior distribution of β with uniform proposal distribution U(0,0.04). For the more complex risk model in eqn (5), Markov Chain Monte Carlo (MCMC) techniques via the Metropolis-Hastings algorithm were applied (Robert and Casella 2005). The approximate multivariate normal distribution eqn (12) was chosen as the proposal distribution. The log acceptance ratio for a proposal sample for Θ ∼ P(Θ| X) was calculated as the difference of the log marginal posterior evaluated at the proposal sample and the current sample. The initial proposal sample for Θ was chosen in proximity of the cohort-specific ML estimate . N = 100,000 samples of the posterior distribution were used for uncertainty assessment after generating 110,000 samples and discarding the first 10,000 samples to account for a burn-in period. Here, the presented 95% LEAR uncertainty intervals are derived by choosing the narrowest interval that covers 95% of the derived LEAR samples, which we refer to as the Highest Posterior Density Interval (HPDI).

Mortality rates
Uncertainties in the mortality rates r0(t), q0(t) are assessed by assuming gamma distributions for all ages t,
with age-dependent shape parameters and rate parameters (SDC Table S6). The parameter estimates are derived from data from the WHO Mortality Database (WHO 2022) with maximum-likelihood (ML) methods. Observations from all available countries from Europe, America, and Asia for females and males from the years 2001, 2006, 2011, 2016, and 2021 are used. The gamma distribution was selected due to its characteristic of being positive, right-skewed, and its strong fit to the observed mortality rates across numerous age groups (SDC Figure S7). The derivation of probability distributions for mortality rates in (10) and (11) from WHO data is independent of the sex-averaged ICRP Euro-American-Asian reference mortality rates and . However, the geographic alignment of the chosen WHO data with the ICRP mortality rates ensures the appropriateness of the data for estimating mortality rate variability. An observation used for fitting models (10) and (11) is characterized as the mortality rate , where d are the number of lung cancer deaths for r0 or all-cause deaths for q0, and n is the mid-year population size. Each observation is uniquely defined by a specific country (out of 153 countries), sex, and calendar year. Only observations with a positive number of individuals at risk n and a positive number of cases d were considered. To obtain uncertainty intervals for lifetime risks, independent samples from the above-described gamma distributions (10) and (11) are drawn (Monte Carlo simulation) analogously to drawing samples in the ANA approach for risk model parameters.

Joint effect of risk model and mortality rate uncertainty
The joint effect of risk model parameter estimate uncertainty and mortality rate uncertainty on LEAR estimates is assessed by simultaneously sampling from multivariate normal distributions for the risk model parameter estimator according to the ANA approach (12) and from gamma distributions for the mortality rates r0(t) for all ages t according to eqn (10), independently. Note that q0(t) variability is here not accounted for, as initial analysis showed a negligible impact (see Results). Analogous to the general Monte Carlo simulation approach, N = 100,000 LEAR samples are calculated from N = 100,000 sets of sampled values for r0(t) for all ages t and . All computations were carried out with the statistical software R (R Development Core Team 2024).

Supplementary and sensitivity analyses
The Supplemental Digital Content (SDC) (http://links.lww.com/HP/A317) contains extensive analyses that provide a more comprehensive understanding of lifetime risk uncertainties. SDC Section A details uncertainties for other lifetime risk measures besides LEAR. Sensitivity analyses (SDC Sections B to D) explore the impact of assuming different probability distributions on lifetime risk uncertainties. Additional insights on risk model parameter uncertainty are shown in SDC Section B.1 and B.3. Furthermore, the Bayesian approach for risk model parameter uncertainty quantification is applied to the simple linear risk model eqn (6) with log-normal priors for β (SDC Section B.2). Analyses for the specific influence of all-cause mortality rate uncertainties are conducted (SDC Section C.1). An alternative Bayesian approach to assess mortality rate uncertainty employing the WHO data is found in the SDC Section C.2. Log-normal distributed mortality rates (SDC Section C.3) and sex-specific uncertainties regarding mortality rate and risk model effects are investigated (SDC Section D). Exposure scenario uncertainty is briefly investigated in the SDC Section E.

RESULTS

Effect of risk model parameter uncertainty on lifetime risks
LEAR estimates derived from full cohort risk models are notably lower than estimates derived with 1960+ sub-cohort models (Kreuzer et al. 2024). This discrepancy is particularly pronounced for the Wismut cohort and primarily comes from the significantly lower estimates for Excess Relative Risk per 100 WLM (ERR/100 WLM) observed in the full cohort (Kreuzer et al. 2023). The variation between these ERR/100 WLM estimates is similarly reflected in the LEAR calculations.

Approximate normality assumption (ANA) approach
Table 1 shows the resulting 95% uncertainty intervals. For risk models fitted on the full cohort, the resulting intervals are comparable, although the model structures differ considerably. The relative uncertainty span (the 95% uncertainty interval span divided by the reference estimate3) is 0.81 and 1.03 for the parametric eqn (4) and the BEIR VI risk model eqn (7), respectively. For the 1960+ sub-cohort models, however, the results vary considerably. Especially risk model eqn (8) implies very wide LEAR uncertainty intervals with a relative uncertainty span of 6.27 and a notable portion of implausible negative LEAR samples. Model eqn (5) implies considerably less uncertainty compared to model eqn (8) with a relative uncertainty span of 1.34, although both models are derived from the 1960+ sub-cohort. Visually (Fig. 1), the empirical distribution of LEAR samples resembles approximately a normal distribution for risk model parameter estimates derived from the full cohort (model eqn 4 and 7). Conversely, the empirical distribution of LEAR samples with parameter estimates derived from the 1960+ sub-cohort inherit considerably heavier tails and a slight right skewness. The 1960+ sub-cohort is smaller and comparably younger with fewer person-years at risk and fewer lung cancer deaths. This results in higher statistical uncertainty, which is reflected in wider LEAR uncertainty intervals.
Note that for the simple linear risk model in eqn (6), we obtain analytically, without sampling, a normal distribution with corresponding 95% uncertainty interval [3.18; 8.22] and a relative uncertainty span of 0.88 for LEAR in %. By definition, this uncertainty interval is proportional to the 95% confidence interval [0.75; 1.93] for 0 × 100 = 1.34 from Kreuzer et al. (2023) by a factor of 4.27 here. The underlying theory is explained in SDC Section B.1.

Bayesian approach
Computational limitations restricted the applicability of this approach to models eqns (5) and (6). For the simple linear model eqn (6), analogous to the ANA approach, LEAR uncertainty is directly proportional to β uncertainty. The posterior mode of P(β|X) with 95% uncertainty intervals (HPDIs) for varying certainty in the prior information is shown in Tables 2 and 3 with corresponding plots in Fig. 2. The resulting HPDI with a uniform prior are comparable to the conventional 95% confidence interval for the estimate . However, the lower bound of the HPDI is notably higher than that of the classical confidence interval for . Further, the relative uncertainty span of the uncertainty interval decreases with increasing certainty in the prior information (increasing a) from 0.94 for a = 2 up to 0.48 for a = 50. The resulting 95% HPDIs for the more complex model eqn (5) incorporating two additional effect-modifying parameters α and ε compared to eqn (6) are shown in Table 2 for the special choice of a standard deviation σ = 0.02 in the prior distributions for α and ε (SDC Section B.3 for more choices for σ and histograms of parameter samples). The corresponding distribution of LEAR samples is shown in Fig. 3 (Plot a). As expected, increasing the prior certainty (increasing shape parameter a) shrinks the uncertainty interval and moves the reference estimate (ICRP 103) of LEAR in % from 6.74 (95% HPDI [2.96; 11.09]) without prior influence toward the prior LEAR in % estimate of 4.30 derived with the Joint Czech+French risk model. For example, a = 50 results in a LEAR in % of 4.99 and a 95% HPDI of [3.14; 7.26]. A very similar effect is observed for decreasing the standard deviation σ for a fixed shape parameter a = 10 (Fig. 3 Plot b). For both models, eqn (5) and eqn (6), the Bayesian approach with a uniform prior demonstrates strong agreement with the ANA approach in terms of point estimates and uncertainty intervals.

Effect of mortality rate uncertainty on lifetime risks
Introductory analyses revealed that all-cause mortality rates q0(t) impose considerably less uncertainty on the LEAR than lung cancer rates r0(t) for all risk models (SDC Section C.1). The empirical distribution of LEAR samples for gamma-distributed q0(t) is considerably narrower compared to the empirical distribution for gamma-distributed r0(t), which is also reflected in the 95% uncertainty intervals. The relative uncertainty span is very similar across all considered risk models with roughly 0.10 for uncertainties in all-cause rates q0(t) and roughly 0.50 for uncertainties in lung cancer rates r0(t) and joint uncertainties r0(t), q0(t) (SDC Table S10). Hence, all-cause mortality rate uncertainty can be neglected reasonably when assessing overall mortality rate uncertainty, as addressed in the subsequent analyses.

Joint effect of mortality rate and risk model parameter uncertainty on lifetime risks
The empirical distribution of resulting LEAR samples, along with the corresponding 95% uncertainty intervals, varies depending on the underlying risk model and its parameters, especially for the 1960+ sub-cohort models (Fig. 1, Table 1, last column). Risk model complexity and cohort size directly influence LEAR uncertainty intervals, as seen in the separate analyses. In contrast, the 95% uncertainty intervals from uncertain r0(t) are overall narrower compared than those from uncertain risk model parameters and consistent across all risk models with a relative uncertainty span of roughly 0.5. The joint effect of uncertain risk model parameters and uncertain lung cancer mortality rates is with a relative uncertainty span of roughly 1, almost similar to the effect with only uncertain risk model parameters. Only the 95% uncertainty interval [−9:91; 23.17] of the joint effect with relative uncertainty span of 5.76 for the BEIR VI 1960+ sub-cohort model eqn (8) remains implausible, although slightly narrower. Overall, accounting for lung cancer mortality rate uncertainty in addition to risk model parameter uncertainty has low impact on the overall LEAR uncertainty interval.

Supplementary results in SDC

Lifetime risk measures
This analysis includes uncertainty intervals for lifetime risk measures such as ELR, REID, and RADS, derived using both ANA and Bayesian techniques (SDC Section A). The relative uncertainty spans across these measures are comparable, with RADS estimates being notably larger than others, suggesting that practical reliance on the LEAR measure is sufficient.

Choice of prior in risk model parameter uncertainty
When using log-normal versus gamma distributed priors for β in the Bayesian approach for the simple linear risk model eqn (6), we observed only minor differences in outcomes (SDC Section B).

Mortality rate uncertainty
Bayesian approaches were applied to mortality rates, using Poisson-distributed numbers of lung cancer deaths, but these provided only minimal additional insights due to the large statistical power of WHO data (SDC Section C). Comparisons between log-normal and gamma distributions for mortality rates yielded similar uncertainty intervals, although log-normal showed slightly larger upper bounds.

Sex-specific lifetime risks
Estimates for sex-specific lifetime risks using WHO mortality data and ICRP reference rates (SDC Section D) show that estimates for males are approximately twice as high as those for females, due to higher baseline mortality in males for all-cause and lung cancer, with similar relative uncertainty spans.

Exposure uncertainty
Analysis of exposure variability (SDC Section E) indicates that the assumptions on annual exposure variability notably affect the distribution of lifetime risk estimates similarly across all risk models and lifetime risk measures. Under conditions of low exposure uncertainty, lifetime risk estimates tend to follow a normal distribution. Since this study investigated the LEAR instead of the LEAR per WLM, the observed variability is to be expected and not contrary to previous findings in (Sommer et al. 2025), which indicated that the LEAR per WLM exhibited only low variation across different exposure scenarios.

DISCUSSION
This work provides a methodological contribution in radiation protection research by successfully deriving uncertainty intervals for radon-induced lifetime lung cancer risk estimates. These intervals are grounded in a sound statistical framework, shifting the approach from solely assessing single lifetime risk estimates to quantifying the uncertainties around these estimates. For a comprehensive assessment of lifetime risks, it is advisable to consider both the point estimate and the uncertainty interval. From this perspective, we investigate two key contributors to uncertainties: risk model parameters and baseline mortality rates. We introduce advanced methods for quantifying these uncertainties, facilitating their application in radiation-related research questions. Our results confirm the BEIR VII report finding that risk model uncertainty is a major driver of overall lifetime risk uncertainty (NRC 2006). Accounting for risk model parameter uncertainty in risk models from the Wismut cohort study as a practical example yields plausible lifetime risk uncertainty intervals that encompass the range of reported lifetime risk estimates from miners studies in the literature, as summarized in Table 4 of Kreuzer et al. (2024).
This study specifically addresses quantifying uncertainties in lifetime lung cancer risks associated with protracted occupational radon exposure. Existing tools like “CONFIDENCE,” “RadRAT,” and “LARisk,” predominantly employ risk models derived from the Atomic Bomb Survivors Life-Span Study and focus on acute radiation exposure (Preston et al. 2007). The latter two build on methods from the BEIR VII report (NRC 2006) to quantify lifetime cancer risk uncertainties. Those methods consider sampling variability in risk model parameters similar to our ANA approach, uncertainties in the risk transfer between populations, and uncertainties in the dose and dose-rate effectiveness factor (DDREF) using the delta method (Doob 1935; Hoef 2012). However, subjective variance inputs are needed in these tools for risk transfer and DDREF uncertainty assessment, which our study avoids. Results from the BEIR VII report are not comparable to the results here due to different study populations and type of exposure.
Moreover, in the context of occupational radon exposure, the present study adds to methods used in the BEIR VI report (NRC 1999), the US EPA report (Chiu and Brattin 2003), and Tomasek (2020), all of which use sampling techniques similar to the ANA approach here to address risk model parameter uncertainty. Although Tomasek (2020) does not provide direct lifetime risk uncertainty intervals, they can be derived from the results therein and align well with our findings (results not shown). Overall, this adds credibility to both approaches. However, the methodology to derive uncertainty intervals has not been extensively introduced, discussed, and compared with other approaches so far, as was done in the present study.

Sources of lifetime risk uncertainty
Following Thomas et al. (1992), three types of uncertainties arise in lifetime risk calculations: sampling uncertainty, uncertainty in choosing and deriving a suitable model structure (model uncertainty), and unspecified uncertainties like data errors and validity of assumptions, that cannot be formally specified with probability distributions. Quantifying uncertainties beyond sampling uncertainty is challenging, as acknowledged elsewhere (NRC 1999; Hoffmann et al. 2021). This analysis focuses on quantifiable sampling uncertainty. For a comprehensive overview of uncertainties, decisions, and potential errors associated with lifetime risk calculations, see NRC (1999).
The risk model defining ERR(t; Θ) is a crucial component of LEAR calculations with inherent uncertainties, arising from factors such as disease misclassification, statistical power limitations, potential confounding, and particularly from exposure assessment—a major challenge in radiation research and risk model derivation (UNSCEAR 2015). Potential measurement errors, especially from early years of uranium mining, are subject of ongoing research; e.g., exploring their effects within the Wismut cohort (Küchenhoff et al. 2018; Ellenbach and Hoffmann 2023). Recently, more focus has been given to more recent periods with low radon exposure, reflecting modern occupational scenarios with higher-quality individual exposure assessments (1960+ sub-cohort) (Laurier et al. 2020; Richardson et al. 2022; Kreuzer et al. 2023). However, assessing uncertainties in risk model derivation is beyond the scope of this study.
Previous sensitivity analyses indicate that many factors (e.g., the choice of the lifetime risk measure, latency time L or the maximum age tmax) have only limited impact on lifetime risk variability (Sommer et al. 2025). Likewise, variability in the annual radon exposure is only briefly investigated (SDC Section E) since the exposure scenario is fixed for most lifetime risk applications, especially for the important dose conversion considerations. Thus, the present study focuses on the most influential components: sampling uncertainty in risk model parameters Θ and baseline mortality rates r0(t), q0(t).

Risk model parameter uncertainty

Approximate normality assumption (ANA) approach
The ANA approach approximates the true underlying likelihood function for estimating risk model parameters using a normal distribution and follows the frequentist approach. The covariance matrix estimate Σˆ0 describes the amount of sampling uncertainty, which decreases as cohort size increases, resulting in narrower uncertainty intervals for lifetime risks. This explains why lifetime risk estimates with risk model parameters derived from the smaller, younger 1960+ sub-cohort have wider uncertainty intervals compared to those derived from the Wismut full cohort. This is particularly evident for models with BEIR VI structure (NRC 1999), with dedicated parameters for miners cohort data at higher age ranges.
The ANA method requires only parameter estimates and their covariance matrix to derive uncertainty intervals via Monte Carlo simulations. This makes the ANA approach practical and efficient, especially when complete access to cohort data is not available. Even though the covariance matrix might not always be included in publications, requesting it from the authors is often easier than requesting the entire cohort data, which typically also raises data protection concerns. While not entirely new, the idea of approximating the underlying likelihood function is rooted in earlier work (Rubin 1984) later coined to the term “Approximate Bayesian Computation (ABC)” (Beaumont et al. 2002; Sisson et al. 2020). A similar method for uncertainty quantification of lifetime thyroid cancer risk related to radiation exposure was used in Xue and Shore (2001).

Bayesian approach
In contrast to the frequentist method, Bayesian statistics incorporates prior knowledge or beliefs about model parameters through probability distributions, providing an alternative perspective on uncertainty.
Since the statistical software Epicure (Preston et al. 1993) does not support Bayesian methods, we implemented a solution in R (R Development Core Team 2024). Although R can be computationally slower than Epicure, which is specifically optimized for fitting ERR(t; Θ) risk model structures, individual R solutions allow greater control over cohort data and model fitting. We calculate the marginal posterior distribution of risk model parameters, inspired by Higueras and Howes (2018). Computing the full posterior is computationally expensive due to numerous baseline stratification parameters, but focusing on the marginal posterior simplifies this by reducing the parameter set to just the risk model parameters. Compared to Higueras and Howes (2018), this technique is here extended to handle more complex risk models and protracted low exposure scenarios. Sampling from the true marginal posterior yields more nuanced uncertainty intervals that do not rely on approximating asymptotic behavior as in the ANA approach. While point estimates of the ANA approach and the Bayesian approach with a uniform prior may differ slightly due to different internal optimization techniques and software, the overall strong agreement in uncertainty intervals confirms the validity of our Bayesian implementation.
A key strength of the Bayesian approach is the ability to integrate prior knowledge (e.g., the results from previous miners studies) through the selection of prior distributions. However, selecting these priors involves subjective judgment (Goldstein 2006). Decisions on prior distributions, in particular the degree of influence on the likelihood, have to be thoroughly considered (Greenland 2006). Non-informative priors, like uniform distributions, have minimal influence on the posterior and lead to uncertainty intervals similar to the ANA approach, especially for large cohorts, while informative priors enable a more adaptive integration of diverse cohort data. Note that labeling estimates from the Joint Czech+French cohort as "prior knowledge" is a structured test of the methodology's adaptability and the effectiveness of combining different cohort information rather than a strict integration of prior results. Users can tailor the approach to their needs, but selecting appropriate priors is a separate consideration beyond this discussion.
Although Bayesian methods offer flexibility and deeper insights, they require full access to cohort data and significant computational resources, especially for complex models with many parameters. While sampling methods like MCMC are efficient in this case due to high acceptance rates by choosing the approximate (asymptotic) normal distribution as a proposal density, the computational challenge lies in calculating the marginal posterior distribution itself.

Mortality rate uncertainty
Mortality rates also introduce notable uncertainty to lifetime risk estimates. Unlike for risk model parameter estimates, which are derived through a rigorous statistical framework (likelihood theory), the sex-averaged ICRP mixed Euro-American-Asian reference mortality rates , are presented as plain numbers sourced from a database. These mortality rates do not result from a statistical estimation process, requiring careful consideration when imposing probability distributions. We assessed this uncertainty by applying gamma distributions to baseline mortality rates r0(t) and q0(t) for all ages t incorporating observed variability in mortality rates across countries in Europe, America, and Asia from WHO data (WHO 2022), which aligns geographically with ICRP reference rates. The gamma distribution was finally chosen as it fits the observed rates for numerous age groups well (see SDC Figure S7). This allowed us to quantify how mortality rate uncertainty influences lifetime risk estimates. In the described method, each rate derived from WHO data is assigned equal weight. Consequently, observations from smaller countries are given the same weight as those from larger countries. The focus was on estimating the variability of mortality rates themselves, irrespective of the population size. Alternatively, a population-weighted approach is applied in SDC Section C.2 leading to estimates with minimal mortality rate uncertainty, largely due to the substantial influence of WHO data.
Our analysis confirmed that uncertainties in all-cause mortality rates have negligible impact (Sommer et al. 2025), while lung cancer mortality rate uncertainties resulted in uncertainty intervals similar to those from full cohort risk model parameter uncertainty.
The derivation of parameter estimates in eqn (10) and (11) from WHO data is independent from ICRP reference rates and . Using a centered distribution with the expected value (mean) set equal to the ICRP rates had low effect on lifetime risk estimates and according uncertainties (SDC Section C.1.2), so we retained the un-centered distribution for Monte Carlo simulations to avoid constraining parameter estimation.
Here, lung cancer mortality rates introduce uncertainty comparable to that of full cohort risk model parameters for the Wismut cohort study. However, including both sources of uncertainty has little effect on the uncertainty intervals compared to just considering risk model parameters. Therefore, focusing solely on risk model parameter uncertainty is sufficient. Although unintuitive, this effect can be explained by acknowledging the product structure of ERR(t; Θ) and lung cancer rates r0(t) in the LEAR calculation. The variance of products of (independent) random variables is not necessarily larger than the variances of single factors (Goodman 1960).

Interpretation of uncertainty intervals
Uncertainty intervals capture the variability in lifetime risk estimates. For risk model effects, uncertainty intervals derived with ANA methods reflect sampling variation in estimated risk model parameters and may be referred to as (approximate) Wald-type classical confidence intervals. A classical confidence interval estimates a range where a population parameter likely lies based on sample data. A 95% confidence interval implies that if a sample from the population was drawn repeatedly, approximately 95% of the resulting intervals would contain the true parameter. Bayesian credible intervals, such as Highest Posterior Density Intervals (HPDIs), represent the probability (e.g., 95%) that the true value lies within the interval, incorporating prior beliefs in risk model parameters.
While both methods provide uncertainty intervals, their interpretations differ. ANA confidence intervals are less interpretable than Bayesian credible intervals. However, this theoretical distinction is practically less relevant for purposes in radiation protection research. In particular, credible intervals with non-informative priors are often similar to confidence intervals derived with the ANA approach.
Mortality rate uncertainty intervals are more difficult to interpret due to their dependence on external data and chosen distributions. They reflect sampling variability by accounting for observed mortality rate variation and can be considered as subjective confidence intervals, similar to NRC (2006). They provide a valuable quantitative sense on mortality rate variability.
The derived uncertainty intervals for lifetime risks like the LEAR reflect the expectable range of potential values that arise from the inherent variability in each calculation component. In particular, as lifetime risks are not directly estimated from data with sampling uncertainty, the intervals should not be interpreted as classical confidence intervals.

Comparison with lifetime risk variation in the literature
The uncertainty in risk model parameters contributes significantly to the overall uncertainty in lifetime risk estimates. Our analysis using the ANA approach applied to risk models from the Wismut cohort study reveals uncertainty intervals that align well with the literature: in the context of the PUMA study (Kreuzer et al. 2024), LEAR values from various studies were recalculated and summarized. Thereby, a range for the LEAR per WLM of 2.50 × 10−4 to 9.22 × 10−4 was reported across all published risk models of uranium miners studies that include time- and age-related effect modifiers. This range translates to an equivalent range for the LEAR in % of 2.35 to 8.65. The 95% uncertainty intervals derived in this study for the parametric 1960+ sub-cohort models, particularly [3.26; 12.28] for the best-fit model (eqn 5) and [3.19; 8:22] for the simple linear model (eqn 6), correspond well to the range of point estimates reported by the PUMA study group. To convert the reported LEAR per WLM values to the total LEAR, each value was multiplied by 94 (total cumulative exposure in WLM). The derived intervals for the full Wismut cohort exhibit a weaker alignment with this range: [2.06; 4.84] for the parametric model (eqn 4) and [1.27; 4.30] for the categorical model (eqn 7). However, the 1960+ sub-cohort Wismut models are preferred to full Wismut cohort models in order to estimate lung cancer risks at low protracted exposures due to high quality exposure assessment (Kreuzer et al. 2023).
Heterogeneity in radiation risk estimates between studies may explain differences in the LEAR and can likely be attributed to diverse factors such as structural differences in cumulative exposure range, duration of employment, and methods in mortality tracking and data analysis (Kreuzer et al. 2024).
While uncertainty intervals depend on the chosen confidence level (here 95%), the close alignment of our derived intervals with literature values supports the reliability and appropriateness of our approach. This is especially remarkable as our results are solely derived from the Wismut miners cohort data as a practical example. The consistency across recognized miner studies confirms the reliability of the ANA methodology in assessing uncertainties in lifetime lung cancer risks from radon exposure. Further follow-up years for miners cohorts will refine our understanding of risk models and lifetime risks for radon-induced lung cancer and associated uncertainties.

Strengths and limitations
Our work benefits from the strengths of the German uranium miners cohort (Wismut cohort), the largest single cohort of uranium miners worldwide, representing roughly half of all miners in the pooled PUMA cohort (Rage et al. 2020). This large cohort (Kreuzer et al. 2009) allows us to achieve reliable estimates for risk model effects on lifetime risk uncertainties.
Furthermore, the present study pioneers the application of the comprehensive WHO mortality database (WHO 2022) to assess mortality rate uncertainties in radon-induced lung cancer lifetime risks. These innovative approaches, along with the implementation of advanced Bayesian techniques, expand the methodological toolbox for uncertainty assessment in this field. To facilitate the application of the Bayesian technique, we developed a R procedure for data grouping and model fitting, overcoming the limitations of existing software for this specific analysis. Transparency in reporting, with detailed descriptions of methods, statistical analysis, and results, facilitates replicability.
Both methods for assessing risk model uncertainties are reliable due to statistically grounded assumptions and show broad applicability. The ANA approach requires minimal assumptions, making it versatile. While the Bayesian approach with the analytically computed marginal posterior distribution offers stronger rigor, it has specific data needs (Poisson-distributed numbers of lung cancer deaths). However, its concept of marginal posterior distributions likely extends beyond lifetime risk assessments. These techniques, given appropriate data, can be applied in radiation epidemiology to analyze various quantitative figures derived from likelihood functions, encompassing different exposures, health outcomes, and risk models. Particularly, an extension of the presented methodologies would allow for the calculation of joint uncertainty bands for dose conversion coefficients of exposure in WLM to effective dose in mSv, offering further insights into potential variability in dose assessments.
However, limitations are acknowledged. The uncertainty intervals depend on the baseline mortality rates applied, the choice of risk model, and the data used to estimate the risk models. The Bayesian approach, while statistically rigorous and contributing to the reliability of results, was here limited in applicability to models with few parameters due to computational constraints. However, generally, better computing power allows for the analysis of more complex models.
Lifetime risk calculations inherit methodological limitations of risk model estimation, such as the incorporation of detailed smoking behavior. In this analysis, smoking is not explicitly accounted for, as the primary aim was to demonstrate the statistical methodology for deriving uncertainty intervals. However, it is important to acknowledge that smoking is the greatest risk factor for lung cancer (IARC 2004; Pesch et al. 2012), and the interaction effect between smoking and radon on lung cancer risk is not yet fully understood (UNSCEAR 2021). Smoking-adjusted (and likewise sex-specific) LEAR uncertainties could be addressed by adjusting exposure scenarios, mortality rates, and risk models and associated parameters appropriately in the presented methodologies.
Future work could explore uncertainty intervals that consider smoking- and sex-specific risk model effects with the methods introduced in this work, provided appropriate risk models. Existing established risk models of uranium miners do not reflect sex-specific risks or other individual characteristics. Data for radon effects on females are sparse. Based on findings from residential radon studies in Darby et al. (2005), we decided to assume the same ERR(t; Θ) for females and males. Particularly, lifetime risk estimates using female mortality rates should be interpreted with caution, as the models are derived from male uranium miners data. Finally, the analysis does not account for uncertainties arising from transferring risk estimates from miners cohorts to the general population (here: multiplicative risk transfer; Ulanowski et al. 2020; UNSCEAR 2021). Aligning mortality rates with the specific characteristics of the cohort data, such as using national mortality rates relevant to the cohort's origin, can help partially mitigate the risk transfer issue. The overall composite nature of a lifetime risk estimate depending on multiple independently conducted analyses limits an all-encompassing uncertainty assessment.

CONCLUSION
Uncertainty quantification is crucial for a comprehensive understanding of lifetime risk estimates. This study demonstrates that uncertainty from risk model parameter estimates explains a substantial fraction of overall lifetime risk uncertainty. Two advanced methods to derive uncertainty intervals were developed and applied. The simple ANA approach proves to be a suitable and reliable uncertainty assessment technique for most cases. The more flexible Bayesian approach offers a more nuanced view of uncertainty; however, the approach is computationally more demanding and requires full access to grouped cohort data, limiting its wider applicability. From a practical perspective, additionally accounting for uncertainties in mortality rates is less critical. The explicit choice of lifetime risk measures is negligible for uncertainty assessment. The uncertainty intervals derived in this study correspond to the range of LEAR values from different miners studies in the literature; thus, uncertainties derived by both methods are mutually confirmed. The introduced methods allow for a more complete comparison of lifetime risk estimates across uranium miners studies. These findings should be accounted for when developing radiation protection policies that are related to lifetime risks.

Acknowledgments
This manuscript was previously published on arXiv: doi:https://doi.org/10.48550/arXiv.2412.06054.

Supplementary Material