 Research article
 Open Access
 Published:
Additivemultiplicative hazards regression models for intervalcensored semicompeting risks data with missing intermediate events
BMC Medical Research Methodology volume 19, Article number: 49 (2019)
Abstract
Background
In clinical trials and survival analysis, participants may be excluded from the study due to withdrawal, which is often referred to as losttofollowup (LTF). It is natural to argue that a disease would be censored due to death; however, when an LTF is present it is not guaranteed that the disease has been censored. This makes it important to consider both cases; the disease is censored or not censored. We also note that the illness process can be censored by LTF. We will consider a multistate model in which LTF is not regarded as censoring but as a nonfatal event.
Methods
We propose a multistate model for analyzing semicompeting risks data with intervalcensored or missing intermediate events. More precisely, we employ the additive and multiplicative hazards model with lognormal frailty and construct the conditional likelihood to estimate the transition intensities among states in the multistate model. Marginalization of the full likelihood is accomplished using adaptive importance sampling, and the optimal solution of the regression parameters is achieved through the iterative quasiNewton algorithm.
Results
Simulation is performed to investigate the finitesample performance of the proposed estimation method in terms of the relative bias and coverage probability of the regression parameters. The proposed estimators turned out to be robust to misspecifications of the frailty distribution. PAQUID data have been analyzed and yielded somewhat prominent results.
Conclusions
We propose a multistate model for semicompeting risks data for which there exists information on fatal events, but information on nonfatal events may not be available due to lost to followup. Simulation results show that the coverage probabilities of the regression parameters are close to a nominal level of 0.95 in most cases. Regarding the analysis of real data, the risk of transition from a healthy state to dementia is higher for women; however, the risk of death after being diagnosed with dementia is higher for men.
Background
In classical timetoevent or survival analysis, subjects are under risk for one fatal event. However, subjects do not fail from just one certain type of event in some applications, but are under risk of failing from two or more mutually exclusive types of events. When an individual is under risk of failing from two different types of event, these different event types are called competing risks. One of the events censors the other, and vice versa, in these competing risks frameworks [1–3]. However, many clinical trials have revealed that a subject can experience both a nonfatal event (e.g., a disease or relapse) and a fatal event (e.g., death), where the fatal event censors the nonfatal event but not vice versa. We call these types of data semicompeting risks data [4–6].
In clinical trials, the occurrence of a nonfatal event can be detected in conjunction with possibly incessant monitoring during periodic followup. For illustration purposes of our methodologies, a dataset named PAQUID (Personnes Agées Quid) is analyzed to investigate the meaningful prognostic factors associated with dementia. These data were initially analyzed by [7] using the conventional Cox model [8, 9]. Complete descriptions of the PAQUID data can be found in Conclusions section. In this paper, we employ a semicompeting risks model where death may occur after dementia has occurred (i.e., been diagnosed), but death censors the disease. An illnessdeath model [10] is perhaps one of the most commonly and frequently used semicompeting risks models. Many studies have been conducted under semicompeting risks frameworks [4, 5, 11].
As shown in the PAQUID data, dementia can be censored informatively by death. Furthermore, an additional informative censoring process can also occur. That is, participants may be excluded from the study due to withdrawal, which is often referred to as losttofollowup (LTF). It is clear that dementia would be censored due to death; however, when an LTF is present it is not guaranteed that dementia has been censored. This makes it important to consider both cases; dementia is censored or not censored. We also note that the illness process (dementia) can be censored by LTF. This forces us to consider a multistate model in which LTF is not regarded as censoring but as a nonfatal event. Considerable studies have utilized this multistate model. For example, [12] proposed a nonparametric method to estimate the survival function associated with disease occurrences, while [6, 13] used the Cox proportional hazards model [8, 9] to estimate regression coefficients.
In the meantime, most nonfatal events are observed periodically. That is, the event time is not observed exactly but lies on an interval of the form (L,R], where L is the last time a subject visited without possessing a disease and R is the first time that the subject was diagnosed with a disease. This type of censoring is called interval censoring. We could emulate what [6] did and assume that a nonfatal event of a subject occurs uniformly on the interval (L,R]. However, using the methods proposed by [14, 15], we instead partition the interval (L,R] into a few subintervals, in which a nonfatal event can occur. Ultimately, different weights can be assigned to each subinterval. Thus, the former method corresponds to an unconditional probability approach with equal weight on all of the subintervals, whereas the latter utilizes a conditional probability approach with a specific weight, depending on the subinterval.
In our study, we use the latter method to deal with nonfatal events that are intervalcensored on an interval. In addition, we propose an additivemultiplicative model by combining the Cox proportional hazards model with the additive risk model of [16], in accordance with a multistate model. The additivemultiplicative model was initially introduced by [17] and has since been developed by a number of researchers. Scheike and Zhang [19] incorporated timevarying covariates for the additive part and timeindependent covariates for the multiplicative part. On the other hand, [18] estimated relevant parameters by considering timevarying covariates for both additive and multiplicative parts. We also consider the frailty effect as a latent variable to incorporate possible connections between events; this is done because each individual is exposed to several events, including the occurrence of illness, LTF, and death.
The rest of the paper is organized as follows. First, we explain notations and procedures for parameter estimation along with the proposed models. Second, extensive simulation studies are presented to investigate the model performances in terms of the relative bias and coverage probability of the proposed estimates. We also provide the results of real data analysis. Finally, we present a summary and concluding remarks, including some of the drawbacks of the proposed models and directions for future research.
Methods
Models
As depicted in Fig. 1, the proposed model in this study consists of five states: healthy (H), nonfatal (NF), fatal (F), losttofollowup (LTF), and unobserved nonfatal (NF(LTF)). Each state is denoted by numbers 0 through 4, respectively. A total of seven possible transitions exists in the model: 0→1, 0→2, 0→3, 1→2, 3→2, 3→4, and 4→2. However, among these transitions, both 3→4 and 4→2 (displayed as dotted lines in Fig. 1 are unobservable and should be regarded as potential transitions.
Let t be the time from study entry. Additionally, S_{t} is defined as the state that each subject can take at t≥0. Then, S_{t}∈{0,1,2,3,4}. Let \(\mathcal {A}=\{(r,s):(r,s)=(0,1), (0,2), (0,3), (1,2), (3,2), (3,4), (4,2)\}\). Also, define λ_{rs}(t) to be the transition intensity from states r to s at time t. That is,
and λ_{rs}(t)=0 for \((r,s)\notin \mathcal {A}\). As mentioned above, the data corresponding to transitions 3→4 and 4→2 are not observable, requiring the following assumptions for λ_{34}(t) and λ_{42}(t):
Assumptions (1) and (2) imply that the transition intensities of H to NF and NF to F may be the same, irrespective of the occurrence of LTF. As mentioned in Background section, given covariates z=(z_{1},z_{2},…,z_{p})^{′} and w=(w_{1},w_{2},…,w_{d})^{′}, along with frailty u, we consider additive and multiplicative models defined as
where θ_{rs}(>0) and γ_{rs}(>0) are the scale and shape parameters of a Weibull distribution, respectively, and β_{rs} and α_{rs} are vectors of the regression coefficients for the additive and multiplicative parts, respectively. Moreover, η= exp(u) is the frailty for a lognormal distribution and u is assumed to follow a normal distribution with a mean of zero and variance σ^{2}. Thus, we use a Weibull distribution as a baseline transition intensity and impose a multiplicative frailty effect on the transition intensities. Since the parameters θ_{34},θ_{42},γ_{34},γ_{42},α_{34},α_{42},β_{34}, and β_{42} related to transitions 3→4 and 4→2 should satisfy the assumptions in (1) and (2), the parameter vector estimated for the model in (3) is ζ=(θ^{∗},γ^{∗},α^{∗},β^{∗},σ^{2})^{′}, where θ^{∗}=(θ_{01},θ_{02},θ_{03},θ_{12},θ_{32}), γ^{∗}=(γ_{01},γ_{02},γ_{03},γ_{12},γ_{32}), \({\boldsymbol {\alpha }}^{*} =\left ({\boldsymbol {\alpha }}_{01}^{\prime }, {\boldsymbol {\alpha }}_{02}^{\prime }, {\boldsymbol {\alpha }}_{03}^{\prime }, {\boldsymbol {\alpha }}_{12}^{\prime }, {\boldsymbol {\alpha }}_{32}^{\prime }\right)\), and \({\boldsymbol {\beta }}^{*} = \left ({\boldsymbol {\beta }}_{01}^{\prime }, {\boldsymbol {\beta }}_{02}^{\prime }, {\boldsymbol {\beta }}_{03}^{\prime }, {\boldsymbol {\beta }}_{12}^{\prime }, {\boldsymbol {\beta }}_{32}^{\prime }\right)\). According to the model in (3), the cumulative hazard functions H_{k}(t_{1},t_{2}) (k=0,1,3,4) for leaving state k between t_{1} and t_{2} can be represented as
Based on Assumptions (1) and (2), we note that β_{34}=β_{01}, β_{42}=β_{12}, α_{34}=α_{01}, α_{42}=α_{12}, θ_{34}=θ_{01}, θ_{42}=θ_{12}, γ_{34}=γ_{01}, and γ_{42}=γ_{12} in the equations of H_{3} and H_{4}.
Parameter estimation
As shown in Fig. 1, a total of six routes can be experienced by a subject from the beginning to the end of the study. These are route 1 (0 → 0), route 2 (0 → 2), route 3 (0 → 1), route 4 (0 → 1 → 2), route 5 (0 → 3), and route 6 (0 → 3 → 2). In particular, route 5 can be classified into two paths, i.e., 0 → 3 and 0 → 3 → 4, depending on whether or not a subject experiences the unobservable NF state. Similarly, route 6 can be classified into two paths: 0 → 3 → 2 and 0 → 3 → 4 → 2. We introduce notations to define the likelihood associated with each route. Consider three random variables R, L, and T, each of which represents a time from the start of the study until the occurrence of a nonfatal event, LTF, and a fatal event, respectively. Furthermore, let \(\mathcal {H}_{0}(s)\) be the set of subjects staying in state 0 at time s. That is,
Let \(\mathcal {H}_{3,f}(s)\) be the set of subjects who have already experienced LTF at time f and remain in state 3 at time s. Then,
Now, let e_{i} be the entry time of study, a_{i} be the last time the i^{th} subject visited before a nonfatal event was observed, and b_{i} be the first time a nonfatal event is observed by the i^{th} subject for i=1,2,…,n. Consider an indicator function I_{ij}, which is 1 if subject i follows route j and zero otherwise for j=1,2,…,6. Let \(\mathcal {B}_{j}=\{i:I_{ij}=1\}\). For subject \(i\in \mathcal {B}_{1}\cup \mathcal {B}_{2}\), we have a_{i},b_{i}≥t_{i}; this is the case because a nonfatal event has not been observed before time t_{i}. For subject \(i\in \mathcal {B}_{3}\cup \mathcal {B}_{4}\), we have a_{i}<b_{i}≤t_{i}; this is the case because a nonfatal event has occurred between a_{i} and b_{i}. When subject i is a member of \(\mathcal {B}_{5}\cup \mathcal {B}_{6}\), LTF has occurred at time a_{i}, which yields a_{i}<t_{i}; however, b_{i}<t_{i} or b_{i}≥t_{i}, depending on whether an unobservable nonfatal event has occurred. Thus, t_{i} would be a censoring time for \(i\in \mathcal {B}_{1}\cup \mathcal {B}_{3}\cup \mathcal {B}_{5}\), whereas it would be a time of death for \(i\in \mathcal {B}_{2}\cup \mathcal {B}_{4}\cup \mathcal {B}_{6}\). Therefore, likelihood functions Q_{1} and Q_{2} can be constructed for routes 1 and 2, respectively. These are given as follows:
Likelihood functions can also be constructed for routes 3 and 4:
Equations (6) and (7) are derived by assuming that a nonfatal event of subject i in the set \(\mathcal {B}_{3}\cup \mathcal {B}_{4}\) can occur uniformly on the interval (a_{i},b_{i}] [6]. However, we partition the interval (a_{i},b_{i}] into several subintervals where nonfatal events could occur and assign different weights to each interval [15].

Let \(R_{i^{\prime }}\in (a_{i^{\prime }},b_{i^{\prime }}]\) be an interval for the occurrence of nonfatal events associated with subjects in routes 3 or 4. Let s_{1} be the smallest value among all \(b_{i^{\prime }}\)’s for subjects in the set \(\mathcal {B}_{3}\cup \mathcal {B}_{4}\). Let s_{2} be the smallest value among all \(b_{i^{\prime }}\)’s corresponding to subjects having \(a_{i^{\prime }}\) greater than or equal to s_{1}. This process is repeated until we have no subjects with \(a_{i^{\prime }}\) greater than or equal to s_{m} (m=1,2,…). Thus, we can have a refined set of time points
$$\begin{array}{@{}rcl@{}} 0=s_{0}< s_{1}< s_{2}<\cdots< s_{l}< s_{l+1}=\infty. \end{array} $$ 
We can define the weight \(w_{i^{\prime }m}\) at time s_{m} (m=1,2,…) for subject i^{′} in the set \(\mathcal {B}_{3}\cup \mathcal {B}_{4}:\)
$$ {\begin{aligned} w_{i^{\prime}m} = \frac{d_{i^{\prime}m}\exp\left\{H_{0}\left(e_{i^{\prime}},s_{m}{\boldsymbol{w}}_{i^{\prime}}, {\boldsymbol{z}}_{i^{\prime}}, u_{i^{\prime}}\right)\right\}\lambda_{01}(s_{m}{\boldsymbol{w}}_{i^{\prime}}, {\boldsymbol{z}}_{i^{\prime}}, u_{i^{\prime}})}{\sum_{m^{\prime}=1}^{l} d_{i^{\prime}m^{\prime}}\exp\left\{H_{0}(e_{i^{\prime}},s_{m^{\prime}}{\boldsymbol{w}}_{i^{\prime}}, {\boldsymbol{z}}_{i^{\prime}}, u_{i^{\prime}})\right\}\lambda_{01}(s_{m^{\prime}}{\boldsymbol{w}}_{i^{\prime}}, {\boldsymbol{z}}_{i^{\prime}}, u_{i^{\prime}})}, \end{aligned}} $$(8)
where \(d_{i^{\prime }m} = I(s_{m}\in (a_{i^{\prime }},b_{i^{\prime }}])\). Subsequently, likelihood functions incorporated with weight \(w_{i^{\prime }m}\) in (8) for routes 3 and 4 are given by
Finally, likelihood functions for routes 5 and 6 are given by
Therefore, based on Eqs. (4)(5) and (9)(12), the likelihood function for the parameter vector ζ is
where ϕ(·) is the probability density function of a normal distribution with a mean of zero and variance σ^{2}. In our analysis, we use the NLMIXED procedure of the SAS software to estimate ζ. For the sake of parameter estimation procedures, we define the marginal likelihood as
Then, we find the value of ζ that minimizes f(ζ)=− logm(ζ), which is referred to as \(\hat {{\boldsymbol {\zeta }}}\). Consequently, the inverse of the Hessian matrix evaluated at \(\hat {{\boldsymbol {\zeta }}}\) is defined as the estimated variancecovariance matrix of \(\hat {{\boldsymbol {\zeta }}}\). Numerical integration is required for the frailty distribution. For this purpose, we use the adaptive importance sampling [20]. Finally, we employ quasiNewton optimization, which utilizes the gradient vector and the Hessian matrix of f(ζ), to achieve the optimal solution of ζ.
Results
Simulation studies
Extensive simulation is performed to investigate the finitesample properties of the estimators proposed in Methods section. As mentioned earlier, we assume a Weibull distribution with a shape parameter of γ_{rs}=1 as the baseline transition intensity and a lognormal distribution for frailty η= exp(u), where u is generated from a normal distribution with a mean of zero and a variance of 0.01. Furthermore, we use a binary covariate for z (generated from a Bernoulli trial with a success probability of 0.5) and a continuous covariate for w (generated from a standard normal distribution). We fix the sample size n at 200 and the censoring time C at 365. A total of 500 replications is used in our simulations. The following presents the details related to the generation of random variates for the i^{th} (i=1,2,…,n) subject.

Step 0: We may allow the total number of occurrences for nonfatal events to be 24 times in a 12month period, such as 15,31,…,349,365 days. However, the actual visiting time of each subject can be different from the designated times. Hence, we add random numbers, generated from a normal distribution with a mean of zero and a variance of 9, to each designated time point. Subsequently, the actual observed time points will be defined as
$$\begin{array}{@{}rcl@{}} 0= l_{0} < l_{1i} < \cdots < l_{23,i} < l_{24}=366. \end{array} $$Let u_{01i}, u_{02i}, and u_{03i} be random numbers generated from a uniform distribution on the interval (0,1). Additionally, let R_{i}, T_{i}, and L_{i} be, respectively, the roots s of the equations:
$$\begin{array}{*{20}l} \Lambda_{01}(sz_{i}, w_{i},u_{i}) + \log(1u_{01i})&=0,\\ \Lambda_{02}(sz_{i}, w_{i},u_{i}) + \log(1u_{02i})&=0,\\ \text{and} \Lambda_{03}(sz_{i}, w_{i},u_{i}) + \log(1u_{03i})&=0, \end{array} $$where
$$\begin{array}{@{}rcl@{}} \Lambda_{0j}(sz_{i}, w_{i},u_{i}) = \eta_{i}\left[\left(\beta_{0j} z_{i}\right) s + \exp\{\alpha_{0j} w_{i}\} \theta_{0j} s^{\gamma_{0j}}\right]\ \text{for}\ j=1,2,3. \end{array} $$ 
Step 1: If C≤R_{i}∧T_{i}∧L_{i}, then the i^{th} subject is defined as being censored without experiencing a nonfatal event, i.e., \(i \in \mathcal {B}_{1}\). If T_{i}=R_{i}∧T_{i}∧L_{i}, then the i^{th} subject is defined as being dead without experiencing a nonfatal event, i.e., \(i \in \mathcal {B}_{2}\). However, if R_{i}=R_{i}∧T_{i}∧L_{i}, proceed to Step 2, and if L_{i}=R_{i}∧T_{i}∧L_{i}, proceed to Step 3.

Step 2: Let u_{12i} be a random number generated from a uniform distribution on the interval (1− exp{Λ_{12}(R_{i}z_{i},w_{i},u_{i})},1), where
$$\begin{array}{@{}rcl@{}} \Lambda_{12}(sz_{i},w_{i},u_{i})=\eta_{i}\left[(\beta_{12} z_{i}) s + \exp\{\alpha_{12}w_{i}\} \theta_{12}s^{\gamma_{12}}\right]. \end{array} $$Redefine T_{i} as the root s of the equation,
$$\begin{array}{@{}rcl@{}} \Lambda_{12}(sz_{i},w_{i},u_{i}) + \log(1u_{12i})=0. \end{array} $$If C≤T_{i}, then the i^{th} subject is defined as being censored after experiencing a nonfatal event, i.e., \(i \in \mathcal {B}_{3}\). Otherwise, the i^{th} subject is defined as being dead at time T_{i} after experiencing a nonfatal event, i.e., \(i \in \mathcal {B}_{4}\). Moreover,
 If R_{i}∈(0,l_{1i}), let a_{i}=0 and b_{i}=l_{1i}. If R_{i}∈(l_{k−1,i},l_{ki}), let a_{i}=l_{k−1,i} and b_{i}=l_{ki} for k=2,3,…,23.
 However, if R_{i}∈(l_{23,i},C), the type of path should be redefined because a nonfatal event for the subject did not occur before the time of the last observation. Thus, if C≤T_{i}, the i^{th} subject is defined as being censored without experiencing a nonfatal event, i.e., \(i \in \mathcal {B}_{1}\). Otherwise, the i^{th} subject is defined as being dead at time T_{i} without experiencing a nonfatal event, i.e., \(i \in \mathcal {B}_{2}\).

Step 3: Let u_{32i} and u_{34i} be random numbers generated from uniform distributions on the intervals (1− exp{Λ_{32}(L_{i}z_{i},w_{i},u_{i})},1) and (1− exp{Λ_{34}(L_{i}z_{i},w_{i},u_{i})},1), respectively, where
$$\begin{array}{@{}rcl@{}} \Lambda_{3j}(sz_{i},w_{i},u_{i})=\eta_{i}\left[(\beta_{3j} z_{i}) s + \exp\{\alpha_{3j}w_{i}\} \theta_{3j}s^{\gamma_{3j}}\right],\ \text{for}\ j=2,4. \end{array} $$Now redefine R_{i} and T_{i} as the roots s of the equations:
$$\begin{array}{*{20}l} \Lambda_{32}(sz_{i},w_{i},u_{i}) + \log(1u_{32i})&=0 \\ \text{and} \Lambda_{34}(sz_{i},w_{i},u_{i}) + \log(1u_{34i})&=0, \end{array} $$respectively. If C≤R_{i}∧T_{i}, the i^{th} subject is defined as being censored without experiencing a nonfatal event after LTF, i.e., \(i \in \mathcal {B}_{5}\). If T_{i}≤R_{i}, then the i^{th} subject is defined as being dead without experiencing a nonfatal event after LTF, i.e., \(i \in \mathcal {B}_{6}\). However, if R_{i}<T_{i}, move to Step 4.

Step 4: Let u_{42i} be a random number generated from a uniform distribution on the interval (1− exp{Λ_{42}(R_{i}z_{i},w_{i},u_{i})},1), where
$$\begin{array}{@{}rcl@{}} \Lambda_{42}(sz_{i},w_{i},u_{i})=\eta_{i}\left[(\beta_{42} z_{i}) s + \exp\{\alpha_{42}w_{i}\} \theta_{42}s^{\gamma_{42}}\right]. \end{array} $$Redefine T_{i} as the root s of the equation,
$$\begin{array}{@{}rcl@{}} \Lambda_{42}(sz_{i},w_{i},u_{i}) + \log(1u_{42i})=0. \end{array} $$If C≤T_{i}, then the i^{th} subject is defined as being censored at time C after experiencing LTF and a nonfatal event, i.e., \(i \in \mathcal {B}_{5}\). Otherwise, the i^{th} subject is defined as being dead at time T_{i} after experiencing LTF and a nonfatal event, i.e., \(i \in \mathcal {B}_{6}\).
In the simulation settings, we consider three types of regression coefficients (i.e., ‘even’, accelerated (‘acc’), and decelerated (‘dec’)) as well as three types of LTF proportions (i.e., ‘low’, ‘moderate’, and ‘high’). For the ‘even’ type, there are no differences in the effects of the covariates on the hazard rate of death before and after experiencing a nonfatal event. That is, α_{02}=α_{12}=0.01 and β_{02}=β_{12}=0.004. Meanwhile, for ‘acc’ and ‘dec’, increasing and decreasing effects are noted on the hazard rates of death, respectively. That is, α_{02}=0.01, α_{12}=0.0125, β_{02}=0.004, and β_{12}=0.005 for ‘acc’, whereas α_{02}=0.02, α_{12}=0.01, β_{02}=0.008, and β_{12}=0.004 for ‘dec’. For the rest of the regression coefficients, we set β_{01}=β_{03}=β_{32}=0.004 and α_{01}=α_{03}=α_{32}=0.01. Moreover, we set θ_{03}=0.00075, θ_{03}=0.002, and θ_{03}=0.004 for the ‘low’, ‘moderate’, and ‘high’ types, respectively. The remaining shape parameters of the baseline transition intensities were set as θ_{01}=θ_{32}=0.002 and θ_{02}=0.001. Tables 1, 2, and 3 provide the relative bias (‘r.Bias’), standard deviation (‘SD’), average of the standard errors (‘SEM’), and coverage probability (‘CP’) of 95% confidence intervals for the regression parameters and the variance estimate of the frailty distribution, respectively, according to the three LTF proportions. For comparison purposes with the proposed approach (‘proposed’), each table also displays the results obtained by simply assuming that a nonfatal event occurred at the end of the right endpoint of the interval (‘imputedbytherightendpoint’). When the type of the regression coefficients is fixed at ‘even’, ‘acc’, or ‘dec’, the CPs of the regression parameters corresponding to the ‘proposed’ case are close to a nominal level of 0.95 irrespective of the LTF proportions, whereas those of the regression parameters such as b_{01} and b_{03}, are much smaller than 0.95 for the ‘imputedbytherightendpoint’ case. For the results based on the ‘proposed’ method, as the proportion of LTF increases, the mean squared error (MSE) for estimates of some regression parameters (e.g., a_{03},a_{32}, and b_{32}) decreases, while the MSE of other regression parameters (e.g., a_{02},b_{02}, and b_{03}) increases, regardless of the type of regression coefficients.
Sensitivity analysis is also conducted to investigate how the estimator of the parameter behaves with different frailty distributions. For simplicity of computation, we consider only the ‘even’ case for the regression parameters and the ‘moderate’ LTF proportion. Three different frailty distributions are used, along with a normal distribution with a mean of zero and a variance of 0.01. These are uniform, double exponential, and gamma distributions with specific parameter value(s), for which the mean and variance of each distribution are the same as those of the normal distribution. Simulation results are provided in Table 4. We compare the results of the three distributions with those of the normal distribution. The uniform and double exponential distributions are symmetric, like the normal distribution. However, the uniform distribution has thinner tails than the normal distribution, while the double exponential distribution has heavier tails than the normal distribution. Alternatively, unlike the normal distribution, the gamma distribution is an asymmetric distribution. Overall, there are no differences in the values of r.Bias and CP between the three distributions and the normal distribution. This implies that the proposed estimators are robust to misspecifications of the frailty distribution.
Illustrative data analysis
PAQUID data were collected to investigate the effects of dementia on mortality. Samples were taken from community residents of two southwestern regions (Gironde and Dordogne) of France [7]. The population consists of elderly people of ages 65 or above, between 1988 and 1990, whose sociodemographic characteristics and mental health status were recorded every two to three years. A total of 3675 persons was selected to participate in the study; among these individuals, 832 (22.6%) were diagnosed with dementia, 639 of whom died. The remaining 2843 participants (77.4%) did not experience dementia but 2298 of them died.
In this article, we performed an analysis based on ‘paq1000’ data, which included 1000 randomly selected observations from the PAQUID data [21]. The paq1000 data consist of several pieces of information, such as the mental health status (diagnosed with dementia or dementiafree), dead or alive status, age (including a n’s age at the start of study, their age at the last dementiafree visit, their age when they were diagnosed with dementia, their age at their time of death, and their age at censoring), gender, and educational background (educated or noneducated in terms of graduation of elementary school, say certificate). When a person who was not diagnosed with dementia at their last visit has not been traced for more than four years, this person is assigned to the LTF category. Table 5 shows summary statistics to briefly grasp the subjects’ characteristics including ages at entry, at dementia diagnosis, at death after dementia, at death without dementia, and at death after LTF. A total of 231 persons was categorized as LTF; among these, 159 (68.8%) died. Moreover, 127 (68.3%) persons out of the 186 who experienced dementia died, and 438 (75.1%) persons out of the 583 dementiafree persons died. Moreover, age at dementia diagnosis is higher for women than for men regardless of the education group (with certificate or without certificate). The same trend is observed both at age after dementia and at age without dementia. Meanwhile, age at death after LTF is a bit higher for women than for men.
First, we check to see whether each covariate satisfies the proportional hazards assumption by using the test procedure of [22] and the Schoenfeld residual [23]. Figure 2 shows diagnostic (scattered) plots of the scaled Schoenfeld residual versus age. In each plot, we mark a spline smoother (solid line) as well as two standard deviation bands (dashed lines). In the curve showing the effect of gender, there is a decreasing trend on transitions 0→1 and 0→2, an increasing trend on 1→2 and 3→2, and a quite steady pattern for 0→3. In the curve showing the effect of certificate, there is a prominent decreasing trend for transition 0→2, while the other transitions show steady patterns. Table 6 provides the pvalues for the test results [11]. For the gender effect, only the pvalue for transition 0→3 is greater than 0.1, which seems to violate the proportional hazards assumption. Alternatively, the pvalues for all transitions for the certificate effect, with the exception of 0→2, are greater than 0.1. This seems to satisfy the proportional hazards assumption. Thus, we put ‘gender’ into the additive side and ‘certificate’ into the multiplicative side for future analyses:
for \((r,s)\in \mathcal {A}.\)
Table 7 shows a summary of the estimation procedures. We provide the estimates of the regression coefficients along with their standard errors and pvalues. For the transition from H to LTF (0→3), the intensities for women are larger than for men, with a value of 0.00849 (849 out of 100000 persons). However, this turned out to be insignificant with a pvalue of 0.439. For the intensity of the 3→2 transition, a reversed outcome was obtained, with a value of 0.00619 for men with a very significant pvalue of less than 0.001. For the 0→1 transition, women showed a larger intensity than men with a value of 0.0156, yielding a nonsignificant result with a pvalue of 0.245. For the 1→2 transition, the intensity for men is similar to that for women, with a value of 0.0101 and a pvalue of 0.961. Finally, for the 0→2 transition, the intensity for men is larger than for women with a value of 0.0295 along with a significant pvalue of 0.038. Meanwhile, all transition intensities of the noneducated group are similar to those of the educated group. Finally, the estimate for the variance σ^{2} on the common frailty is 0.999 with a highly significant pvalue less than 0.001, showing nonhomogeneity between clusters classified by the age at entry. Figure 3 shows five transition intensities over age by gender and certificate and estimated normal frailties of each cluster. As presented in Fig. 3, the transition intensities of 0→1 and 0→3 are higher for women than for men over age regardless of the education group, while these trends are reversed for the 0→2 and 3→2 transitions. For the transition 1→2, no difference is observed between women and men, but there is a significantly monotone increasing trend over age. These are consistent with the results in Table 7.
Conclusions
We considered a multistate model for semicompeting risks data, for which there exists information on fatal events, but information on nonfatal events may not be available due to lost to followup. More precisely, we proposed an additive and multiplicative random effect model by combining the additive risk model [16] with the proportional hazards model [8, 9] in order to derive the conditional likelihood function. An adjusted importance sampling method was used to compute the marginal likelihood function, where the MLEs for the regression coefficients were obtained by an iterative quasiNewton algorithm. The proposed model was illustrated using PAQUID data and yielded several promising results. The risk of transition from a healthy state to dementia is higher for women. However, the risk of death is higher for men regardless whether a subject is diagnosed with dementia or not. Meanwhile, the risk of transition from a healthy state to dementia is higher for the educated group. The risk of death after being diagnosed with dementia is higher for the educated group; however, a reversed result is observed for nondiagnosed subjects. Furthermore, we conducted simulations with finitesample sizes to investigate the efficiency of the proposed estimators. In particular, we considered nine combinations of three different types of regression coefficients and three different types of LTF proportions. In general, the coverage probabilities of the regression parameters are close to a nominal level of 0.95 in most cases. Moreover, according to a referee’s suggestion, we investigated influence of parameter estimation when LTF is omitted (not censoring) in our simulation studies. Finally, we performed simulations with the same realizations generated from each configuration included in Tables 1, 2, and 3. Based on the results not reported here, the CP of parameter β_{01} is extremely lower than a nominal level of 0.95 for all configurations. In addition, when the type of regression coefficients is ‘dec’, the CP of parameter β_{02} is much smaller than 0.95 regardless of types of LTF proportions. Moreover, compared to each table, namely, Tables 1, 2, and 3, the relative bias of β_{01} increases around ten times for all configurations. This is the reason omitting LTF results in route changes of subjects included in routes 5 and 6 at the data generation stage, i.e., from route 5 to 1 (when the fatal event is censored) or from route 6 to 2 (when the fatal event is observed) according to the fatal status of subjects.
The proposed model has some drawbacks. At the initial stage of this research, as developed in [8, 9, 16], we intended to consider a semiparametric model incorporating five transition intensity models. However, it was quite difficult to handle nonparametric estimation procedures for the baseline transition intensity associated with the nuisance parameter because the total number of parameters is proportional to the number of subjects. Rather, we assumed a Weibull distribution for the baseline transition intensity; further research should be carried out to avoid the use of this specific distribution. To circumvent arbitrariness, it is necessary to calculate the NelsonAalen estimators for the cumulative baseline transition intensity [8, 9, 16] and extend this method to semicompeting risks models based on the profile likelihood function. Another plausible remedy would be to apply a spline smoothing method on the baseline transition intensity proposed by [24, 25]. In the additive and multiplicative hazards model with frailty, one could use a semiparametric Bayesian approach by assuming a prior distribution on the lognormal frailty. Subsequently, conventional Markov Chain Monte Carlo computation would be proceeded for the full conditional distribution on the frailty [26, 27].
Abbreviations
 CP:

Coverage probability
 F:

Fatal
 H:

Healthy
 LTF:

Losttofollowup
 MSE:

Mean squared error
 NF:

Nonfatal
 PAQUID:

Personnes Agées Quid
 r.Bias:

Relative bias
 SD:

Standard deviation
 SEM:

Average of the standard errors
References
Putter H, Fiocco M, Geskus RB. Tutorial in biostatistics: competing risks and multistate models. Stat Med. 2007; 26:2389–430.
Lau B, Cole SR, Gange SJ. Competing risk regression models for epidemiologic data. Am J Epidemiol. 2009; 170:244–56.
Berry SD, Ngo L, Samelson EJ, Kiel DP. Competing risk of death: an important consideration in studies of older adults. J Am Geriatr Soc. 2010; 58:783–7.
Fine JP, Jiang H, Chappell R. On semicompeting risks data. Biometrika. 2001; 88:907–19.
Xu J, Kalbfleisch JD, Tai B. Statistical analysis of illnessdeath processes and semicompeting risks data. Biometrics. 2010; 66:716–25.
Barrett JK, Siannis F, Farewell VT. A semicompeting risks model for data with intervalcensoring and informative observation: An application to the MRC cognitive function and aging study. Stat Med. 2011; 30:1–10.
Helmer C, Joly P, Letenneur L, Commenges D, Dartigues JF. Mortality with dementia: results from a French prospective communitybased cohort. Am J Epidemiol. 2001; 154:642–8.
Cox DR. Models and lifetables regression. J R Stat Soc Series B. 1972; 34:187–220.
Cox DR. 1975. Biometrika; 62:269–76.
Andersen PK, Borgan Ø, Gill RD, Keiding N. Statistical Models Based on Counting Processes. New York: Springer; 1993.
Day R, Bryant J, Lefkopoulou M. Adaptation of bivariate frailty models for prediction, with application to biological markers as prognostic indicators. Biometrika. 1997; 84:45–56.
Frydman H, Szarek M. Nonparametric estimation in a Markov illnessdeath process from interval censored observations with missing intermediate transition status. Biometrics. 2009; 65:143–51.
Siannis F, Farewell VT, Head J. A multistate model for joint modeling of terminal and nonterminal events with application to Whitehall II. Stat Med. 2007; 26:426–42.
Lindsey JC, Ryan LM. Tutorial in biostatistics: Methods for intervalcensored data. Stat Med. 1998; 17:219–38.
Collett D. Modeling Survival Data in Medical Research, 3rd ed. London: Chapman and Hall/CRC; 2015.
Lin DY, Ying Z. Semiparametric analysis of the additive risk model. Biometrika. 1994; 81:61–71.
Andersen PK, Værth M. Simple parametric and nonparametric models for excess and relative mortality. Biometrics. 1989; 45:523–35.
Scheike TH, Zhang M. An additivemultiplicative CoxAalen regression model. Scand J Stat. 2002; 29:75–88.
Martinussen T, Schieke TH. A flexible additive multiplicative hazard model. Biometrics. 2002; 89:283–98.
Pinheiro JC, Bates DM. Approximations to the loglikelihood function in the nonlinear mixedeffects model. J Comput Graph Statist. 1995; 4:12–35.
Touraine C, Joly P, Gerds TA. SmoothHazards: Fitting illnessdeath model for intervalcensored data. 2014. https://CRAN.Rproject.org/package=SmoothHazard.
Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika. 1994; 81:515–26.
Schoenfeld D. Partial residuals for the proportional hazards regression model. Biometrika. 1982; 69:239–41.
Joly P, Commenges D, Letenneur L. A penalized likelihood approach for arbitrarily censored and truncated data: application to agespecific incidence of dementia. Biometrics. 1998; 54:185–94.
Leffondré K, Touraine C, Helmer C, Joly P. Intervalcensored timetoevent and competing risk with death: is the illnessdeath model more accurate than the Cox model?. Int J Epidemiol. 2013; 42:1177–86.
Duchateau L, Janssen P. The Frailty Model. New York: Springer; 2008.
Ibrahim JG, Chen MH, Sinha D. Bayesian Survival Analysis. New York: Springer; 2001.
Acknowledgements
The authors would like to thank two reviewers for their constructive comments and suggestions that greatly contributed to improving the earlier version of the paper. They would also like to thank the Editors for their generous comments and support during the review process.
Funding
This work was supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (NRF2017R1D1A1B03028535). The funding body had a role in the design of the study, analysis, and interpretation of ‘paq1000’ dataset, and in writing the manuscript as well.
Availability of data and materials
The datasets generated and/or analysed during the current study are available in the R ‘SmoothHazard’ package repository, https://CRAN.Rproject.org/package=SmoothHazard.
Author information
Authors and Affiliations
Contributions
JHK developed the estimation procedure for the proposed model and designed the simulation settings. JYK performed computational works related to simulations and real data analysis. SWK organized and wrote this manuscript and discussed the theoretical methodologies with JHK. All authors have read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Kim, J., Kim, J. & Kim, S.W. Additivemultiplicative hazards regression models for intervalcensored semicompeting risks data with missing intermediate events. BMC Med Res Methodol 19, 49 (2019). https://doi.org/10.1186/s128740190678z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s128740190678z
Keywords
 Additive and multiplicative hazards model
 Interval censoring
 lognormal frailty
 Missing intermediate event
 Multistate model
 Semicompeting risks data