Chapter 9 Multivariate Survival Models
Sometimes survival data come in clusters, and multivariate, or frailty, models are appropriate to use.
Ever since the paper by Vaupel, Manton, and Stallard (1979), the concept of frailty has spread in even wider circles of the research community. Although their primary purpose was to show various consequences of admitting individual frailties (“individuals age faster than cohorts”, due to the selection effect), the effect was that people started to implement their frailty model in Cox regression models.
9.1 An Introductory Example
Let us assume that in a follow-up study, the cohort is not homogeneous but instead consists of two equally sized groups with differing hazard rates. Assume further that we have no indication of which group an individual belongs to, and that members of both groups follow an exponential life length distribution: \[\begin{equation*} \begin{split} h_1(t) &= \lambda_1 \\ h_2(t) &= \lambda_2 \\ \end{split} \qquad t > 0. \end{equation*}\] This implies that the corresponding survival functions \(S_1\) and \(S_2\) are \[\begin{equation*} \begin{split} S_1(t) &= e^{-\lambda_1 t} \\ S_2(t) &= e^{-\lambda_2 t} \\ \end{split} \qquad t > 0, \end{equation*}\] and a randomly chosen individual will follow the “population mortality” \(S\), which is a mixture of the two distributions: \[\begin{equation*} S(t) = \frac{1}{2} S_1(t) + \frac{1}{2} S_2(t), \quad t > 0. \end{equation*}\] Let us calculate the hazard function for this mixture. We start by finding the density function \(f\): \[\begin{equation*} f(t) = -\frac{dS(x)}{dx} = \frac{1}{2}\left(\lambda_1 e^{-\lambda_1 t} + \lambda_2 e^{-\lambda_2 t} \right), \quad t > 0. \end{equation*}\] Then, by the definition of \(h\) we get \[\begin{equation} h(t) = \frac{f(t)}{S(t)} = \omega(t) \lambda_1 + \big(1 - \omega(t)\big) \lambda_2, \quad t > 0, \tag{9.1} \end{equation}\] with \[\begin{equation*} \omega(t) = \frac{e^{-\lambda_1 t}}{e^{-\lambda_1 t} + e^{-\lambda_2 t}} \end{equation*}\] It is easy to see that \[\begin{equation*} \omega(t) \rightarrow \left\{ \begin{array}{ll} 0, & \lambda_1 > \lambda_2 \\ \frac{1}{2}, & \lambda_1 = \lambda_2 \\ 1, & \lambda_1 < \lambda_2 \end{array} \right. , \quad \mbox{as } t \rightarrow \infty, \end{equation*}\] implying that \[\begin{equation*} h(t) \rightarrow \min(\lambda_1, \lambda_2), \quad t \rightarrow \infty, \end{equation*}\] see Figure 9.1.
FIGURE 9.1: Population hazard function (solid line). The dashed lines are the hazard functions of each group, constant at 1 and 2.
The important point here is that it is impossible to tell from data alone whether the population is homogeneous, with all individuals following the same hazard function given by equation (9.1), or if it in fact consists of two groups, each following a constant hazard rate. Therefore, individual frailty models like \(h_i(t) = Z_i h(t), \quad i = 1, \ldots, n\), where \(Z_i\) is the “frailty” for individual No. \(i\), and \(Z_1, \ldots, Z_n\) are independent and identically distributed (iid) are less useful.
A heuristic explanation to all this is the dynamics of the problem: We follow a population (cohort) over time, and the composition of it changes over time. The weaker individuals die first, and the proportion stronger will steadily grow as time goes by.
Another terminology is to distinguish between individual and population hazards. In Figure 9.1 the solid line is the population hazard, and the dashed lines represent the two kinds of individual hazards present. Of course, in a truly homogeneous population, these two concepts coincide.
9.2 Frailty Models
Frailty models in survival analysis correspond to hierarchical models in linear or generalized linear models. They are also called mixed effects models. A general theory, with emphasis on using R, of mixed effects models can be found in Pinheiro and Bates (2000).
9.2.1 The simple frailty model
Vaupel, Manton, and Stallard (1979) described an individual frailty model, \[\begin{equation*} h(t; \mathbf{x}, Z) = h_0(t) Z e^{\boldsymbol{\beta} \mathbf{x}}, \quad t > 0, \end{equation*}\] where \(Z\) is assumed to be drawn independently for each individual. Hazard rates for ``random survivors’’ are not proportional, but converging (to each other) if the frailty distribution has finite variance. Thus, the problem may be less pronounced in AFT than in PH regression. However, as indicated in the introductory example, with individual frailty the identification problems are large, and such models are best avoided.
9.2.3 Parametric frailty models
It is possible utilize the connection between Poisson regression and the piecewise constant proportional hazards model discussed in Chapters 7 and 8 to fit parametric frailty models. We look at the fertility data again. The standard analysis without frailty effects is shown in Table 9.4.
fit0 <- pchreg(Surv(next.ivl, event) ~ parity + ses,
cuts = 0:13, data = fe)| Covariate | level | W_mean | Coef | HR | SE | LR_p |
|---|---|---|---|---|---|---|
| parity | 4.242 | -0.1306 | 0.8775 | 0.0048 | 0.0000 | |
| ses | 0.0000 | |||||
| farmer | 0.491 | 0 | 1 | |||
| unknown | 0.186 | 0.0949 | 1.0995 | 0.0291 | ||
| upper | 0.018 | 0.0866 | 1.0904 | 0.0829 | ||
| lower | 0.305 | -0.1508 | 0.8600 | 0.0255 | ||
| Max Log | Likelihood | -14068.4 | ||||
| Restricted | mean | 2.299 |
Note the use of the function survival::survSplit (comment # 1 in the code),
the transformation to
factor for the slices of time (years) created by
survSplit, and the creation of an offset (# 3). This is described in
detail in Chapter 8.
For testing the presence of a random effect over women, an
alternative is to use the function glmmML in the package
with the same name.
library(glmmML)
fe13 <- survSplit(fe, end = "next.ivl", event = "event",
cut = 1:13, episode = "years",
start = "start") # 1
fe13$years <- as.factor(fe13$years) # 2
fe13$offs <- log(fe13$next.ivl - fe13$start) # 3
fit1 <- glmmML(event ~ parity + ses + years + offset(offs),
cluster = id, family = poisson, method = "ghq",
data = fe13, n.points = 9)
out <- with(fit1, cbind(coefficients, coef.sd))
colnames(out) <- c("Coef", "se(Coef)")
round(out[1:5, ], 3) Coef se(Coef)
(Intercept) -3.729 0.090
parity -0.274 0.006
sesunknown 0.086 0.065
sesupper -0.033 0.171
seslower -0.262 0.054The estimated scale parameter in the mixing distribution is 0.851 with standard error 0.024,
so the conclusion is very much the same as with the nonparametric
(coxme) approach: The clustering effect of mother is very
strong and must be taken into account in the analysis of birth intervals.
The nonparametric approach is easier to use and recommended, but see Section 9.3.
9.3 Stratification
A simple way to eliminate the effect of clustering is to stratify on the clusters. In the birth intervals example, it would mean that intervals are only compared to other birth intervals from the same mother. The drawback with a stratified analysis is that it is not possible to estimate the effect of covariates that are constant within clusters. In the birth intervals case, it is probable thatses, socio-economic status, would
vary little within families. On the other hand, the effect of
birth order or mother’s age would be suitable to analyze in a stratified
setting, see Table 9.5.
| Covariate | level | W_mean | Coef | HR | SE | LR_p |
|---|---|---|---|---|---|---|
| parity | 4.242 | -0.3287 | 0.7199 | 0.0078 | 0.0000 | |
| prev.ivl | 2.423 | -0.0536 | 0.9478 | 0.0157 | 0.0006 | |
| Max Log | Likelihood | -9557.9 |
| Covariate | level | W_mean | Coef | HR | SE | LR_p |
|---|---|---|---|---|---|---|
| parity | 4.242 | -0.0838 | 0.9196 | 0.0048 | 0.0000 | |
| prev.ivl | 2.423 | -0.3192 | 0.7267 | 0.0108 | 0.0000 | |
| Max Log | Likelihood | -70390.9 |
Note how the effect of parity is diminished when aggregating
comparison over all women, while the effect of the length of the previous
interval is enlarged. Try to explain why this result is expected!
Under certain circumstances it is actually possible to estimate an effect of a covariate that is constant within strata, but only if it is interacted with a covariate that is not constant within strata. See the example about maternal and infant mortality in Chapter 10.