A Bayesian Nonlinear Mixed-effects Disease Progression Model

A nonlinear mixed-effects approach is developed for disease progression models that incorporate variation in age in a Bayesian framework. We further generalize the probability model for sensitivity to depend on age at diagnosis, time spent in the preclinical state and sojourn time. The developed models are then applied to the Johns Hopkins Lung Project data and the Health Insurance Plan for Greater New York data using Bayesian Markov chain Monte Carlo and are compared with the estimation method that does not consider random-effects from age. Using the developed models, we obtain not only age-specific individual-level distributions, but also population-level distributions of sensitivity, sojourn time and transition probability. Journal of Biometrics & Biostatistics J o u rn al of Bio metrics & Bistatis t i c s


Introduction
The disease progression model introduced by Zelen and Feinleib [1] assumes that disease progresses through three states: where So corresponds to the disease-free state, Sp the preclinical disease state when an asymptomatic individual unknowingly has the disease that a screening exam can detect, and Sc the clinical state when the disease manifests itself in clinical symptoms. In addition, if one enters the preclinical state (Sp) at time t_1 and the clinical state (Sc) at time t 2 , the time difference (t 2 -t 1 ) is called sojourn time. The length of the time (t 2 -t) is called lead time if one is offered a screening exam at time t within the preclinical state (t 1 ≤ t ≤ t 2 ) and the disease is diagnosed.
Screening aims to detect disease in the preclinical state before symptoms appear, which may greatly increase the chances for effective treatment. The disease progression model has been used to analyze the screening data with three key components [2,3]: sensitivity of the screening modality, the sojourn time distribution, and the transition probability density from the disease-free to the preclinical state. These three parameters are building blocks for screening modeling, because all other parameters of interest can be expressed as a function of these key parameters. In particular, the sensitivity of the screening modality is critical for evaluating the predictive performance of a screening program. Note that the sensitivity is defined as the conditional probability that a screening test is positive when one is in the preclinical state (Sp).
Screening sensitivity may depend on a variety of factors, such as position, location and size of the tumor, experience of the radiologist, age, etc. Wu et al. [2] modeled the sensitivity as a function of age with an age-dependent transition probability density and then applied the model to the Health Insurance Plan of Greater New Yorker (HIP), a breast cancer screening study, and later to the Mayo Lung Project data, a lung cancer screening study [3]. However, sensitivity of mammogram increases as a woman's age increases, while screening sensitivity of chest X-ray for lung cancer does not depend on age. In later developments, sensitivity was modeled as a function of the time spent in the preclinical state along with age at diagnosis [4], but sensitivity was influenced by the proportion of time in the preclinical state to the sojourn time more than by age for breast cancer screening. For this reason, Kim and Wu [5] recently modeled sensitivity as a function of the sojourn time and time spent in the preclinical state and then applied the model to the Johns Hopkins Lung Project data. Nevertheless, it still remains unclear which of these models describes sensitivity more appropriately.
Besides, the progress of disease can vary by age. In particular, it is well known that human cancer incidence depends on age, and the risk of being diagnosed with cancer increases with age [6]. Therefore, it is desirable to estimate the results of disease (e.g., cancer) screening for overall summaries as well as for age-specific summaries. To our knowledge, disease progression models have not considered variation in age.
The main objectives of this study were to resolve the two aforementioned issues: (i) finding a proper sensitivity model and (ii) estimating the disease progression models by considering variation in age. To this end, we generalize sensitivity as a function of age at diagnosis, sojourn time, and time spent in the preclinical state and a nonlinear mixed-effects model is developed for disease progression models from a Bayesian framework. We then applied our models to the Johns Hopkins Lung Project (JHLP) data and the Health Insurance Plan for Greater New York (HIP) data. All simulations were run by using the statistical software R (R Development Core Team), and the algorithms described in this study can be obtained upon request from the authors.
The remainder of the paper is organized as follows. In Section 2, we introduce a disease progression model and our generalized sensitivity model. A Bayesian nonlinear mixed-effects model is developed based on a trinomial distribution in Section 3. In Section 4, the developed models are applied to JHLP and HIP data. Concluding remarks are found in Section 5.

( )
where κ and are positive scale and location parameters in the loglogistic family. For more detailed descriptions of the log-logistic distribution, refer to Kim and Wu [5].
The disease progression model with the proposed sensitivity model has the parameter θ = (α,β,γ,,σ 2 ,κ,ρ). For the cohort aged t i,0 at study entry, the conditional likelihood function is where N i is the number of screenings. As a result, the overall likelihood function for all the study groups is the product of the age-specific contributions across all age groups The conventional parameter estimation ignores the age-specific effect since its likelihood function is based on Equation (3). One approach to incorporating the age-specific effect into a model is to use a mixed-effects model. In the next section, a nonlinear mixed-effects model is developed with the age-specific effect as a random-effect.

A Bayesian Nonlinear Mixed-Effects Model
The three-stage hierarchical nonlinear mixed-effects model is developed for disease progression models from a Bayesian framework. The first-stage model has the form of ( ) ( ) ij ij i ij ij ij i i s , r |~Tri n , D , I | ,1 i K;1 j N , where D i,j is the probability of an individual correctly diagnosed at the jth scheduled exam given at the ith age group; I i,j is the probability of an interval case in (t i,j-1 ,t i,j ) at the ith age group; n ij is the total number of the individuals examined at t i,j-1 at the ith age group; S ij : the number of cases diagnosed at t i,j-1 at the ith group; r ij : the number of interval cases within the interval (t i,j-1 ,t i,j ) at the ith age group; and θ i is a p-dimensional age-specific parameter vector, where, in this project, p=7 and Note that Tri stands for a trinomial distribution. The second-stage is structured based on a multivariate normal distribution (MVN) and is given by where θ is the population-average parameter vector, and ∑ p is the between-age covariance matrix. The third stage of the model describes the prior distributions for given c, C p , , and R p , where W is a Wishart distribution. The posterior distribution for (θ 1 ,… θ k , θ, ∑ p ) is proportional to The previous published results in Kim represent N i ordered screening exams starting at age, t i,0 where D i is the fixed follow-up time after the last examination. The jth screening interval is (t i,j-1 ,t i,j ), j = 1,2,…N i at the ith age group. The following notation is used: the jth annual screening exam occurs at age t i,j-1 = t i,0 + j-1, for j = 1,2,…N i , and t i,-1 = 0; n ij is the total number of individuals examined at t i,j-1 ;s ij is the number of cases diagnosed and confirmed at t i,j-1 ; and r ij is the number of interval cases within the interval (t i,j-1 ,t i,j ).
For each age group, we model sensitivity to vary with three factors, which are screening age t, sojourn time T, and time spent in the preclinical state s, by Where t is the average age at entry in the entire study group and s ∈ [0,T] is the time spent in S p . This is an extension of the sensitivity of Kim and Wu [5] where the sensitivity depends on the sojourn time and the time spent in the preclinical state. Note that sojourn time T is a random variable in this model. Here, in general, the parameters α and γ are responsible for the maximum value and for the rate of the sensitivity, respectively, while the parameter β explains how the behavior of the sensitivity changes with age. Namely, the maximum sensitivity increases as the parameter α increases. When s/Tis close to zero, sensitivity increases rapidly if γ < 1, while sensitivity increases gradually if > 1. Sensitivity is an increasing function of age when the parameter β is positive.
Let D ij be the probability of an individual correctly diagnosed at the jth scheduled exam, given at t i,j-1 and started the screening exam at age t i,0 (i.e. the ith age group), and I ij the probability of an interval case in (t i,j-1 ,t i,j ). These two probabilities, for j = 1,2,…N i , are: is the probability density function for transition from S o to Sp at age X and is modeled as a sub-density of a log-normal distribution, where 20% was selected based on the previous analysis on Lung cancer screening [3]; q(t) is the probability density function (pdf) Because of the nonlinear functions D ij and I ij , the conditional distribution of θ 1 is non-standard and known up to a normalizing constant, where L i (θ 1 ) denote the conditional log-likelihood for the ith age of Equation (2).

Applications
The developed methods are applied to both the Johns Hopkins Lung Project (JHLP) data [5,7] and the Health Insurance Plan for Greater New York (HIP) data [8]. Since no analytical formulas were available for the likelihood function, we used the Markov chain Monte Carlo (MCMC) approach to estimate the posterior distribution for each parameter. In details, let the full conditional posterior distribution of a vector of parameters, θ, be [ | ] θ ⋅ , that is, the posterior distribution of θ given all other quantities in the model. Since our model is nonlinear, it leads to a form for [ | ] θ ⋅ that is non-standard and is known only up to a normalizing constant. For this reason, we updated the parameters via Metropolis-Hastings-within-Gibbs steps and chose random-walk Metropolis for updating the θ 1 parameters, which is a natural choice of Metropolis-Hastings [9]. For setting of the prior distribution (hyperparameters), we employed the previous results in Kim and Wu [5] and Wu et al. [2] as stated before.
For comparison, we analyzed the JHLP and HIP data using both the conventional and developed approaches. Here the conventional approach means that all the parameters were estimated based on the overall likelihood, which is Equation (3), without considering the variation in different age groups. For convenience, we call the developed algorithm the mixed-effects disease model (ME-DM) and the conventional algorithm the fixed-effects disease model (FE-DM).
For both ME-DM and FE-DM, we generated one MCMC chain for each of JHLP and HIP data until the MCMC chain reached at least 20,000 iterations. We then obtained the MCMC chain with the size of 1,000 using the last 10,000 iterations for each chain with burnin of at least 10,000 and a thinning of every 10 steps. Supplementary  Information (Figures S1 and S2 in Appendix) are the trace plots of the simulated Markov chains only for population-level parameters. The solid and dotted horizontal lines represent the mean and the median of each parameter, respectively.

The jhlp and hip data
In the Baltimore metropolitan area from 1973 to 1978, the Johns Hopkins Lung Project (JHLP) trials enrolled 10,386 men aged 45 years and older who smoked at least one pack of cigarettes per day (or who had smoked this much within one year of enrollment) and who had no prior history of respiratory cancer. All participants were then randomized to either chest X-ray only or a dual screen (chest X-ray and sputum cytology) groups, resulting 5,160 men in the chest X-ray only arm and 5,226 in the dual-screen arm. Participants in the chest X-ray arm received chest X-ray screening test annually, for 8 consecutive years. If any of the tests were positive, the screen was considered positive and a definitive work-up exam, such as biopsy, was done. In this study, we used the chest X-ray arm, including the total number of participants in each screening exam, the number of detected and confirmed cancer cases in each screening exam, and the number of interval cases. These data were stratified by age at study entry from 45 to 88 years old.
The Health Insurance Plan of Greater New York (HIP) study began at the end of 1963 and was the first randomized clinical trial for regular screening exams that include mammography as a screening test for breast cancer [8]. The study enrolled asymptomatic women aged 40 to 64 years who had no history of breast cancer. The participants were randomized into study and control arms, with about 31,000 women in each arm. The screening program for the study arm specified up to four annual breast exams with both a mammogram and a clinical breast exam, while the control arm received usual care. Data from the study arm was used for this study, where the data were stratified by age at study entry from 40 to 64 years old. Table 1 displays the empirical means, standard deviations, and 95% credible intervals (CIs) of the posterior distributions of parameters. As for ME-DM, the fixed-effect estimates (i.e., population-level estimates) are considered to compare with the estimates of FE-DM. We can see the large absolute difference in log (α) and log (γ) between two estimation methods. However, all 95% CIs of estimates of FE-DM overlap with these of population-level estimates of ME-DM except for the estimate of JHLP and HIP and the estimate ñ of HIP (Table 1).

Results
Estimates of the variance-covariance matrix ∑ of ME-DM are shown in Table 2. Since only log (α), β, log (γ), and m are considered as random-effects, the size of ∑ is four by four. For both JHLP and HIP data, there is greater variation in the parameters log (α) and log (γ) than these in other parameters, indicating that sensitivity is influenced by age at diagnosis. Forest plots of each individual-level estimate of ME-DM are plotted in Figure 1. In case of the parameters β and µ, the empirical means of the individual-level estimates are very close to that of the population-level estimate for both JHLP and HIP data. On the other hand, we can see a larger variation of the individual-level estimates of log (α) and log (γ). These imply that the parameters α and β have a large influence on age, so does sensitivity. The individuallevel estimates of each age at diagnosis can be found in Supplementary Information (Tables S1 and S2 In appendix).
The developed sensitivity models depend on age at diagnosis, the time spent in the preclinical state and the sojourn time, resulting in a function of age and the proportion of time spent in the preclinical state to the sojourn time. Note that the average age in Equation (1) is globally set to 55 years for all age groups in both JHLP and HIP data. ( ) logˆ0 γ < ), while these of ME-DM are greater than one (i.e., ( ) logˆ0 γ > ), for both JHLP and HIP data ( Table  1). As a result, the population-level posterior sensitivities estimated by FE-DM are larger than these by ME-DM when s/T is near zero, based on the subplots between η and s/T in Figure 2. Interestingly, many individual-level posterior estimates of γ are less than one (i.e., ( ) logˆ0 γ < ), although the population-level estimates are larger than one (i.e., ( ) logˆ0 γ > ), as can be seen in Figure 1 and Supplementary  Information Tables 1 and 2. On the other hand, the population-level sensitivities of JHLP are an increasing function in age regardless of estimation methods, while these of HIP are a decreasing function in age in ME-DM , based on the subplots between age and s/T in Figure 2. This is because the posterior estimates of JHLP and HIP are positive and negative, respectively. In general, both JHLP and HIP data show large differences in sensitivity between FE-DM and ME-DM.
The individual-level posterior sensitivities are shown in Supplementary Information (Figures S3 and S4 in appendix). In particular, these predicted sensitivities show significant variations among age groups, which might be resulted from the large variations in parameters log (α) and log (γ) in Table 2.

Concluding Remarks
We propose a generalized sensitivity model which depends on age at diagnosis, time spent in the preclinical state and sojourn time, and the developed sensitivity model is applied to a novel nonlinear mixedeffects model for a disease progression model in a Bayesian framework.
As for JHLP data, FE-DM along with the developed sensitivity model estimates the parameters using the same data as used in Kim and Wu [5]. The main difference is the sensitivity model, and their sensitivity model is a special case of our model. That is, when the parameter b is equal to zero in Equation (1), our model becomes same as their model. Generally, our results are consistent with these of Kim and Wu [5]. However, our estimate of γ is smaller than their estimate, 0.01 vs. 0.13, respectively (Table 3 in Appendix), resulting that our sensitivity increases much faster than theirs. This might be due to ageeffect in the sensitivity.
When the population-level sensitivity of JHLP data is compared with that of HIP data, we can see a different trend from each other in the sense that the sensitivity increases as men get older in lung cancer, while the sensitivity decreases as women get older. In other words, the probability to detect lung cancer is higher in old than in young, and     the breast cancer might be detected more in young than in old by a screening exam. This can be explained by the fact that lung cancer does not have a reliable early detection test compared with other cancers [10].
The main advantage of ME-DM over FE-DM is to incorporate variation in age into the disease progression model, resulting in estimates with better precision. By doing so, we can obtain not only population-level estimates but also individual-level estimates. Accurate estimates are critical for policy makers to predict the performance of a screening exam. In this regard, the proposed ME-DM along with a generalized sensitivity model might provide a more accurate assessment of screening for policy makers [11].