2.1. A three-state time-inhomogeneous Markov process

Following previous literature (see, e.g. Olivieri & Pitacco, Reference Olivieri and Pitacco2001; Rickayzen & Walsh, Reference Rickayzen and Walsh2002; Fong et al., Reference Fong, Shao and Sherris2015; Shao et al., Reference Shao, Sherris and Fong2017), we assume that individuals’ health transitions can be modelled as a multi-state Markov process, where the conditional probability distribution of future states of the process (conditional on both past and present values) only depends on the state presently occupied and is independent of the process history. We define a three-state Markov process as in Figure 1. The process has two transient states, “N” (non-disabled) and “F” (functionally disabled), and one absorbing state, “D” (dead). It allows for three transitionsFootnote ^{1 ,} Footnote ² :

• ∙ σ: N→F, the intensity to become functionally disabled.

• ∙ μ: N→D, the mortality intensity for a healthy person.

• ∙ ν: F→D, the mortality intensity for a disabled person.

Figure 1 Three-state Markov process.

We assume that health transitions follow a time-inhomogeneous Markov process, where the transition probability depends on the time at which the transition takes place:

(1) $$P_{{i,j}} (x,t,h){\equals}{\rm Pr}(S(x{\plus}h,t{\plus}h){\equals}j\,\mid\,S(x,t){\equals}i)$$

(2) $$\alpha _{{i,j}} (x,t){\equals}\mathop{{{\rm lim}}}\limits_{{h\to0^{{\plus}} }} P_{{i,j}} (x,t,h)\,/\,h$$

where x represents age, t the time with h≥0. S(x,t) denotes the stochastic health status of an individual at age x and time t, and i,j∈{N,F,D}. P _i,j (x,t,h) denotes the transition probability from state i at age x and time t to state j at age x+h and time t+h. α _i,j (x,t) denotes the corresponding transition intensity at age x and time t.

2.2. Model specification

Following earlier works of Renshaw & Haberman (Reference Renshaw and Haberman1995) and Fong et al. (Reference Fong, Shao and Sherris2015), we model the transition intensities using a GLM approach. Separate GLMs are estimated for each of the three transition intensities σ, μ and ν. The models are specified by three components: the link function, the linear predictor and the probability distribution.

Link function: We adopt a log link function g(⋅):

(3) $$g(\alpha _{{i,j}} (x,t)){\equals}{\rm ln}(\alpha _{{i,j}} (x,t)){\equals}\eta _{{x,t}} $$

where α _i,j (x,t) are the respective transition intensities σ _x,t , μ _x,t or ν _x,t for age x at time t. η _x,t is a linear predictor of regressors.

Linear predictor: As our primary interest is to explore time trends in health transitions, we introduce time effects as additional covariates besides the age factors considered in Fong et al. (Reference Fong, Shao and Sherris2015) and allow for age–time interactionsFootnote ³ . The linear predictor is given by:

(4) $$\eta _{{x,t}} {\equals}\beta _{0} {\plus}\beta _{1} x{\plus}\beta _{2} x^{2} {\plus}\beta _{3} t{\plus}\beta _{4} tx{\plus}\beta _{5} tx^{2} $$

where x represents age, t the time and the β _j are unknown coefficients that need to be estimated. The model allows for some of the values of β _j to be 0, allowing for flexibility in the functional form.

We include age factors up to the quadratic effect in agreement with the findings of Fong et al. (Reference Fong, Shao and Sherris2015). This is also in line with common practice in mortality modelling (see, e.g. Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009). Since the CLHLS data only allow us to compute at most five transition intensities per individual, we focus on a linear time trend in this study. The inclusion of age–time interaction effects has been an important feature in recent developments in mortality modelling (Cairns et al., Reference Cairns, Blake and Dowd2006; Li et al., Reference Li, O’Hare and Zhang2016; Plat, Reference Plat2009). It ensures that the improvement in mortality has a non-trivial correlation structure across different age groups. Moreover, several studies have recognised the benefits of including quadratic age–time effects for model fitting and forecasting (Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009; Dowd et al., Reference Dowd, Cairns, Blake, Coughlan, Epstein and Khalaf-Allah2010). Therefore, to build on recent developments in mortality and health transition modelling, and to keep the model parsimonious and interpretable, we consider the aforementioned age, time and age–time factors in this study.

Probability distribution: The Poisson distribution has been widely used in the actuarial literature to model the number of deaths (see, e.g. Brouhns et al., Reference Brouhns, Denuit and Vermunt2002; Cairns et al., Reference Cairns, Blake and Dowd2006; Haberman & Renshaw, Reference Haberman and Renshaw2009). In this paper, assuming that each transition intensity described in section 2.1 is constant for each 1-year age group in a given time interval, we assume that the number of health transitions follows an independent Poisson distribution between survey waves. For illustrative purpose, we use the mortality rate μ _x,t from the healthy state as an example in the following. Let n _x,t be the number of transitions from state N to D at age x and time t:

(5) $$n_{{x,t}} \,\sim\,{\rm Poisson}\,(e_{{x,t}}^{H} {\rm }\mu _{{x,t}} ){\rm }\forall {\rm }x,t$$

where $e_{{x,t}}^{H} $ is the central exposure to risk in healthy state at age x and time t.

The Poisson assumption implies that the dispersion parameter equals 1, which means that the mean and variance of transition counts should be the same. However, several recent mortality studies have found that death data has an “overdispersed” feature in many countries (see, e.g. Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009; Li et al., Reference Li, O’Hare and Vahid2015), implying that the variance of the number of deaths is much larger than the mean. Heterogeneity is often considered as a key reason for overdispersion in insurance data. Moreover, strong dependence in the observed sample can also lead to overdispersion. Therefore, we test the dispersion parameter in each of the transition counts before estimating the GLMs. We found that in most of the cases the dispersion parameter is close to 1Footnote ⁴ . Only for cases where the dispersion parameter is significantly different from one, we relax the restriction on the value of the dispersion parameter and estimate this parameter based on the underlying data. However, the estimates for parameters in the model will remain unchanged upon the introduction of an dispersion parameter in the Poisson distribution (for discussions, see Hinde & Demetrio, Reference Hinde and Demetrio1998; McCullagh & Nelder, Reference McCullagh and Nelder1989).

2.3. Estimation and model selection

Maximum likelihood estimation (MLE) is used to obtain estimates of the proposed GLM models. We define Φ as the parameter set. Using the mortality rates of healthy individuals again as an example, the log-likelihood function is given by:

(6) $$l({\rm \Phi }\,{\hbox {;}}\,n,e){\equals}\mathop \sum\limits_x \mathop \sum\limits_t \left\{ {n_{{x,t}} {\rm ln}[e_{{x,t}}^{H} \,\mu _{{x,t}} ({\rm \Phi })]{\minus}e_{{x,t}}^{H} \,\mu _{{x,t}} ({\rm \Phi }){\minus}{\rm ln}(n_{{x,t}} \,!\,)} \right\}$$

We select the model specification for each transition intensity by comparing the Bayesian information criterion (BIC) for all possible combinations of the six terms in equation (4). That is, we do not impose a single model structure for all health transition rates. Our modelling approach is data-driven and more flexible than relevant earlier works (e.g. Haberman & Renshaw, Reference Haberman and Renshaw2009; Fong et al., Reference Fong, Shao and Sherris2015). For each type of transition, we tailor the model design to only include those age, time and age–time effects that are important and relevant.

We choose the BIC for model selection because it is widely used in statistics and is proven to be consistent (Schwarz, Reference Schwarz1978). The BIC penalises the number of parameters estimated in the model as follows:

(7) $${\rm BIC}{\equals}{\minus}2{\rm }l({\rm \hat{\Phi }}){\plus}k\,{\rm ln}(N)$$

where l( $$\hat{\Phi }$$ ) is the log-likelihood based on the MLE estimators, N the total number of observations and k the number of parameters in the model. The model with the smallest BIC value is selected as the preferred model. We also provide the results of a stepwise comparison of several nested model variants based on the BIC in section 4.2.

4.1. Selected models: estimation results

We estimated the GLM described in section 2.2 separately for the three transition intensities σ, μ and ν for each sub-population in our sample. Table 3 gives the estimation results of the selected linear predictor for each case based on a comparison of all possible model variants as described in section 2.2. The selected models all have highly significant parameters and the results are interpretable.

Table 3 Selected model: parameter estimates.

Note: The functional form of the linear predictor is η _x,t =β ₀+β ₁ x+β ₂ x ²+β ₃ t+β ₄ tx+β ₅ tx ².

BIC=Bayesian information criterion.

*p<0.05; **p<0.01; ***p<0.001.

The selected models for the disability rate σ _x,t in each subpopulation all include a positive linear and a negative quadratic age terms, implying that disability rates increase with age, but at a decreasing rate. Moreover, the negative age–time interaction effects in all four subpopulations except for urban males show that there has been an overall improvement in disability rates over time. It also shows that the rate of improvement in disability rates over time differs across age groups.

For all subpopulations, the mortality rates for healthy individuals μ _x,t increase with age but again there is a deceleration for higher age groups. We note that there are no significant time effects or age–time interaction effects in any of the four sub-populations. This agrees with the fact that the plots of μ _x,t show fairly stables pattern throughout the sample period (see Figure 3).

The models for the mortality rates of the disabled ν _x,t all include positive linear age effects and negative time/age–time effects. For urban and rural males and urban females, a linear negative mortality trend is found for all age groups. The speed of mortality decline over time is similar in these three sub-populations. The model for rural females includes a negative quadratic age–time effect, indicating that mortality decline for this subpopulation is more rapid for older age groups.

Overall, these results show that both age and time effects are important factors explaining patterns in health transitions at higher ages in China. In addition, several of the models rely on age–time interactions which capture time trends that differ by age. Our results also suggest that the improvements in mortality rates of older Chinese aged 65+ are largely driven by the decline of mortality rates of functionally disabled elderly, rather than by the mortality rates of healthy individuals. China has experienced rapid economic growth since the beginning of its market reforms in 1978. Over our sample period 1998–2012, the average annual disposable income increased 3.5 times for urban households and 2.6 times for rural households (National Bureau of Statistics of China, 2013). This economic development allowed households to afford higher living standards and better medical care. Furthermore, several policy reforms have improved people’s access to health services and increased health insurance coverage in China (Meng et al., Reference Meng, Xu, Zhang, Qian, Cai, Xin, Gao, Xu, Boerma and Barber2012). As a result, China has seen major changes in the causes of death (Zhou et al., Reference Zhou, Wang, Zhu, Chen, Wang, Liu, Li, Wang, Liu, Yin and Liu2016) and in the survival and disability of the elderly (Zeng et al., Reference Zeng, Feng, Hesketh, Christensen and Vaupel2017). We argue that these factors can explain our results.

Figures A.1–A.3 in the Appendix show the residuals for the 12 selected models, computed as the difference between the crude and estimated transition rates. The errors fluctuate around 0 and show no systematic patterns. We conclude that the selected models effectively capture age and time patterns in the data.

4.2. Stepwise model selection

The previous section discussed the selected model specifications for each transition intensity and each subpopulation that were identified by comparing all possible model variants. It is interesting to compare these results with those from a stepwise model selection process where additional terms are added in each step. Table 4 gives the BIC values for six nested model variants and compares these with the BIC values for the selected models identified in Table 3.

Table 4 Stepwise model selection: goodness-of-fit of nested models.

Note 1: The table gives the Bayesian information criterion (BIC) for several nested model variants.

Note 2: Bold font indicates minimum BIC values. Selected model refers to the model identified in Table 3.

We note that in six of the 12 cases the stepwise model selection identifies models with higher BIC values (representing a worse model fit) than the models selected in section 4.1. This is the case for all four models for the mortality rate from functionally disabled state ν _x,t , and for two of the models for the disability rate σ _x,t . The limitations of stepwise model selection algorithm are widely recognised in statistics (Hurvich & Tsai, Reference Hurvich and Tsai1990; Grafen et al., Reference Grafen, Hails, Hails and Hails2002; Whittingham et al., Reference Whittingham, Stephens, Bradbury and Freckleton2006). One of the weaknesses of this method is the fact that model selection is very sensitive to factors such as the order of parameter entry and whether we choose to use forward selection algorithm or backward elimination algorithm (Derksen & Keselman, Reference Derksen and Keselman1992). Therefore, to avoid these limitations, in this paper we have considered all possible model designs for the three health transition models. The preferred models turn out to be parsimonious and have good fitting performances.

Nevertheless, the detailed analysis of the stepwise model comparison confirms that including time effects and age–time interaction terms improves the model fit for most health transition models except for the mortality rates μ _x,t of the non-disabled.

4.3. Likelihood ratio tests (LRTs)

In order to further verify the significance of the time effects in the three transition intensities, we conduct the LRT on the selected models given in Table 3. For cases where time effects are included in the selected models, we drop these effects to construct the null model and then perform the LRT. For cases where the selected model does not contain time effects, we construct an alternative model with all the selected age terms and at least one time effect. We select the alternative model using the BIC. The test statistic based on the deviance D is defined as:

(8) $${\rm \Delta }D{\equals}D_{{null}} {\minus}D_{{alternative}} {\equals}{\minus}2{\times}{\rm log}(L_{{null}} \,/\,L_{{alternative}} )$$

with

(9) $$D_{{null}} {\equals}{\minus}2{\times}{\rm log}(L_{{null}} \,/\,L_{{saturated}} )$$

(10) $$D_{{alternative}} {\equals}{\minus}2{\times}{\rm log}(L_{{alternative}} \,/\,L_{{saturated}} )$$

where L _null represents the likelihood of the null (simpler) model and L _alternative the likelihood of the alternative (larger) model. The saturated model has the same number of parameters as the sample size and L _saturated is the corresponding likelihood.

Under the null hypothesis that the simpler model is preferred, ΔD follows a χ ² distribution with degrees of freedom equal to the difference in the number of parameters between the null model and the alternative model. We conduct the LRT at the 5% level of significance. The test results given in Table 5 show that for all selected models which include time effect, adding these time effects significantly improves the model fit. The results also show that adding time effects does not significantly improve the fit for the selected models without time effect. The only exception is the model for the mortality rate μ of healthy urban females, where the LRT indicates that including time effects would improve the model fit. We note this result but decide to keep the more parsimonious model that was selected in Table 3 based on the BIC. Apart from this exception, the results of the LRTs confirm our model choices based on the BIC.

Table 5 Results of the likelihood ratio test.

Note: D denotes the deviance and ΔD is the test statistic defined in section 4.3.

4.4. Life expectancy and healthy life expectancy

We use the models identified in section 4.1 to compute estimates for life expectancy and healthy life expectancy. Table 6 shows the estimated life expectancies at age 65 and 75 conditional on the initial health status, where “healthy” is defined as having at most one ADL limitation (see section 3.1). We provide estimates for the first and last time point of the investigation period (1998 and 2011) and out-of-sample forecasts for 2020. A death time dispersion measure for each subgroup, which is calculated as the standard deviation of death time distributions, is also shown in Table 6. This dispersion measure captures idiosyncratic mortality risk. As death time distribution for older age groups can be very different from a normal distribution, this standard deviation-based measure is an approximate indicator of the death time dispersions around life expectancy.

Table 6 Life expectancy and death time dispersion conditional on health status at age 65 and 75.

The estimated life expectancies vary in plausible ways: life expectancies of urban residents are higher than those of rural residents, females have higher life expectancies than males and healthy individuals have higher life expectancies than disabled ones. Our models include time trends in three of the four models for disability rate and in all models for the disabled mortality rates. These trends are reflected in the life expectancies which increase over time for all population subgroups and show the largest improvements for disabled individuals. When comparing the computed life expectancies in Table 6 with several related studies, we find consistencies in the results. For example, Luo et al. (Reference Luo, Wong, Lum, Luo, Gong and Kendig2016) report for the age group 65–69 in 2011 a remaining life expectancy of 15.0 years for males and 18.7 years for females (Luo et al., Reference Luo, Wong, Lum, Luo, Gong and Kendig2016, Table 2). From Table 6, we can also see that death time dispersions for the disabled subgroups increase over time, a trend which is potentially driven by the mortality improvement of the disabled individuals. On the other hand, the death time dispersions of the healthy subgroups are generally stable over time.

Table 7 gives the estimated healthy life expectancies at age 65 and 75Footnote ⁸ . Female urban residents have the highest healthy life expectancy and male rural residents have the lowest healthy life expectancy. Healthy life expectancies of females improve faster over time than those of males. We find that the ratios of healthy life expectancy to life expectancy are quite stable over the period 1998–2020, indicating a dynamic equilibrium where both life expectancy and healthy life expectancy shift to the right – a finding which agrees with the results of several related studies on China (see, e.g. Liu et al., Reference Liu, Chen, Song, Chi and Zheng2009; Guo, Reference Guo2017).

Table 7 Healthy life expectancy at age 65 and 75.

Overall, our results show persistent health differences between urban and rural China. For life expectancy, we find that the existing urban–rural gaps increase over time for healthy males and for healthy and disabled females. For disabled males, the gap seems to be slowly decreasing (see Table 6). For healthy life expectancy, our results suggest convergence between urban and rural males, but divergence for females.