Modeling absolute differences in life expectancy with a censored skew-normal regression approach

Data

The Swiss National Cohort (SNC) is a longitudinal study of the entire resident population of Switzerland, based on the 1990 and 2000 national censuses (Bopp et al., 2009). Deterministic and probabilistic record linkage (Fellegi & Sunter, 1969) were performed using the Generalized Record Linkage System (Fair, 2004) to link census records to a death record or an emigration record, based on a set of key variables that are available in both data sets (sex, date of birth, place of residence, marital status, religion, nationality, profession, date of birth of partner and date of birth of children). Mortality patterns and life expectancy patterns are of major interest in the SNC (Moser et al., 2014; Spoerri et al., 2014; Spoerri et al., 2006). Initial SNC mortality linkage was successful for 94% of the deaths (Bopp et al., 2009). The 6% not linked deaths bias the calculation of absolute age-specific mortality rates which in turn would bias estimates of life expectancy (Schmidlin et al., 2013). This bias was removed when including the 6% deaths using pragmatic linkages matching for age and sex only. For analyses presented here, we used the SNC data from the 2000 census onwards with mortality follow-up up to end of 2008 including the pragmatically linked deaths. We investigated 4,098,675 individuals aged ≥35 years or older from the 2000 census. Of those, 481,157 individuals died between 5th December 2000 (date of census) and 31th December 2008. We investigated associations of gender, civil status and education on individual’s life span, using parametric survival regression and censored skew-normal and censored Gaussian regression approaches.

For a second analysis and the simulation study we used data from the Human Mortality Database (HMD) which contains death rates and life tables from various populations including Switzerland (Human Mortality Database, 2014). Data are provided from national statistical offices or other sources. We used death rates for one year age intervals from age 35 to 105. For a hypothetical population of N = 100, 000 men and women we estimated the number of deaths for each one year age interval from age 35 to age 105, I_[x,x+1] ≔ I_x, x ∈ 35, …, 105, as follows. The death rate for an age interval I_x, x ∈ 35, …, 105, is denoted as m_x. The number of persons alive at the start of the age interval I_x is denoted as l_x, x ∈ 36, …, 105. l_x is equal to 100,000 minus all the deaths in all age intervals before I_x. The number of deaths n_x for each one year age interval I_x was then calculated as n_x = m_x ⋅ l_x, x ∈ 35, …, 105. Data were downloaded and analyzed using the R-package demography (Hyndman, 2012).

Parametric survival models

Parametric survival regression assumes that, for each individual, the time T elapsed from a starting time point (e.g., at age 35) to an event (e.g., death) has a cumulative distribution function (cdf) F(t) and a probability density function (pdf) f(t) from certain classes of distribution functions. Often one assumes a distribution function from an extreme value type, e.g., Weibull, or from a lifetime type, e.g., Gompertz. The survival function is defined as S(t) = 1 − F(t) = ℙ(T > t).

One basic concept of survival regression is the hazard function, defined as (1)${\displaystyle \lambda \left(t\right)=\underset{u\to 0}{lim}\frac{\mathbb{P}\left(t\le T<t+u|T\ge t\right)}{u}=\frac{f\left(t\right)}{1-F\left(t\right)}=\frac{f\left(t\right)}{S\left(t\right)}.}$ The survival function can then be expressed in terms of the hazard function as ${\displaystyle S\left(t\right)=exp\left(-{\int }_0^t\lambda \left(y\right)dy\right).}$

If one assumes that the observed survival time of each individual is a realization from the same distribution of T (without covariates representing differences between the groups of individuals) and T comes from a Weibull distribution function, then ${\displaystyle \lambda \left(t\right)=\frac{\gamma }{\alpha }{\left(\frac{t}{\alpha }\right)}^{\gamma -1},S\left(t\right)=exp\left(-{\left(\frac{t}{\alpha }\right)}^{\gamma }\right),t>0,}$ where α > 0 is a scale parameter and γ ≥ 0 a shape parameter of the Weibull distribution, or when T is Gompertz distributed with scale parameter α > 0 and shape parameter γ > 0, then ${\displaystyle \lambda \left(t\right)=\alpha exp\left(\gamma t\right),S\left(t\right)=exp\left(-\frac{\alpha }{\gamma }\left[exp\left(\gamma t\right)-1\right]\right),t>0.}$

For the survival regression problem one often assumes that for a set of covariates X = (X₁, …, X_k)^⊤ the following equation holds ${\displaystyle \lambda \left(t|\mathbf{X}\right)=\lambda \left(t\right)exp\left({\mathbf{X}}^{\top }\mathbf{\beta }\right),}$ where β = (β₁, …, β_k)^⊤ is a vector of regression coefficients. This model specification is known as proportional hazard (PH) model and λ(t) is then often called an underlying baseline hazard function. In case that T is assumed to be Weibull distributed one gets ${\displaystyle \lambda \left(t|\mathbf{X}\right)=\frac{\gamma }{\alpha }{\left(\frac{t}{\alpha }\right)}^{\gamma -1}exp\left({\mathbf{X}}^{\top }\mathbf{\beta }\right),}$ or with the assumption of a Gompertz distribution ${\displaystyle \lambda \left(t|\mathbf{X}\right)=\alpha exp\left(\gamma t\right)exp\left({\mathbf{X}}^{\top }\mathbf{\beta }\right).}$ Note that in such a model representation the relationship between the covariates and the hazard function are linear on the log hazard scale. Thus, the effect of increasing a continuous covariate, say X₁, by d, holding all other variables constant, is to increase the hazard of the event by a factor of exp(β₁d) at all time point in time, assuming X₁ is linearly related to logλ(t) (Harrell, 2001). Often the above described effect is reported in terms of hazard ratios, i.e., one compares the ratio of hazard rates of an individual with predictor value d compared to one with a value of 0. The hazard ratio is then $\lambda \left(t\right)exp\left({\beta }_{1}d\right)/\lambda \left(t\right)=exp\left({\beta }_{1}d\right)$. Unlike the interpretation of a multivariable linear Gaussian regression model, where the regression coefficients β are reflecting an increment in the expected value of a response variable by a one unit change in the predictor variable, the interpretation of regression coefficients from survival regression models is not easily translated to differences in mean survival time.

Another often used type of survival model is the accelerated failure time (AFT) model (see e.g., Harrell (2001)), where one assumes that (2)${\displaystyle log\left(T\right)={\mathbf{X}}^{\top }\mathbf{\beta }+\sigma ϵ,}$ with σ, a scale parameter, and ϵ is assumed to come from a survival distribution function ϑ.

Common choices of ϑ(u) are the logistic distribution ϑ(u) = [1 + exp(u)]⁻¹ (log-logistic model) or the Gaussian distribution ϑ(u) = 1 − Φ(u) (log-normal model). Both distributions fail to fit lifetime data adequately, because of their positive skewed distribution. To address a negative skewed distribution, an underlying Weibull distribution is possible. Note from (2) that in the AFT model specification the interpretation of the estimated parameters are in terms of T = exp(X^⊤β + σϵ). The effect of increasing a continuous covariate, say X₁, by d, holding all other variables constant, is to increase the survival time T by a factor of exp(β₁d), assuming Eq. (2). Similarly to reporting hazard ratios in PH models, one reports time ratios in AFT models, i.e., the ratio of survival times of an individual with predictor value d compared to one with a value of 0. The time ratio is then exp(β₁d + σϵ)/exp(σϵ) = exp(β₁d). Also in this type of survival modeling, interpretation of regression coefficients is not straightforward in terms of mean survival time. Note that the Weibull PH model and the Weibull AFT are equivalent (Harrell, 2001).

Life expectancy from survival models

It is well known that life expectancy (or expected survival time) conditional on covariates 𝔼(T|X) is related to the conditional survival function S(t|X) through ${\displaystyle \mathbb{E}\left(T|\mathbf{X}\right)={\int }_0^{\infty }S\left(t|\mathbf{X}\right)dt={\int }_0^{\infty }{\left\{S\left(t\right)\right\}}^{exp\left({\mathbf{X}}^{\top }\mathbf{\beta }\right)}dt,}$ where S(t) is the underlying survival distribution, see e.g., Harrell (2001). Hence, the expected survival time for a reference individual is $\mathbb{E}\left(T|\mathbf{X}\right)={\int }_0^{\infty }S\left(t\right)dt$. For Weibull, Gompertz and log-normal distribution functions closed-form expressions exist, i.e., for the Weibull distribution ${\displaystyle \mathbb{E}\left(T|\mathbf{X}\right)=\alpha \Gamma \left(1+1/\gamma \right),}$ for the Gompertz distribution ${\displaystyle \mathbb{E}\left(T|\mathbf{X}\right)=\frac{1}{\gamma }exp\frac{\alpha }{\gamma }{\int }_{\alpha /\gamma }^{\infty }{x}^{-1}exp\left(-x\right)dx,}$ and for the log-normal distribution ${\displaystyle \mathbb{E}\left(T|\mathbf{X}\right)=exp\left(\mu +{\alpha }^2/2\right),}$ where μ, α, γ are location, scale and shape parameters from the corresponding distribution functions, see e.g., Johnson, Kotz & Balakrishnan (1994) and Missov & Lenart (2011). 95% confidence intervals (CI) were approximated by sampling procedures from multivariate Gaussian vectors using the covariance matrices of the parameter estimates.

Censored skew-normal regression

The skew-normal distribution function generalizes the Gaussian distribution function allowing for skewness in its shape. We start with a recapitulation of the definition of a Gaussian distributed random variable. A random variable X is Gaussian (normal) distributed with location parameter μ ∈ ℝ and scale parameter σ > 0 if it has the pdf ${\displaystyle f\left(x;\mu ,\sigma \right)=\frac{1}{\sigma \sqrt{2\pi }}{exp}^{-\frac{1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2}.}$ One writes X ∼ N(μ, σ²) if X is Gaussian distributed with mean μ and variance σ². The pdf of standard Gaussian distributed random variable Z ∼ N(0, 1) is in the following written as φ(⋅).

The definition of a skew-normal distributed random variable is as follows, see e.g., Genton (2004).

The skew-normal regression problem is similar to Gaussian linear regression. One assumes for a set of covariates X = (X₁, …, X_k)^⊤ and regression coefficients β = (β₁, …, β_k)^⊤ that ${\displaystyle Y={\mathbf{X}}^{\top }\mathbf{\beta }+ϵ,}$ where ϵ ∼ SN(0, σ², ψ) and is assumed to be independent across individuals. The effect of increasing a continuous covariate, say X₁, by d, holding all other variables constant, is to increase the survival time T by a shift of β₁d. For the conversion of DP to CP (or vice versa) in the regression problem only the distributional parameters need to be transformed accordingly (Azzalini & Capitanio, 1999).

Parameter estimation

Parameters for a censored skew-normal regression can be estimated by maximum likelihood estimation. Let the underlying measurement scale be individual’s age at end of study Y_i or censoring C_i, i.e., min{Y_i, C_i}. Censoring is assumed to be non-informative (i.e., censoring is independent of the underlying time scale), or from type I censoring (i.e., follow-up time ends at predetermined time). We write δ_i = 𝕀(Y_i ≤ C_i) for the indicator variable for an observed death, where 𝕀 is the indicator function. Under censoring the likelihood is (6)${\displaystyle L\left(\xi ,{\sigma }^2,\psi \right)=\prod _{i\le n:\underset{i}{\delta }=1}g\left(y_i;\xi ,{\sigma }^2,\psi \right)\prod _{i\le n:\underset{i}{\delta }=0}S\left(c_i;\xi ,{\sigma }^2,\psi \right)=\prod _{i\le n:\underset{i}{\delta }=1}\frac{2}{\sigma }\phi \left(\frac{y_i-\xi }{\sigma }\right)\Phi \left(\psi \frac{y_i-\xi }{\sigma }\right)\prod _{i\le n:\underset{i}{\delta }=0}S\left(c_i;\xi ,{\sigma }^2,\psi \right)=\prod _{i\le n}{\left[\frac{2}{\sigma }\phi \left(\frac{y_i-\xi }{\sigma }\right)\Phi \left(\psi \frac{y_i-\xi }{\sigma }\right)\right]}^{{\delta }_i}{\left[S\left(c_i;\xi ,{\sigma }^2,\psi \right)\right]}^{1-{\delta }_i},}$ where g(⋅; ξ, σ², ψ) is the pdf given in (3) and S(⋅; ξ, σ², ψ) is the survival function assuming F(⋅) = SN(⋅; ξ, σ², ψ). Choosing the relevant time scale is a crucial decision in modeling survival times. Often not the total risk time starting from time origin (e.g., date of birth) is observed, but only the risk time from study entry (the date at which a person entered the study and came under observation) until the end of the study. This concept is called delayed entry or left-truncation. In this case one has to consider the conditional distribution of Y_i given that Y_i ≥ D_i, where D_i is a given time point or individual’s age at study entry. Note that the likelihood of an uncensored individual i ≤ n : δ_i = 1 is then ${\displaystyle g^{\ast }\left(y_i;\xi ,\sigma ,\psi |Y_i\ge D_i\right)=\frac{g\left(y_i;\xi ,{\sigma }^2,\psi \right)}{S\left(D_i;\xi ,{\sigma }^2,\psi |Y_i\ge D_i\right)},-\infty <y_i<\infty ,i=1\le n,}$ and for a censored individual i ≤ n : δ_i = 0, ${\displaystyle S^{\ast }\left(c_i;\xi ,\sigma ,\psi |Y_i\ge D_i\right)=\frac{S\left(c_i;\xi ,{\sigma }^2,\psi \right)}{S\left(D_i;\xi ,{\sigma }^2,\psi |Y_i\ge D_i\right)},-\infty <c_i<\infty ,i=1\le n.}$ To obtain the likelihood of all individuals one replaces in (6) g(y_i; ξ, σ², ψ) by g^∗(y_i; ξ, σ, ψ|Y_i ≥ D_i), and S(c_i; ξ, σ, ψ) by S^∗(c_i; ξ, σ, ψ|Y_i ≥ D_i), respectively. It has been mentioned in e.g., (Azzalini, 1985) that maximizing the negative log likelihood of (6) has singularity problems if ψ = 0, and yield convergence problems in the MLE. To overcome this problem it has been suggested to use the CP approach, which removes the singularity at ψ = 0, and has further advantages in faster convergence and improved likelihood shape over the DP (see e.g., Azzalini (1985) and Azzalini & Capitanio (1999)). We used the CP for the numerical derivation of the estimates by MLE using R Version 3.1.1 and Stata Version 13.1. Program code and used functions are provided as Supplemental Information 1.

Model assessment

For all types of survival models, the model should adequately fit the data in order to obtain correctly interpretable estimates and correct coverage of confidence intervals. For PH models the relationship between the covariates and the log hazard should be linear. Further, the covariates affect the underlying distribution of the time variable by adding X^⊤β to the log hazard function. The effect of the covariates is assumed to be the same at all time points. For AFT models each covariate affects log(T) linearly. Further, the underlying variance σ in Eq. (2) is a constant, independent of the covariates (Harrell, 2001). Assessing the model fit of parametric survival models is often restricted to e.g., graphical assessment by stratified predictor levels or stratified Kaplan–Meier estimates of the distribution of residuals. For the skew-normal regression problem Y = X^⊤β + ϵ, where ϵ ∼ SN(0, σ², ψ), certain assumptions must be validated (Harrell, 2001): residuals should have no systematic trend in central tendency against any predictor, they should have the same dispersion and they should have a skew-normal distribution in the predictor-space. These assumptions can be checked by median and lower and upper quartiles of the residuals, stratified by intervals of the predicted outcome. Distributional model assumptions of a skew-normal regression can be visually checked by a comparison of quantiles of estimated residuals and quantiles from a theoretical skew-normal distribution in a quantile–quantile plot.

We graphically compared the model fit from a Gompertz model with a skew-normal model by a log-hazard plot. In general from Eq. (1) one obtains that logλ(t) = logf(t) − logS(t). Note that for a Gompertz distribution the log-hazard is linear in the time scale t, i.e., logλ(t) = logα + γt, t ≥ 0. Thus, any non-linearity in the log-hazard from a skew-normal model would indicate a deviation from the Gompertz model fit.

Goodness-of-fit using official life expectancy estimates was assessed by Pearson’s chi-squared test statistic $X^2=\sum _i{\left(O_i-E_i\right)}^2/E_i$ (Agresti, 2013), where O_i denotes the observed number of deaths from offical estimates in one year age intervals. We calculated expected number of deaths E_i for one year age intervals from the underlying regression models. A larger value of X² indicates a greater difference of O_i − E_i (Agresti, 2013).

Model setting

SNC data were analyzed by survival regression models with remaining age at 35 years as the underlying time scale and delayed entry date as the 5th of December 2000. Assumed underlying survival distribution function was either Weibull or Gompertz for a PH model, or Weibull and log-normal for an AFT model. Individuals were censored if they were alive after 31th December 2008. For a direct comparison using the censored skew-normal regression approach we investigated age at 31th December 2008 (censoring information) or age at death as the dependent variable with a delayed entry date of 5th December 2000. To compare with a Gaussian linear regression approach we used an author programmed MLE function for censored Gaussian regression with left-truncated observations, not further described here. HMD data were analyzed using the same regression models but without delayed entry or censoring.

Simulation study for model distribution misspecification

Parametric modeling assumes that a given underlying distribution function is the true distribution of the outcome variable. We performed a simulation study to assess whether life expectancy estimates of a skew-normal regression and a Gompertz survival regression are biased in case of model distribution misspecification. For that purpose we first fitted Gompertz and skew-normal models to Swiss data of the Human Mortality Database to get location, scale and shape parameters for each distribution, as close to real data as possible. Second, using these parameter estimates, we built random samples with different samples sizes (i.e., 100, 1,000, and 10,000) from Gompertz and skew-normal distribution functions. Third, for each sample we fitted a Gompertz or a skew-normal model and reported estimated life expectancy $\hat{\mu }$ and corresponding confidence intervals. As a third distribution function we combined the Gompertz and skew-normal distribution functions to get a mixture distribution functions with mixing proportions δ = {0.1, 0.25, 0.5, 0.75, 0.9}, where the mixing pdf is defined as m(x) = δf_G(x) + (1 − δ) g(x) with f_G, the Gompertz pdf, and g(x) the skew-normal pdf, defined as in Eq. (3). Thus, with the mixture distribution we mimic a situation where the study sample consists of samples from different underlying distribution functions, with different mixing proportions. For each sample we reported coverage of the true mean life expectancy μ₀ and the bias $\hat{\mu }-{\mu }_0$. We did 1,000 simulation repetitions.