Flexible multi-state models for interval-censored data: specification, estimation, and an application to ageing research

TL;DR

This paper introduces a flexible estimation method for multi-state survival models with interval-censored data, combining parametric and semi-parametric hazard specifications, demonstrated through ageing research data.

Contribution

It develops a novel scoring algorithm for estimating multi-state models with flexible hazard shapes, including semi-parametric P-spline approaches.

Findings

Effective estimation of multi-state models with interval censoring.

Application to ageing data demonstrates practical utility.

Flexible hazard modeling improves understanding of health transitions.

Abstract

Continuous-time multi-state survival models can be used to describe health-related processes over time. In the presence of interval-censored times for transitions between the living states, the likelihood is constructed using transition probabilities. Models can be specified using parametric or semi-parametric shapes for the hazards. Semi-parametric hazards can be fitted using $P$-splines and penalised maximum likelihood estimation. This paper presents a method to estimate flexible multi-state models which allows for parametric and semi-parametric hazard specifications. The estimation is based on a scoring algorithm. The method is illustrated with data from the English Longitudinal Study of Ageing.

Tables

Table 1. Table 1: State table for the ELSA data: number of times each pair of states was observed at successive observation times. The four living states are defined by number of words remembered

	To
From	10-7 words	6-5 words	4-2 words	1-0 words	Dead
10-7 words	164	150	49	12	8
6-5 words	156	440	303	48	40
4-2 words	52	336	616	151	85
1-0 words	11	35	114	149	72

Table 2. Table 2: Comparison between models for the ELSA data with N = 1000 𝑁 1000 N=1000 , where -2LL stands for -2 times the (penalised) loglikelihood function evaluated at its maximum

Model	Baseline hazards	#Parameters	-2LL	AIC
Intercept-only	Exponential	10	8089.5	8109.5
$t$	Gompertz	15	7750.5	7780.5
$t$, $sex$, $education$	Gompertz	22	7636.3	7680.3
$t$, $sex$, $education$	Gompertz for living	22	7644.7	7688.7
	and Weibull for death
$t$, $sex$, $education$	Weibull	22	7685.5	7729.5
$t$, $sex$, $education$	Weibull for living	22	7675.7	7719.7
	and Gompertz for death
$t$, $sex$, $education$	$P$-splines I	38	7626.0	7678.2
	for $2\to 3$ and $3\to 4$
$t$, $sex$, $education$	$P$-splines II for $3\to 4$	30	7630.9	7678.2

Table 3. Table 3: Results for sex, education and time for the five-state P 𝑃 P -splines II model for the ELSA data. Estimated standard errors in parentheses

$sex$		$education$		$t$
${\beta }_{12.1}$	0.552  (0.138)	${\beta }_{12.2}$	-0.281  (0.146)	${\xi }_{12}$	0.030  (0.010)
${\beta }_{23.1}$	0.178  (0.101)	${\beta }_{23.2}$	-0.836  (0.103)	${\xi }_B$	-0.031  (0.006)
${\beta }_{34.1}$	0.141  (0.145)	${\beta }_{34.2}$	-0.445  (0.160)	${\xi }_D$	0.042  (0.009)
${\beta }_D$	0.477  (0.151)

Equations

$$p_{rs}(t)=P\big{(}Y(t)=s|Y(0)=r\big{)},$$

$$q_{rs}=\lim_{\Delta\downarrow 0}\frac{P\big{(}Y(t+\Delta)=s|Y(t)=r\big{)}}{\Delta},$$

$$\displaystyle q_{rs}(t)=q_{rs.0}(t)\exp\big{(}\bm{\beta}_{rs}^{\top}\mathbf{x}\big{)},$$

$$\text{mbox}Weibull:q_{rs.0}(t)$$

$$\text{mbox}Gompertz:q_{rs.0}(t)$$

$$q_{rs.0}(t)=\exp \left(\sum _{k=1}^K{\alpha }_{rs.k}{B}_k\right),$$

$$P\left(Y_n=y_n,...,Y_2=y_2|Y_1=y_1,\mathbf{\theta },\mathbf{t},\mathbf{X}\right),$$

$$\prod _{j=2}^{n}P\left(Y_j=y_j|Y_{j-1}=y_{j-1},\mathbf{\theta },t_{j-1},{\mathbf{x}}_{j-1}\right),$$

$$\displaystyle L_{ij}=\left\{\begin{array}[]{lll}P\left(Y_{j}=y_{j}|Y_{j-1}=y_{j-1},{\bm{\theta}},t_{j-1},\mathbf{x}_{j-1}\right)&\mbox{for}&j=2,...,n-1\\ C(y_{n}|y_{{n}-1})&\mbox{for}&j=n\ .\end{array}\right.$$

$$C(y_n|y_{n-1})=\sum _{s=1}^{D-1}P\left(Y_n=s|Y_{n-1}=y_{n-1}\right)q_{sD}(t_{n-1}).$$

$$\ell (\mathbf{\theta })=\sum _{i=1}^N\log L_i=\sum _{i=1}^N\sum _{j=2}^{n_i}\log L_{ij},$$

$$\displaystyle\ell_{p}(\mbox{\boldmath$\theta$})$$

$$\mathbf{P}(t_1,t_n)=\mathbf{P}(t_1,t_2)\times \cdots \times \mathbf{P}(t_{n-1},t_n),$$

$$\mathbf{P}(t_1,t_2)=\mathbf{P}\left(t_1,u_{j_1+1}\right)\mathbf{P}\left(u_{j_1+1},u_{j_1+2}\right)\times \cdot \cdot \cdot \times \mathbf{P}\left(u_{j_2},t_2\right),$$

$$\mathbf{P}(t)=\mathbf{A}\text{mbox}diag\left(e^{b_{1}t},...,e^{b_{D}t}\right){\mathbf{A}}^{-1}.$$

$$\frac{\partial }{\partial {\theta }_k}\mathbf{P}(t)=\mathbf{A}{\mathbf{V}}_k{\mathbf{A}}^{-1},$$

$$\displaystyle\left\{\begin{array}[]{ll}g_{lm}^{(k)}\left[\exp(b_{l}t)-\exp(b_{m}t)\right]/(b_{l}-b_{m})&l\neq m\\ \\ g_{ll}^{(k)}t\exp(b_{l}t)&l=m,\end{array}\right.$$

$$\sum _{i=1}^N\sum _{j=2}^{n_i}\frac{\partial }{\partial {\theta }_k}\log L_{ij}.$$

$$\sum _{i=1}^N\sum _{j=2}^{n_i}\frac{\partial }{\partial {\theta }_k}\log L_{ij}\frac{\partial }{\partial {\theta }_l}\log L_{ij}.$$

$${\mathbf{S}}_p(\mathbf{\theta })$$

$$\text{vspace}4cm{\mathbf{M}}_p(\mathbf{\theta })$$

$$\displaystyle\bm{\theta}^{(v+1)}=\bm{\theta}^{(v)}+\mathbf{M}_{p}\big{(}\bm{\theta}^{(v)}\big{)}^{-1}\mathbf{S}_{p}\big{(}\bm{\theta}^{(v)}\big{)}\,.$$

$$\displaystyle\mbox{AIC}(\mbox{\boldmath$\lambda$})=-2\ell_{p}+2df.$$

$$\displaystyle df(\mbox{\boldmath$\lambda$})=\mbox{tr}[\mathbf{M}(\mathbf{M}+\mathbf{J}(\mbox{\boldmath$\lambda$}))^{-1}],$$

$$\displaystyle q_{rs}(t)=\exp\big{(}\beta_{rs.0}\big{)},$$

$$\displaystyle q_{rs}(t)=\exp\big{(}\beta_{rs.0}+\xi_{rs}t\big{)},$$

$$\displaystyle q_{rs}(t)=\exp\big{(}\beta_{rs.0}+\xi_{rs}t+\beta_{rs.1}sex+\beta_{rs.2}education\big{)},$$

$$\begin{array}{rl}{q}_{23}(t)& =\exp \left(\sum _{k=1}^K{\alpha }_{23.k}{B}_k+{\beta }_{23.1}sex+{\beta }_{23.2}education\right)\\ q_{34}(t)& =\exp \left(\sum _{k=1}^K{\alpha }_{34.k}{B}_k+{\beta }_{34.1}sex+{\beta }_{34.2}education\right).\end{array}$$

$$q_{34}(t)$$

$$\displaystyle\widehat{\mathbf{P}}\left(t_{1}=11,t_{2}=13\Big{|}\begin{tabular}[]{l}$sex=1,$\\ $education=1\!\!$\end{tabular}\!\!\right)={\small\left(\begin{array}[]{lllll}0.330&0.488&0.154&0.010&0.018\\ 0.171&0.531&0.253&0.023&0.022\\ 0.083&0.391&0.429&0.066&0.031\\ 0.034&0.219&0.410&0.291&0.046\\ 0&0&0&0&1\end{array}\right),}$$

$$\displaystyle\widehat{\mathbf{P}}\left(t_{1}=11,t_{2}=13\Big{|}\begin{tabular}[]{l}$sex=1,$\\ $education=1\!\!$\end{tabular}\!\!\right)={\small\left(\begin{array}[]{lllll}0.330&0.484&0.153&0.011&0.021\\ 0.170&0.529&0.253&0.024&0.024\\ 0.083&0.390&0.427&0.069&0.032\\ 0.034&0.219&0.407&0.292&0.048\\ 0&0&0&0&1\end{array}\right),}$$

Full text

††journal: ArXiv

1 Introduction

Multi-state models are routinely used in research where change of status over time is of interest. In epidemiology and medical statistics, the models are used to describe health-related processes over time, where status is defined by a disease or a condition. In social statistics and in demography, the models are used to study processes such as region of residence, work history, or marital status. A multi-state model which includes a dead state is called a multi-state survival model.

The specification of a multi-state survival model depends partly on the study design which generated the longitudinal data that are under investigation. An important distinction is whether or not exact times are observed for transitions between the states. This paper considers study designs where death times are known exactly (or right censored) and where transition times between the living states are interval censored. Many applications in epidemiology and medical statistics have this property as it is often hard to measure the exact time of onset of a disease or condition. Examples are dementia, cognitive decline, disability in old age, and infectious diseases.

A multi-state survival model describes change in a discrete longitudinal outcome variable and attrition due to death. If a longitudinal outcome variable can be adequately described by a set of states, then a multi-state survival model is an alternative to so-called joint models. An example of the latter is the shared-parameter joint model which consists of combining a survival model for the event time with a mixed-effects model for the longitudinal outcome [1].

This paper defines continuous-time multi-state survival models by specifying transition-specific hazard models. Time-dependency of the process is defined by using parametric and semi-parametric formulations in the specification of the baseline hazard functions. The semi-parametric specification is made with $P$-splines [2], which are $B$-splines with penalties on the difference of adjacents splines [2]. Using $B$-splines is a general method for smoothing [2, 3].

Because the definition of the functional form of the hazard can be transition specific, model specification can cover a wide range of multi-state survival processes. The methodology presented in this research is for multi-state models with at least one transition hazard specified with $P$-splines. Estimation is carried out with penalised maximum likelihood, where the maximisation is undertaken by using a Fisher scoring algorithm. This algorithm is an extension of the work by Jennrich and Bright [4] and Kalbfleisch and Lawless [5].

The models are formulated in a Markov process framework. Time-dependency is approached by using a piecewise-constant approximation and defining a series of time-homogeneous processes. For each of these homogeneous processes, the solution to the Kolmogorov forward equations (a first-order differential equation) is computed using eigenvalue decomposition. The method can be applied to multi-state models with any number of states and specification of hazard transitions can vary across transitions.

Semi-parametric multi-state models for interval-censored data have been discussed in the literature. Titman [6] uses a numerical approximation to calculate the transition probabilities. The advantage is that there is no need to define a grid for a piecewise-constant approximation, but it is computationally more demanding than using eigenvalue decomposition—especially in the case of continuous-scale covariates. Also, even though $B$-splines are used to model transitions intensities, the log-likelihood is maximised without penalisation. Joly and Commenges [7] use a penalised approach for a progressive three-state model. Estimation is performed with an algorithm which uses derivatives of the penalised log-likelihood. The smoothing parameters of the model are selected using a grid search with cross-validation. Joly et al. [8] use the same approach for an illness-death model. The method used in both papers requires explicit expressions for the transition probabilities. Calculating these formulas can be intractable for more complex models, such as models with more than four states and recovery [9].

Sennhenn-Reulen and Kneib [10] developed an estimation procedure for multi-state models based on a structured lasso penalisation. The aim of their research is to identify covariate effects coefficient equal to zero. Baseline transition intensities are specified with piecewise-constant models or unspecified and equal across all transitions. Their method is not defined for interval-censored data. Therefore, their work is different from ours in scope and methodology.

To illustrate the statistical modelling and the penalised maximum likelihood estimation, longitudinal data on survival and change of cognitive function in older population is analysed. The data stem from the English Longitudinal Study of Ageing (ELSA, www.elsa-project.ac.uk) and the longitudinal response variable is the number of words remembered in a recall from a list of ten. Of interest is the effect of age and gender on cognitive change over time when controlling for education. Four states are defined by the number of words an individual can remember, see Figure 1. The dead state is the fifth state.

The transition times between the living states are interval-censored, but death times are known. We acknowledge that these data are also analysed in Van den Hout [11]. This paper extends that analysis by using $P$-splines for hazard specification.

2 Models

For a continuous-time Markov chain $Y(t)$ on finite state space $\mathcal{S}$, time-homogeneous transition probabilities are given by

[TABLE]

for $r,s\in\mathcal{S}$ and $t\geq 0$. Transition matrix $\mathbf{P}(t)$ contains these probabilities such that the rows sum up to 1. The Chapman-Kolmogorov equation is $\mathbf{P}(u+t)=\mathbf{P}(u)\mathbf{P}(t)$. The transition intensities (or hazards) are given by

[TABLE]

for $r\neq s$. The matrix with off-diagonal entries $q_{rs}$ and diagonal entries $q_{rr}=-\sum_{r\neq s}q_{rs}$ is the generator matrix $\mathbf{Q}$. Given $\mathbf{Q}$, the solution for $\mathbf{P}(t)$ subject to $\mathbf{P}(0)=\mathbf{I}$ is $\mathbf{P}(t)=\exp(t\mathbf{Q})$, see, e.g., Norris [12]. In general, the computation of the exponential of a square matrix is not straightforward, see Moler and Van Loan [13] for a discussion of methods and efficiency.

A time-dependent hazard regression model for transition intensities combines baseline hazards with log-linear regression and is given by

[TABLE]

where $\mathbf{x}$ is a covariate vector without an intercept. Transition-specific time dependency can be introduced via baseline hazards. Parametric examples are

[TABLE]

Semi-parametric models can be defined in a similar way. An example is using $P$-splines to allow for flexible modelling of the time-dependency. Applications to multi-state models can be found in Kneib and Hennerfeind [14]. The semi-parametric formulation with $P$-splines of the baseline hazard is

[TABLE]

where—in this case—the choice of the number of knots $K$ is the same for all transition, but the $\alpha$s are not. Flexible multi-state models can be defined by $P$-splines or a combinations of the hazard specifications above.

3 Penalised maximum likelihood estimation

3.1 Likelihood function

Given a multi-state survival model, maximum likelihood inference can be used to analyse longitudinal data. In the presence of interval censoring, the likelihood function is constructed using transition probabilities. Let the state space be $\mathcal{S}=\{1,2,..,D\}$, with $D$ the dead state.

Consider a series of states $Y_{1},...,Y_{n}$ observed at times $t_{1},...,t_{n}$, respectively. The inference is conditional on the first observed state. For $Y_{2},...,Y_{n}$, the distribution is

[TABLE]

where ${\bm{\theta}}$ is the vector with the model parameters, $\mathbf{t}=(t_{1},...,t_{n})^{\top}$, and the $n\times p$ matrix $\mathbf{X}$ contains the values of the $p$ covariates at each of the $n$ time points. A conditional first-order Markov assumption is used to define the distribution (4) of $Y_{2},...,Y_{n}$ as

[TABLE]

where $\mathbf{x}_{j-1}$ is the $(j-1)^{th}$ row in $\mathbf{X}$.

Next consider an individual $i$ with observed values $y_{1},...,y_{n-1}\in\mathcal{S}\backslash D$, and a last observation $y_{n}$ which is either a value in $\mathcal{S}$ or a code for right-censoring. The likelihood contribution for this individual is $L_{i}=\prod_{j=2}^{n}L_{ij}$, where

[TABLE]

If a living state at $t_{n}$ is observed, then $C(y_{n}|y_{n-1})=P\left(Y_{n}=y_{n}|Y_{n-1}=y_{n-1}\right)$, where part of the conditioning is ignored in the notation. If the state is right censored at $t_{n}$, then $C(y_{n}|y_{n-1})=\sum_{s=1}^{D-1}P\left(Y_{n}=s|Y_{n-1}=y_{n-1}\right)$. If the state at $t_{n}$ is $D$, then known time of death is taken into account by defining

[TABLE]

Given $N$ individuals, the log-likelihood function is given by

[TABLE]

where $n_{i}$ is the number of observation times for individual $i$.

Above definition of the likelihood function can also be found in Jackson [9]. Including time-dependency as defined by the models in Section 2, does not affect the basic structure of the likelihood function. Similar expressions of the likelihood function can be found in Kalbfleisch and Lawless [5], Kay [15], and Gentleman et al. [16].

3.2 Penalised log-likelihood function

For the semi-parametric multi-state model, at least one baseline hazard function is specified with $P$-splines. If a transition is defined by $P$-splines, we specify a large set of equidistant knots. To control the smoothness of the estimated curve, a penalty based on finite differences of the coefficient of adjacent $P$-splines is imposed to the log-likelihood function. Without loss of generality, suppose there are $K$ knots for each smoothed hazard. Let $\beta$, $\alpha$ and $\xi$ represent the vector of parameters associated to the parametric, semi-parametric and covariates components of a multi-state model, respectively. Let \mbox{\boldmath\theta}^{\top}=(\mbox{\boldmath\beta}^{\top},\mbox{\boldmath\alpha}^{\top},\mbox{\boldmath\xi}^{\top}) be the full set of parameters and l(\mbox{\boldmath\theta}) be the log-likelihood of a semi-parametric multi-state model. The penalised log-likelihood function is

[TABLE]

where \mbox{\boldmath\alpha}_{j}=(\alpha_{j1},\ldots,\alpha_{jK})^{\top}, $\lambda$ is the vector of smoothing parameters, $\mathbf{D}$ is the matrix representation of the difference operator $\Delta$ of adjacent $P$-splines [2] and \mathbf{J}(\mbox{\boldmath\lambda}) is the penalty matrix. \mathbf{J}(\mbox{\boldmath\lambda}) is a block diagonal matrix with blocks $\lambda_{j}\mathbf{D}^{\top}\mathbf{D}$ for penalising $P$-splines parameters and zeros elsewhere [17].

3.3 Piecewise-constant hazards

Let $\mathbf{P}(t,t+\Delta)$ denote the transition matrix for any time interval $(t,t+\Delta]$. Time-dependency in hazard model (1), implies that $\mathbf{P}(t_{1},t_{1}+\Delta)\neq\mathbf{P}(t_{2},t_{2}+\Delta)$ for $t_{1}\neq t_{2}$.

Time-dependency of the hazard can be taken into account by using a piecewise-constant approximation. Given consecutive times $t_{1},t_{2},\dots t_{n}$, define the transition matrix for $(t_{1},t_{n}]$ by

[TABLE]

where the matrices at the right-hand side are derived using generator matrices $\mathbf{Q}(t_{1}),\mathbf{Q}(t_{2}),\dots,\mathbf{Q}(t_{n-1})$, respectively.

In longitudinal data for continuous-time models, follow-up times often vary across individuals. If that is the case, the individual-specific follow-up times can be used to define the piecewise-constant approximation for the individual likelihood contributions. This implies that a transition probability such $P\left(Y_{j}=y_{j}|Y_{j-1}=y_{j-1}\right)$ is derived by using $\mathbf{Q}(t_{j-1})$ to compute $\mathbf{P}(t_{j-1},t_{j})$.

Instead of letting the data determine the grid for the piecewise-constant approximation, it is also possible to impose a grid, which is the same for all individual likelihood contributions [18]. In this case, time intervals in the data are embedded in the grid. For example, say the grid is defined by $u_{1},...,u_{M}$. For an observed time interval $(t_{1},t_{2}]$, determine $j_{1}$ and $j_{2}$ such that $u_{j_{1}}<t_{1}\leq u_{j_{1}+1}$ and $u_{j_{2}}<t_{2}\leq u_{j_{2}+1}$. The transition matrix for $(t_{1},t_{2}]$ is then defined by

[TABLE]

using generator matrices $\mathbf{Q}(u_{j_{1}}),\mathbf{Q}(u_{j_{1}+1}),\dots,\mathbf{Q}(u_{j_{2}})$, respectively. For this approach covariate values are needed at all grid points $u_{1},...,u_{M}$. For a covariate with a stochastic time-dependency, these values may not be available in the data.

3.4 Scoring algorithm

Given a piecewise-constant approximation to the time-dependency in the hazard model (1), a scoring algorithm can be used to maximise the logarithm of the likelihood function (9). A scoring algorithm solves maximum likelihood equations numerically by iteratively estimating a root of the first-order derivative of the log-likelihood function. The first-order derivative of the log-likelihood function is called the score function.

As in Section 3.1, let ${\bm{\theta}}=(\theta_{1},...,\theta_{q})$ be the vector with model parameters. The crucial step is to derive $\partial\mathbf{P}(t_{1},t_{2})/\partial\theta_{k}$ for a given time interval $(t_{1},t_{2}]$. The important aspects of the scoring algorithm are

• (i)

Because of the piecewise-constant approximation, the basic formulas for the time-homogeneous case in Kalbfleisch and Lawless [5] apply to the constituent intervals with constant hazards in the likelihood function.

• (ii)

By using an eigenvalue decomposition of a generator matrix $\mathbf{Q}(t)$, only the derivatives $\partial\mathbf{Q}(t)/\partial\theta_{k}$ are needed [4].

• (iii)

For the Weibull and the Gompertz hazard models, and for a model with $P$-splines, derivatives $\partial\mathbf{Q}(t)/\partial\theta_{k}$ are straightforward to derive.

• (iv)

The likelihood contributions for exact death times and right-censoring are made up of transition probabilities and transition hazards and can be dealt with by using (i) – (iii).

To specify the scoring algorithm, the derivative of a transition matrix is presented first. Given piecewise-constant intensities, the likelihood contribution for an observed time interval $(t_{1},t_{2}]$ is defined using a constant generator matrix $\mathbf{Q}=\mathbf{Q}(t_{1})$. For the eigenvalues of $\mathbf{Q}$ given by $\mathbf{b}=(b_{1},...,b_{D})$, define $\mathbf{B}=\mbox{diag}(\mathbf{b})$. Given matrix $\mathbf{A}$ with the eigenvectors as columns, the eigenvalue decomposition is $\mathbf{Q}=\mathbf{A}\mathbf{B}\mathbf{A}^{-1}$. The transition probability matrix $\mathbf{P}(t)=\mathbf{P}(t_{1},t_{2})$ for elapsed time $t=t_{2}-t_{1}$ is given by

[TABLE]

As described in Kalbfleisch and Lawless [5], the derivative of $\mathbf{P}(t)$ can be obtained as

[TABLE]

where $\mathbf{V}_{k}$ is the $D\times D$ matrix with $(l,m)$ entry

[TABLE]

where $g_{lm}^{(k)}$ is the $(l,m)$ entry in $\mathbf{G}^{(k)}=\mathbf{A}\partial\mathbf{Q}/\partial\theta_{k}\mathbf{A}^{-1}$.

For the parametric and semi-parametric time-dependent hazard models in Section 2, matrix $\partial\mathbf{Q}(t_{1})/\partial\theta_{k}$ is straightforward to derive.

The scoring algorithm can now be defined as follows. Let the $q\times 1$ vector $\mathbf{S}(\bm{\theta})$ denote the score function. The $k$th entry of $\mathbf{S}(\bm{\theta})$ is given by

[TABLE]

The expected observed information matrix is called the Fisher information and is given by $\bm{\mathcal{I}}(\bm{\theta})=I\!\!E\hskip 0.85355pt\left[\mathbf{S}(\bm{\theta})\mathbf{S}(\bm{\theta})^{\top}\right]$, which can be estimated by defining the $q\times q$ matrix $\mathbf{M}({\bm{\theta}})$ with $(k,l)$ entry

[TABLE]

The penalised score $\mathbf{S}_{p}(\bm{\theta})$ and estimated penalised Fisher information matrix $\mathbf{M}_{p}(\bm{\theta})$ are given by

[TABLE]

Given starting values $\bm{\theta}^{(0)}$, the scoring algorithm is given for $v=1,2,3\dots$ by

[TABLE]

Let $\bm{\mathcal{I}}_{p}(\bm{\theta})$ represents the penalised Fisher information matrix. The asymptotic covariance matrix of the penalised maximum likelihood estimate $\widehat{\bm{\theta}}$ is equal to $\bm{\mathcal{I}}_{p}(\bm{\theta})^{-1}$. Hence, after convergence, the covariance matrix of the penalised maximum likelihood estimate $\widehat{\bm{\theta}}$ is estimated by $\mathbf{M}_{p}(\widehat{\bm{\theta}})^{-1}$ [19].

4 Estimation of smoothing parameter

Estimating optimal value for the smoothing parameters $\lambda$ is crucial for fitting models with splines [20]. A common method for choosing smoothing parameters is the Akaike Information Criterion (AIC). The AIC definition is equivalent to

[TABLE]

The degrees of freedom $df$ is a measure of model complexity. For parametric models, the degrees of freedom are equal to the number of independent parameters in the model. For semi-parametric models with $P$-splines, these can be defined as

[TABLE]

where $\mathbf{M}$ is the (estimated) Fisher information matrix and $\mathbf{J}$ is the penalty function [17]. Small values of $\lambda$ lead to wiggly functions, while large values lead to more conservative estimated functions that tend to a straight line.

5 Prediction

Once a multi-state model is fitted using a parametric and semi-parametric hazard model, estimated model parameters can be used for prediction. Typically this concerns computing transition matrices as a function of the penalised maximum likelihood estimate. The covariance of a function of model parameters can be estimated by Monte Carlo simulation or by using the multivariate delta method, see also Titman [6]. Because transition probabilities are restricted to $[0,1]$, using simulation is recommended as the default method. The delta method does not take the restriction into account and this can have a substantial knock-on effect on long-term prediction. This paper focuses on the simulation method.

Let $\widehat{\mathbf{V}}_{\bm{\theta}}$ denote the estimated covariance matrix of the penalised maximum likelihood estimate $\widehat{\bm{\theta}}$. Of interest is the estimation of $\mathbf{P}(t_{1},t_{2})$ for arbitrary $t_{1}$ and $t_{2}>t_{1}$. In the case of a time-dependent model, let the grid for the piecewise-constant approximation be defined by $u_{j+1}=u_{j}+h$ for $j=1,...,M$ such that $u_{1}=t_{1}$ and $u_{M}=t_{2}$. Given this grid, matrix $\mathbf{P}(t_{1},t_{2})$ is estimated by $\mathbf{P}(u_{1},u_{2})\times\cdot\cdot\cdot\times\mathbf{P}(u_{M-1},u_{M})$.

For Monte Carlo simulation, parameter vectors ${\bm{\theta}}^{(b)}$ are drawn from $N(\widehat{\bm{\theta}},\widehat{\mathbf{V}}_{\bm{\theta}})$, for $b=1,...,B$, and for each sampled ${\bm{\theta}}^{(b)}$, $\mathbf{P}(t_{1},t_{2})$ is calculated. Summary statistics such as mean and covariance can be derived easily from the $B$ realisations of $\mathbf{P}(t_{1},t_{2})$.

Sampling from a $K$-variate normal distribution $N(\bm{\mu},\bm{\Sigma})$ is possible by using the Cholesky decomposition $\bm{\Sigma}=\mathbf{L}\mathbf{L}^{\top}$. First $K$ draws are taken independently from the standard normal and collected in the $K\times 1$ vector $\mathbf{z}$. A multivariate draw from $N(\bm{\mu},\bm{\Sigma})$ is then given by $\bm{\mu}+\mathbf{L}\mathbf{z}$.

6 Application

6.1 English Longitudinal Study of Ageing (ELSA)

To illustrate the methodology, longitudinal data are analysed from the English Longitudinal Study of Ageing (ELSA). The ELSA baseline (1998-2001) is a representative sample of the English population aged 50 and older. Data from ELSA can be obtained via the Economic and Social Data Service (www.esds.ac.uk). There are 11932 individuals in the ELSA baseline.

Of interest for the current analysis is the change of cognitive function in older population. For the current analysis, a random sample of size $N=1000$ is taken from ELSA. Of these 1000 individuals, 205 died during the follow-up with age at death available. Because ELSA data are publicly available, measures have been taken by the data provider to prevent identification of the individuals. One of those measures is the censoring of ages above 90 years. In the sampling of the subset of $N=1000$, individuals who were 90 years or older at baseline were ignored. The sample has 544 women and 456 men.

Highest educational qualification is dichotomised for the current analysis according to years of formal education: fewer than ten versus ten or more. There are 558 individuals with fewer than ten years of education.

6.2 A five-state model for remembering words

This application focuses on the number of words remembered in a delayed recall from a list of ten. The score on this test is equal to the number of words remembered, i.e., score $\in\{0,1,2,\cdots,10\}$. The top graph in Figure 2 provides information on the number of words remembered at baseline. Most people remember 4 or 5 words, and the data show that remembering 9 or 10 words is exceptional. The bottom graph in Figure 2 depicts the change of number of words remembered over time for a random subset of 30 individuals. The 30 trajectories illustrate that the delayed recall is a noisy process. Nevertheless, already in these trajectories there is some evidence of a decline in cognitive function as people get older. The statistical modelling in this section aims to explore the effect of age and gender on cognitive change over time when controlling for education.

Four living states are defined by the number of words an individual can remember: state 1, 2, 3, and 4, for the number of words $\{7,8,9,10\}$, $\{6,5\}$ $\{4,3,2\}$, and $\{1,0\}$, respectively. An additional state 5 is defined as the dead state, see also Figure 1.

The interval-censored multi-state process is summarised by the frequencies in Table 1. Note that the sum of the transitions into the dead state is equal to the number of deaths in the sample, i.e., 205. Table 1 also shows that the process is mainly progressive in the sense that the main trend over time is towards the higher states.

In what follows, model estimation is undertaken by using the scoring algorithm. Let ${\bm{\theta}}^{\top}=(\theta_{1},...,\theta_{q})$ be the vector with model parameters, where $q$ depends on the chosen model. The convergence criterion for the algorithm is to stop at iteration $v+1$ when $\sum_{k=1}^{q}|\theta_{k}^{(v)}-\theta_{k}^{(v+1)}|<10^{-6}$.

Model selection is bottom-up starting with the time-homogeneous exponential hazard model given by

[TABLE]

for the transitions $r\rightarrow s$ depicted in Figure 1. This intercept-only model with 10 parameters has AIC = 8109.5. Convergence of the scoring algorithm was reached after 14 iterations, using starting values $\beta_{rs.0}=-3$ for all the parameters.

For the process at hand, age is the most suitable time scale. Age in the ELSA data is transformed by subtracting 49 years. This results in 1 being the minimal age in the sample.

Even though the sample size is not small, Table 1 shows that mortality information is limited because only about 20% of the individuals end up in the dead state during follow-up. In what follows, model (11) is extended by adding parameters with parameter equality constraints. The first extension is a Gompertz model given by

[TABLE]

where $\xi_{21}=\xi_{32}=\xi_{43}=\xi_{B}$ and $\xi_{15}=\xi_{25}=\xi_{35}=\xi_{45}=\xi_{D}$. That is, the effect of time is the same for all backwards transitions and for transitions into the dead state. In the estimation, the grid for the piecewise-constant approximation is defined by individually observed follow-up times in the data. This model has 15 parameters, and needs 16 scoring iterations when using starting values $\beta_{rs.0}=-3$ and $\xi_{rs}=0$ for all the relevant $r,s$-combinations. The model has AIC = 7780.5.

Subsequently, covariate information is added for the transitions of interest, i.e., those transitions that represent a decline in cognitive function. For this, model (12) is extended to

[TABLE]

where $sex$ is 0/1 for women/men, and $education$ is 0/1 for fewer than ten years/ten years of more of education. For the transitions into the dead state, the constraints on the coefficients for $sex$ are $\beta_{15.1}=\beta_{25.1}=\beta_{35.1}=\beta_{45.1}$, and for $education$ are $\beta_{r5.2}=0$ for $r=1,2,3,4$. This model has 22 parameters, needs 16 iterations, and has AIC = 7680.3.

It is worthwhile to investigate alternative time-dependent models. First, in model (13), the Gompertz baseline models for the transitions into the dead state are replaced by Weibull models. Starting values for the transitions into the dead state are $\beta_{r5.0}=-10$, $\tau_{15}=\exp(0.5)$, and for the remaining parameters the values are as given above. This yields AIC = 7688.7 after 20 iterations.

Next, all baseline hazards definitions in model (13) are replaced by Weibull models, which results in AIC = 7729.5 after 28 iterations. Alternatively, model (13) is defined with Gompertz baseline models for the transitions into the dead state and Weibull models for progression through the living states. This yields AIC = 7719.7 after 25 iterations.

Semi-parametric models with $P$-splines can be used to model non-linear functional forms and to check shapes specified by parametric models. To illustrate this, the Gompertz hazard for transition $2\rightarrow 3$ and $3\rightarrow 4$ in model (13) are replaced by

[TABLE]

For this model, the number of $P$-splines bases for both hazard functions is $K=10$ and the vector of smoothing parameters is \mbox{\boldmath\lambda}^{\top}=(\lambda_{1},\lambda_{2}). The initial grid is given by all pairs of combinations of $\log_{10}\lambda_{1}=\{-3,-2,-1,0,1,2,3\}$ and $\log_{10}\lambda_{2}=\{-3,-2,-1,0,1,2,3\}$. A possible graphical representation of the AIC results is to plot its values when one smoothing parameter is fixed. Figure 4(c) illustrates the resulting AIC for different values of $\lambda_{2}$ with fixed $\lambda_{1}=10^{-3}$. The value which minimises the AIC is $\lambda_{2}=10$. It happens for all values of $\lambda_{1}$. The search for the optimal values of $\lambda_{1}$ is less straightforward as $\lambda_{1}\rightarrow\infty$. Figure 4(a) shows the AIC for several values of $\lambda_{1}$ with fixed $\lambda_{2}=10$. The AIC decreases quickly for small values of $\lambda_{2}$; however, it gets constant for large values. This result indicates that the functional form of the hazard for transitions $2\rightarrow 3$ is $\log$-linear. Because both AIC and parameter estimates do not change much for sufficiently large values of $\lambda_{1}$, it is possible to set $\lambda_{1}=10^{7}$. In this case, the best model ($P$-splines I) according to the AIC is obtained with smoothing parameter \widehat{\mbox{\boldmath\lambda}}^{\top}=(10^{7},10). This model has 30 parameters and 26.1 degrees of freedom.

The fitted hazards for transition $2\rightarrow 3$ for the Gompertz (13) and $P$-splines I (14), for men with ten or more years of education are illustrated in Figure 4(b). The functional forms of both models are very similar for this transition; however, the functional forms for transition $3\rightarrow 4$ are quite different, as indicated in Figure 4(d). Model (14) has AIC = 7678.2 indicating that it performs better than the Gompertz model with AIC = 7680.3.

The functional form of hazard for transition $2\rightarrow 3$ in model (14) indicates that a Gompertz specification can be reasonable for this transition. Therefore, in model (13), only the hazard for transition $3\rightarrow 4$ is specified with $P$-splines:

[TABLE]

The number of $P$-splines bases is $K=10$ and the grid search is made on the values $\log_{10}\lambda=\{-3,-1,0,1,3\}$. The resulting AIC values are illustrated in Figure 4(e). The minimum AIC with value 7678.2 is obtained at $\lambda=10$. That is the same AIC value as for model (14); however, the degrees of freedom is slightly smaller $df=23.65$. As model (15) ($P$-splines II) is easier to estimate if compared to model (14), it is considered the best model among all illustrated in this paper. Table 2 summarises the comparison of the investigated models. Figure 4(f) illustrates the fitted hazard for transition $3\rightarrow 4$ in model (15) for men with ten or more years of education. As expected, there is an increase of risk of progression to a decline of cognitive function over the years.

Model validation is hampered by the interval censoring of the transitions between the living states. But given that death times are available, it make sense to compare survival as estimated by the model with Kaplan-Meier curves [16]. Of course, this will only check part of the fitted model. Figure 3 depicts baseline-specific survival as estimated by the model and as described by the Kaplan-Meier curves. For survival given baseline state 3, there is some discrepancy between model-based mean survival and the Kaplan-Meier curve, but overall the fit is reasonably good. Although this is not a proper goodness-of-fit test, the comparison shows that the model is able to capture the attrition due to death during the follow-up.

Table 3 illustrates the parametric component estimates for the $P$-splines II model. Most of the point estimates are according to expectation. For example, the effect of getting older is associated with decline of cognitive function, i.e., $\widehat{\xi}_{12}$, and with a decreasing hazard of remembering more words, i.e., $\widehat{\xi}_{B}<0$. For transitions $1\rightarrow 2$, $2\rightarrow 3$, and $3\rightarrow 4$ more years of education is associated with a lower risk of moving.

6.3 Predicting cognitive function

Although parameters for the transition intensities help to understand the estimated model, interpretation is more straightforward when transition probabilities are considered. Firstly, consider a short time interval for which we assume that the intensities are constant. For men aged 60 with ten or more years of education, the two-year transition probabilities are estimated at

[TABLE]

where $t$ denotes age transformed by subtracting 49 years. The diagonal entries in this matrix dominate. But there are some large off-diagonal entries as well. For example, if a man aged 60 is in state 3, then he has a 39% chance of being in state 2 two years later. This high chance is an illustration of the noisiness of the process under investigation: it is quite likely that a 60 year old man moves between states 2 and 3 within the next two years.

Next we illustrate the estimation of standard errors and 95% confidence intervals for transition probabilities. Using simulation with $B=1000$, we obtain

[TABLE]

with estimated standard errors

[TABLE]

The 95% confidence intervals for the first row are given by

[TABLE]

Next, ten-year transition probabilities are estimated for men aged 60 with ten or more years of education. The grid is defined by $h=1/2$ years. The estimation is shown in Figure 5.

Figure 5 concurs with the expectations. For example, given the progressive trend of the process, it is to be expected that probability of being in state 3 decreases over time, as moving to states 4, and 5 becomes more likely due to increased age.

7 Discussion

Specification and estimation of continuous-time multi-state survival models are presented and shown to be a flexible framework for statistical modelling of time-dependency processes. By defining transition-specific parametric and semi-parametric hazard models, a wide range of multi-state processes can be investigated. Penalised maximum likelihood estimation is undertaken by a scoring algorithm using a piecewise-constant approximation to time-dependent hazards. The Akaike information criterion is used to select the optimal value for the smoothing parameters.

The Markov process formulation to semi-parametric multi-state models extends the method described in Joly and Commenges [7] and Joly et al. [8]. This is an important methodology to medical studies as backwards transitions occur naturally in many applications [21, 22]. Furthermore, using the piecewise-constant approximation is an alternative to the method introduced by Titman [6] which handles the time-dependency by using numerical solutions to the non-linear differential equations which are defined directly by the time-dependency of the Markov process. As stated by Titman, computation using the non-linear differential equations can become prohibitively slow when adding continuous covariates. This is not a problem when using the piecewise-constant approximation and the scoring algorithm.

The scoring algorithm is implemented in R in such a way that it is easy to vary transition-specific choices for parametric and semi-parametric shapes. An example of such a model is explored in the application, where $P$-splines are used for transitions $2\rightarrow 3$ and $3\rightarrow 4$, and Gompertz hazards are defined for the other transitions. The eigenvalue decomposition in the algorithm is computed with the function eigen in R, which uses the LAPACK routine [23]. $P$-spline bases are computed using the code in the appendix in Eilers and Marx [2].

If prediction of a time-dependent process beyond the time range in the data is of interest, hazard models with $P$-splines can be used to validate the parametric choices which underlie the prediction. This was illustrated in the application with the ELSA data in which age range is from 50 to 90 years. If risk factors are the main focus of the research, $P$-splines can be used to capture non-parametric shapes of time-dependency.

The choice of the type of spline is not essential. $P$-splines were used in this paper, but any other spline function with a first-order derivative can be handled within the current framework. The same holds for parametric shapes other than the Gompertz and the Weibull. The specification and estimation of continuous-time survival model is very general and does not pose restrictions on the number of states, scale of covariates, or number of transitions.

Acknowledgements

The ELSA data were made available through the UK Economic and Social Data Service. ELSA was developed by a team of researchers based at the National Centre for Social Research, University College London and the Institute for Fiscal Studies. The developers and funders of ELSA do not bear any responsibility for the analyses or interpretations presented here. This research was supported by CNPq - Brazil [249308/2013-4].