In identifying the potential risk of exposure to ionizing radiation, the excess relative risk (ERR) is a common measure to quantify the relationship between radiation exposure and the risk of cancer incidence or mortality. The crucial research on radiation-associated cancer risk assessment using the ERR have been conducted using the life span study (LSS) cohort data of Japanese atomic bomb survivors in Hiroshima and Nagasaki (Grant et al., 2017). The LSS cohort is a long-term follow-up data and has been used to investigate health-related effects, such as mortality and the incidence of radiation-associated cancers and cardiovascular diseases.

Understanding the dose–response relationship is a main goal for research on radiation-associated health effects. The dose–response function represents this relationship; thus, estimating this function has been a primary concern. In general, parsimonious models are preferred as simple parametric forms, such as linear non-threshold, linear-quadratic, and quadratic functions in radiation epidemiological studies. However, the estimates from these models are sometimes unable to reflect the uncertainty at low doses. In addition, the risk at higher doses is more influential than that at lower doses (Furukawa et al., 2016), emphasizing the need for more sophisticated, flexible, and nonparametric models to examine the low-dose effect.

One alternative for understanding the dose–response relationship is to use piecewise functions which divides the domain into several equally spaced intervals and applies a different subfunction to each interval. Furukawa et al. (2016) proposed a piecewise linear dose–response model, setting the similarity among adjacent dose intervals as the prior distribution of coefficients. Park and Kim (2020) employed the same model, but applied a Gaussian process prior to adjust rapid changes in slope at each knot. Similar Bayesian approaches have been studied by Buscot et al. (2017), Kern et al. (2005), and Holmes and Mallick (2001).

Splines defined by a sum of piecewise polynomials are more flexible methods. When using splines, the number and location of knots should be specified in advance. The important tuning parameters for splines are the number and location of knots. Kauermann and Opsomer (2011) proposed data-driven selection of the number of spline basis functions in the penalized spline regression based on a likelihood criterion. Dung and Tjahjowidodo (2017) investigated the identification of the number and location of knots in non-uniform space using B-spline functions. In addition, Dimatteo et al. (2001) used the reversible-jump Markov chain Monte Carlo to capture the number or location of knots for an exponential family. However, the radiation dose distribution is highly skewed (Grant et al., 2017) and is not a member of an exponential family. Thus the choice of equally spaced knots may be inappropriate. Therefore, we propose a knot selection method based on a Dirichlet process mixture model (DPMM), which is one of the most widely used models in nonparametric Bayesian statistics (MacEachern and Müller, 1998). McAuliffe et al. (2006) applied the DPMM to nonparametric empirical Bayes problems, and Da Silva (2007) analyzed brain magnetic resonance imaging (MRI) tissue using the DPMM. Straub et al. (2015) modeled directional data on a spherical domain utilizing the DPMM, and Ngan et al. (2015) used the DPMM for the detection of outliers.

We assume that the dose distribution follows an infinite mixture of normal distributions; thus, we can select knots as the overlapping points of two normal distributions. The DPMM is an infinite mixture model with a countably infinite number of clusters inferred from data (Teh et al., 2006; Teh, 2011). This flexibility makes the DPMM more promising than finite mixture models. After selecting the knots, the coefficients of the spline are estimated using the Bayesian inference. The Bayesian method can compute credible intervals more easily than the frequentist approach, although the computational time may be higher due to the large size of LSS data.

The remainder of this paper is organized as follows. Section 2 describes the ERR model and properties of the radiation dose distribution. In Section 3, we explains the knot selection method based on the DPMM, and then an estimation method of the spline function in a Bayesian framework. Section 4 presents an example using real data to illustrate the proposed knot selection method based on the Bayesian approach, and Section 5 provides concluding remarks.

The ERR describes the proportional risk increase above the baseline rate λ(t, d = 0, z) defined as

where RR is the relative risk, d is the exposed dose, t is the event time, and z is a vector of covariates, such as age at exposure, birth year, attained age, city etc. In particular, the ERR (2.1) is relevant for modeling the dose–response relationship with effect modification. Thus, the ERR (2.1) can be fitted to individual failure time data using partial likelihood or to person-year data using Poisson rate regression (Furukawa et al., 2016). The person-year data consists of strata of subjects on radiation exposure dose groups, sex, age at exposure etc. We count the event occurrences for exposed and non-exposed subjects, respectively and compute the sum of each person-year on each stratum (NRC, 2006). In general, ERR models (2.2) tend to be nonlinear in the parameter coefficients using a person-year data. Grant et al. (2017) stated that the ERR model for the LSS cohort data (person-year data) at the radiation effects research foundation (RERF), which has conducted health-associated studies for atomic bomb survivors in Hiroshima and Nagasaki, Japan for more than 70 years, is defined as

The ERR model consists of the dose-response function ρ(d) and the effect modification term ɛ(z) as ERR(d, z) = ρ(d)ɛ(z). For example, the effect modification ɛ(z) for the LSS data is defined in (3.7). The baseline incidence rate λ(t, d = 0, z) is an exponential function of sex, age at exposure, birth year, and other factors. RERF estimated the ERR model using maximum likelihood estimation in general estimating the baseline risk and Furukawa et al. (2016) proposed the Bayesian approach using piecewise linear dose–response model. We apply the maximum likelihood estimation for the baseline risk, find the number and locations of knots in the nonparametric Bayesian knot selection method, and then estimate the ERR in the Bayesian inference. Hence,

where Y_i is the number of event occurrences at each stratum, PY_i is the sum of person-years in the i^th row, and e^η_i is the baseline incidence rate, λ(t_i, d_i = 0, z_i). η_i for the LSS data will be described in Section 3.2. Common parametric dose–response functions for ρ(d) are as follows,

To determine a more stable dose–response relationship than the above parametric dose-response functions, we consider a spline function as smooth piecewise polynomials. This allows us to apply a polynomial function to each subinterval, and connect the pieces to construct a smooth function. The piecewise quadratic spline function is defined as,

However, it is difficult to define the correct number and location of knots, δ = {δ₁, . . . , δ_C}. The number of knots is usually determined by cross-validation. One of clear ways to determine the location of knots is to divide the domain into subintervals with an equal length. However, as depicted in Figure 1, the dose distribution is highly skewed, and has a long right tail. Thus, equally spaced subintervals may not be suitable for this ERR model.

3.1. Knot selection based on Dirichlet process mixture model

The goal of this subsection is to divide a spline domain by clustering dose observations. The DPMM and infinite mixture of Gaussian models are widely used to estimate a density function (Escobar and West, 1995; Müller et al., 1996; Rasmussen, 1999). Although it can be used for density estimation, the DPMM has also been applied to a clustering algorithm (Dahl, 2006; Yu et al., 2010; Reich and Bondell, 2011). The most significant difference between the DPMM and common clustering methods (e.g., K-means or Gaussian finite mixture model) is that the DPMM can estimate the appropriate number of clusters from data. Thus, the DPMM has an advantage because the shape of the spline curve is greatly affected by the number of knots.

Let the density function be f (x) = ∫ f (x|θ)G(θ), dθ where G=Σk=1Nπkδθk is a mixing distribution (Orbanz and Teh, 2010). This simply means that f (x) has a parameter θ, and the assigned probability for each θ_k is π_k. Then, it is written as a finite mixture model. We can extend it to an infinite mixture model with G=Σk=1∞πkδθk. Bayesian assigns the Dirichlet process (DP) as a prior for G, and f (x) is then the DPMM. The DP is a stochastic process with its marginal distribution as a Dirichlet distribution, denoted as G ~ DP(α,G₀), where α is a concentration parameter and G₀ is the base distribution. α controls π_k and θk∼iidG0. Since the DP is on the infinite dimension, instead of writing it as an explicit formula, many studies have focused on developing construction schemes. Three famous schemes are the Chinese restaurant process, the Polya urn scheme, and the stick-breaking construction. Because the computation time is shorter in the stick-breaking construction for large data, we select this scheme to generate samples from the DP. When the mixing distribution G is a normal distribution, f (x) is a DP Gaussian mixture model (DPGMM). For simplicity, f(x)=Σk=1∞πkN(μk,τk-1) for θk=(μk,τk-1). Using the stick-breaking construction, it can be summarized as follows,

By properties of the beta distribution, as α is close to zero, only the first several clusters are important. In contrast, as α increases, the assigned probabilities for each cluster become equal. We assume that the density function of the dose distribution is fd(·)=Σk=1∞πkN(μk,τk-1). After estimating (μk,τk-1) and the number of clusters, we select a knot as an overlapping point of two normal distributions that has a maximum likelihood value. We repeat it to obtain reasonable knots. This is demonstrated in Figure 2.

Let d_i be the i^th dose observation. Then, the model structure to cluster the dose observations can be written as a stick-breaking construction,

Note that the finite number N is introduced as the number of cluster components. This is because the modeling phase cannot handle infinity. Thus, by using a sufficiently large number N instead of infinity, the probabilities can be calculated. Choosing the exact value of N depends on the concentration parameter, α. Further details about this problem can be found in Ishwaran and James (2002). For this model (3.1) to be the DPGMM, f must be a normal density function with θj=(μj,τj-1). In clustering analysis, the cluster indicator variables K = (K₁, . . . , K_n) and Z = (Z₁, . . . , Z_N) indicate that K_i = j signifies ZKi=θj=(μj,τj-1). With these variables, the model (3.1) can be rewritten as

To estimate the parameters, we must obtain values from the posterior distribution of (Z,K, π|D). This can be achieved by the block Gibbs algorithm (Ishwaran and James, 2001), which sequentially updates each parameter in the following order,

In each iteration, K*={K1*,…,Km*} denotes the set of m unique values of K. The procedure to estimate the parameters is described by the following steps,

• Sample Z = (Z₁, . . . , Z_N),

  • (a) Sample Z*=(ZK1*,…,ZKm*) from the full conditional distribution f(ZKj*∣K,D),

    f(ZKj*∣K,D)∝G0(ZKj*)∏{i:Ki=Kj*}f(di∣ZKj*),         j=1,…,m

  • (b) For each k∈K-{K1*,…,Km*}, draw Z_k from the base distribution G₀.

  More specifically, we can rewrite this as an explicit formula. Let hKj*=#{i:Ki=Kj*} and G₀ be a normal-gamma distribution. The normal-gamma distribution can be used as a conjugate prior when we do not know either the mean or variance of a normal distribution (Görür and Rasmussen, 2010). There are four hyperparameters, μ₀, λ₀, α₀, and β₀, and we can write it as NG(μ₀, λ₀, α₀, β₀). Then, (a) and (b) can simply be achieved by the following procedure,

  • 1) Draw τZKj*~Gamma(α0+(hKj*)/2, β0+1/2(hKj*s+(λ0 hKj*(d¯-μ0)2)/(λ0+hKj*))) where d¯=Σ{i:Ki=Kj*}di/hKj* and s=Σ{i:Ki=Kj*}(di-d¯)2/hKj*

  • 2) Draw μZKj*~N((λ0μ0+hKj*d¯)/(λ0+hKj*), 1/{(λ0+hKj*)τZKj*})

  • 3) For each k∈K-{K1*,…,Km*},

    • Draw τ_Zk ~ Gamma(α₀, β₀)

    • Draw μ_Zk ~ N(μ₀, 1/(λ₀τ_Zk ))

• Sample K = (K₁, . . . , K_n),

  (Ki∣Z,π,α,D)∼ind∑1Nπk,i*δk(·),         i=1,…,n,

  where πk,i*∝πkf(di∣Zk). This is equivalent to sampling K_i from Multinomial(N,πi*), where πi*=(π1,i*,…,πN,i*)

• Sample π = (π₁, . . . , π_N),

  Uj~Beta (1+Mj,α+∑j+1NMk),π1=U1,         πj=Uj∏ℓ=1j-1(1-Uℓ),         for j=1,…,N,

  where M_j = #{d_i : K_i = j} and πN=1-Σj=1N-1πj. To effectively define a measure of α, U_N = 1 (Ishwaran and James, 2001). In this paper, we use a fixed α of 0.5. We iterate the above steps until the chain of (Z,K, π|D) converges. The number of clusters n is determined by choosing the mode of the posterior samples.

3.2. Bayesian inference of spline coefficients

After choosing the knots, we estimate the coefficients β of ρ(d) (2.3) using the Bayesian paradigm. We must estimate not only β, but also the coefficients of λ(t, d = 0, z) (2.2) and the effect modification term ɛ(z). For simplicity, however, we focus only on estimating β, and fix other parameter estimates from maximum likelihood estimation. We do this because our priority is to make an inference about the shape of the dose-response curves, and inferences about other parameters remain secondary. Bayesian updates the parameters via MCMC. The Metropolis algorithm is used to draw samples from the multivariate distribution of β = {β₁, . . . , β_k+c}, f (β|Y, α, γ, ϕ). The prior distribution of β is a multivariate normal distribution with mean 0 and variance 100 considering a noninformative prior. The posterior distribution of β is defined as,

where λ_i = PY_ie^η_i(1 + ρ_s(d_i)ɛ(z_i)). Here, the subscript s represents a sex-specific parameter, and ρ_s(d) is the spline defined above. For clarity, explicit formulas for e^η_i and ɛ(z_i) are defined as follows,

where a is the attained age, e is the age at exposure, yr is the calendar year, b is the birth year defined as b = (floor(yr – a) – 1915)/10, and NIC is the “not in city” indicator. We ran 30,000 MCMC and treated 5,000 samples for a burn-in period. Due to the high correlation between iterations, every 100th sample remained, and four chains were used. Thus, a total of 250 × 4 = 1,000 samples approximated the posterior distribution. In addition, we used Gelman-Rubin statistics to diagnose the convergence of the chains. As mentioned, the coefficients α = (α_1s, . . . , α₁₂) of λ₀(d) and (γ, ϕ_s) of ɛ(z) were fixed as the maximum likelihood estimates for simplicity.

4.1. Life span study

The LSS cohort pertains to atomic bomb survivors in Hiroshima and Nagasaki, Japan, and provides the risk estimates of cancers related to radiation exposure. The LSS cohort is a main source of radiation-associated risk assessment for humans. The report in Grant et al. (2017) contained information for a 52-year follow-up period, and a total of 3,079,484 person-years. The total number of subjects was 105,444, and several covariates, including city, gender, age at exposure, and attained age, were considered. The target outcome was all types of solid cancers, and more details are provided by Grant et al. (2017).

4.2. Knot selection results

Table 1 presents the estimated number of cluster components. The mode was used as the appropriate number of clusters for males and females: 7 for males and 6 for females. The estimates of the means and variances of each cluster are presented in Table 2. As illustrated in Figure 2, we selected each red point as a knot that had the maximum likelihood in the overlapping region of two normal distributions, and the results are presented in Table 3. It should be noted that the selected knots in Figure 2 appear to reflect the highly skewed dose distribution in Figure 1. Thus, as expected, the knot selection using the DPMM can be an effective method to draw data-driven knots.

4.3. Sensitivity to the choice of the concentration parameter

Since the concentration parameter α of the DP plays a vital role in determining the number of cluster components, it was important to appropriately determine this parameter. Although it can also be estimated by its posterior distribution, we fixed it to 0.5 for simplicity. However, it was necessary to verify its effect on the final clustering results described in Figure 2. We thus varied its value from 0.1 to 10. Whereas the estimated number of cluster components did not change when α had a value close to 0.5, it increased to 10 and 11 as α became 5 and 10, respectively. The variances of the added clusters, however, were much higher than those of the initially estimated clusters with α = 0.5. Therefore, the number of cluster components may not be sensitive to the choice of α.

4.4. Excess relative risk model estimates

The chains converged because the Gelman-Rubin statistics were close to 1 in Table 4. There was also the estimated β, which was represented as the posterior mean. Based on this result, we compared the ERR calculated by several models for age 70 after exposure at age 30, as illustrated in Figure 3. Overall, the estimated ERR curves for females increased more sharply than those for males. In comparing the individual curves, the spline curve (black, solid) for males was fairly similar to the linear-quadratic (green, dotted) and quadratic curves (red, dash-dotted) but different from the linear curve (blue, dashed). However, it estimated the ERR to be smaller than the other two curves. This was also a prominent feature in the curve estimation for females. In particular, the estimated ERR at doses of < 0.2 Gy is presented in Figure 4. The spline curve traced the linear-quadratic curve rather than the quadratic curve. An interesting result is that the interval estimates of the spline included zero at doses of < 0.08 Gy, indicating that there was no difference between the male exposed and unexposed group when the doses were < 0.08 Gy, whereas there was no such interval for females. The quadratic curves estimated the ERR to be very small; however, the ERR increased rapidly as the doses increased. This was due to characteristics of a quadratic function. Except for the quadratic curves the spline curves also estimated the ERR to be smaller than the other curves at lower doses.

To demonstrate the importance of the knot selection procedure, we also compared two spline curves using DPMM-based knots (black) and equally spaced knots (red) in Figures 5 and 6, respectively. In this case, the number of knots remained unchanged, and only the location of the knots was altered. Over the entire dose range, there appeared to be little difference between the two curves. At low doses, however, the red curve estimated the ERR to be smaller than the black curve. The clearest distinction was that the interval estimates of the DPMM-based spline contained zero at doses of < 0.08 Gy, whereas the interval estimates of the spline with equally spaces knots did not include zero over the entire dose range for both males and females.

In this paper, we proposed spline curve fitting of the dose–response function for the excess relative risk model commonly used in radiation epidemiology. Since the choice of the number and location of knots is crucial in splines, we proposed the Dirichlet process mixture model in selecting data-driven knots, treating knot selection as a clustering problem. When the dose distribution is highly skewed, it is particularly effective due to its flexibility in the number of clusters. The chosen knots appeared to successfully reflect the dose distribution. Then, based on these knots, we estimated the Excess Relative Risk model in the Bayesian framework using the Lifetime Span Study cohort, and compared the spline curve to other parametric dose-response functions. For both males and females, the spline curves estimated the Excess Relative Risk model to be smaller than other curves; however, the Excess Relative Risk model became similar to that of the other curves as the dose increased. The estimation of the dose-response curves was greatly affected by observations at higher doses, but the use of the spline can alleviate this problem. Although we assumed a noninformative prior for the coefficients of the spline, a Gaussian process prior with a covariance matrix depending on the knots can be used as an alternative. The limitation of this research only provided that a Bayesian knot selection is conduced separately from fitting the Bayesian excess relative model.