Additive-multiplicative hazards regression models for interval-censored semi-competing risks data with missing intermediate events

Background In clinical trials and survival analysis, participants may be excluded from the study due to withdrawal, which is often referred to as lost-to-follow-up (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 multi-state model in which LTF is not regarded as censoring but as a non-fatal event. Methods We propose a multi-state model for analyzing semi-competing risks data with interval-censored or missing intermediate events. More precisely, we employ the additive and multiplicative hazards model with log-normal frailty and construct the conditional likelihood to estimate the transition intensities among states in the multi-state 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 quasi-Newton algorithm. Results Simulation is performed to investigate the finite-sample 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 multi-state model for semi-competing risks data for which there exists information on fatal events, but information on non-fatal events may not be available due to lost to follow-up. 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 time-to-event 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 In clinical trials, the occurrence of a non-fatal event can be detected in conjunction with possibly incessant monitoring during periodic follow-up. 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 semi-competing 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 semi-competing risks models. Many studies have been conducted under semi-competing 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 lost-to-follow-up (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 multi-state model in which LTF is not regarded as censoring but as a non-fatal event. Considerable studies have utilized this multi-state 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 non-fatal 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 non-fatal event can occur. Ultimately, different weights can be assigned to each sub-interval. Thus, the former method corresponds to an unconditional probability approach with equal weight on all of the sub-intervals, whereas the latter utilizes a conditional probability approach with a specific weight, depending on the sub-interval.
In our study, we use the latter method to deal with non-fatal events that are interval-censored on an interval. In addition, we propose an additive-multiplicative model by combining the Cox proportional hazards model with the additive risk model of [16], in accordance with a multi-state model. The additive-multiplicative model was initially introduced by [17 ] and has since been developed by a number of researchers. Scheike and Zhang [19] incorporated time-varying covariates for the additive part and time-independent covariates for the multiplicative part. On the other hand, [18] estimated relevant parameters by considering time-varying 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.

Models
As depicted in Fig. 1, the proposed model in this study consists of five states: healthy (H), non-fatal (NF), fatal (F), lost-to-follow-up (LTF), and unobserved non-fatal (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 (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) / ∈ 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 log-normal 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.

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 non-fatal event, LTF, and a fatal event, respectively. Furthermore, let H 0 (s) be the set of subjects staying in state 0 at time s. That is, Let 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 non-fatal event was observed, and b i be the first time a non-fatal 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 B j = {i : I ij = 1}. For subject i ∈ B 1 ∪ B 2 , we have a i , b i ≥ t i ; this is the case because a non-fatal event has not been observed before time t i . For subject i ∈ B 3 ∪B 4 , we have a i < b i ≤ t i ; this is the case because a non-fatal event has occurred between a i and b i . When subject i is a member of B 5 ∪ 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 non-fatal event has occurred. Thus, t i would be a censoring time for i ∈ B 1 ∪ B 3 ∪ B 5 , whereas it would be a time of death for i ∈ B 2 ∪ B 4 ∪ 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 B 3 ∪ B 4 can occur uniformly on the interval (a i , b i ] [6]. However, we partition the interval (a i , b i ] into several sub-intervals where non-fatal events could occur and assign different weights to each interval [15]. • Let R i ∈ (a i , b i ] be an interval for the occurrence of non-fatal events associated with subjects in routes 3 or 4. Let s 1 be the smallest value among all b i 's for subjects in the set B 3 ∪ B 4 . Let s 2 be the smallest value among all b i 's corresponding to subjects having a i greater than or equal to s 1 . This process is repeated until we have no subjects with a i greater than or equal to s m (m = 1, 2, . . .). Thus, we can have a refined set of time points • We can define the weight w i m at time s m (m = 1, 2, . . .) for subject i in the set B 3 ∪ B 4 : where . Subsequently, likelihood functions incorporated with weight w i 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 (ζ ) = − log m(ζ ), which is referred to asζ . Consequently, the inverse of the Hessian matrix evaluated atζ is defined as the estimated variance-covariance matrix ofζ . Numerical integration is required for the frailty distribution. For this purpose, we use the adaptive importance sampling [20]. Finally, we employ quasi-Newton optimization, which utilizes the gradient vector and the Hessian matrix of f (ζ ), to achieve the optimal solution of ζ .

Simulation studies
Extensive simulation is performed to investigate the finite-sample 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 log-normal 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 non-fatal events to be 24 times in a 12-month 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 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: Step 3. • Step 2: Let u 12i be a random number generated from a uniform distribution on the interval Redefine T i as the root s of the equation, If C ≤ T i , then the i th subject is defined as being censored after experiencing a non-fatal event, i.e., i ∈ B 3 . Otherwise, the i th subject is defined as being dead at time T i after experiencing a non-fatal event, , the type of path should be redefined because a non-fatal 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 non-fatal event, i.e., i ∈ B 1 . Otherwise, the i th subject is defined as being dead at time T i without experiencing a non-fatal event, i.e., i ∈ B 2 . • Step 3: Let u 32i and u 34i be random numbers generated from uniform distributions on the intervals Now redefine R i and T i as the roots s of the equations: respectively. If C ≤ R i ∧ T i , the i th subject is defined as being censored without experiencing a non-fatal event after LTF, i.e., i ∈ B 5 . If T i ≤ R i , then the i th subject is defined as being dead without experiencing a non-fatal event after LTF, i.e., i ∈ 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 Redefine T i as the root s of the equation, If C ≤ T i , then the i th subject is defined as being censored at time C after experiencing LTF and a non-fatal event, i.e., i ∈ B 5 . Otherwise, the i th subject is defined as being dead at time T i after experiencing LTF and a non-fatal event, i.e., i ∈ 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 non-fatal 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 non-fatal event occurred at the end of the right endpoint of the interval ('imputed-by-the-right-endpoint'). 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 'imputed-by-the-right-endpoint' 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 socio-demographic 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 dementia-free), dead or alive status, age (including a n's age at the start of study, their age at the last dementia-free 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 dementia-free 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 p-values for the test results [11]. For the gender effect, only the p-value for transition 0 → 3 is greater than 0.1, which seems to violate the proportional hazards assumption. Alternatively, the p-values 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: λ rs (t|z = gender, w = certificate, u) = exp(u) β rs z + exp(α rs w)θ rs γ rs t γ rs −1 for (r, s) ∈ A. Table 7 shows a summary of the estimation procedures. We provide the estimates of the regression coefficients along with their standard errors and p-values. 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 p-value 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 p-value 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 non-significant result with a p-value 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 p-value 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 p-value of 0.038. Meanwhile, all transition intensities of the non-educated 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 p-value less than 0.001, showing non-homogeneity 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 multi-state model for semi-competing risks data, for which there exists information on fatal events, but information on non-fatal events may not be available due to lost to follow-up. 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 A B C D E Fig. 3 Five transition intensities over age by gender and certificate: 0 → 1, 0 → 2, 1 → 2, 0 → 3, 3 → 2 transitions and estimated normal frailties of each cluster classified by age at entry 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 quasi-Newton 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 non-diagnosed subjects. Furthermore, we conducted simulations with finite-sample 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 semi-parametric 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 Nelson-Aalen estimators for the cumulative baseline transition intensity [8,9,16] and extend this method to semi-competing 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