Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

Chapter Contents
Chapter Introduction
NAG Toolbox

# NAG Toolbox Chapter IntroductionG12 — Survival Analysis

## Scope of the Chapter

This chapter is concerned with statistical techniques used in the analysis of survival/reliability/failure time data.
Other chapters contain functions which are also used to analyse this type of data. Chapter G02 contains generalized linear models, Chapter G07 contains functions to fit distribution models, and Chapter G08 contains rank based methods.

## Background to the Problems

### Introduction to Terminology

This chapter is concerned with the analysis on the time, t$t$, to a single event. This type of analysis occurs commonly in two areas. In medical research it is known as survival analysis and is often the time from the start of treatment to the occurrence of a particular condition or of death. In engineering it is concerned with reliability and the analysis of failure times, that is how long a component can be used until it fails. In this chapter the time t$t$ will be referred to as the failure time.
Let the probability density function of the failure time be f(t)$f\left(t\right)$, then the survivor function, S(t)$S\left(t\right)$, which is the probability of surviving to at least time t$t$, is given by
 ∞ S(t) = ∫ f(τ)dτ = 1 − F(t) t
$S(t)=∫t∞f(τ)dτ=1-F(t)$
where F(t)$F\left(t\right)$ is the cumulative density function. The hazard function, λ(t)$\lambda \left(t\right)$, is the probability that failure occurs at time t$t$ given that the individual survived up to time t$t$, and is given by
 λ(t) = f(t) / S(t). $λ(t)=f(t)/S(t).$
The cumulative hazard rate is defined as
 t Λ(t) = ∫ λ(τ)dτ, 0
$Λ (t)=∫0tλ (τ)dτ ,$
hence S(t) = eΛ(t)$S\left(t\right)={e}^{-\Lambda \left(t\right)}$.
It is common in survival analysis for some of the data to be right-censored. That is, the exact failure time is not known, only that failure occurred after a known time. This may be due to the experiment being terminated before all the individuals have failed, or an individual being removed from the experiment for a reason not connected with effects being tested in the experiment. The presence of censored data leads to complications in the analysis.

### Rank Statistics

There are a number of different rank statistics described in the literature, the most common being the logrank statistic. All of these statistics are designed to test the null hypothesis
• H0 : S1 (t) = S2 (t) = = Sg (t) , t τ ${H}_{0}:{S}_{1}\left(t\right)={S}_{2}\left(t\right)=\cdots ={S}_{g}\left(t\right),\forall t\le \tau$
where Sj${S}_{j}$ is the survivor function for group j$j$, g$g$ is the number of groups being tested and τ$\tau$ is the largest observed time, against the alternative hypothesis
• H1 : ${H}_{1}:$ at least one of the Sj (t) ${S}_{j}\left(t\right)$ differ, for some t τ $t\le \tau$.
A rank statistics T$T$ is calculated as follows:
Let ti ${t}_{i}$, for i = 1 , 2 , , nd $i=1,2,\dots ,{n}_{d}$, denote the list of distinct failure times across all g$g$ groups and wi${w}_{i}$ a series of nd${n}_{d}$ weights.
Let dij${d}_{ij}$ denote the number of failures at time ti${t}_{i}$ in group j$j$ and nij${n}_{ij}$ denote the number of observations in the group j$j$ that are known to have not failed prior to time ti${t}_{i}$, i.e., the size of the risk set for group j$j$ at time ti${t}_{i}$. If a censored observation occurs at time ti${t}_{i}$ then that observation is treated as if the censoring had occurred slightly after ti${t}_{i}$ and therefore the observation is counted as being part of the risk set at time ti${t}_{i}$.
Finally let
 g g di = ∑ dij   and  ni = ∑ nij . j = 1 j = 1
$di = ∑ j=1 g d ij and ni = ∑ j=1 g n ij .$
The (weighted) number of observed failures in the j$j$th group, Oj${O}_{j}$, is therefore given by
 nd Oj = ∑ wi dij i = 1
$Oj = ∑ i=1 nd wi d ij$
and the (weighted) number of expected failures in the j$j$th group, Ej${E}_{j}$, by
 nd Ej = ∑ wi ( nij di )/(ni) i = 1
$Ej = ∑ i=1 nd wi n ij di ni$
and if x$x$ denote the vector of differences x = (O1E1,O2E2,,OgEg) $x=\left({O}_{1}-{E}_{1},{O}_{2}-{E}_{2},\dots ,{O}_{g}-{E}_{g}\right)$
 nd Vjk = ∑ wi2 (( di (ni − di) (nin i k Ijk − nijnik) )/( ni2 (ni − 1) )) i = 1
$V jk = ∑ i=1 nd w i 2 ( di ( ni - di ) ( ni n i k I jk - n ij n ik ) n i 2 ( ni - 1 ) )$
where Ijk = 1 ${I}_{jk}=1$ if j = k$j=k$ and 0$0$ otherwise, then the rank statistic, T$T$, is calculated as
 T = x V− xT $T = x V- xT$
where V${V}^{-}$ denotes a generalized inverse of the matrix V$V$.
Under the null hypothesis, T χν2 $T\sim {\chi }_{\nu }^{2}$ where the degrees of freedom, ν$\nu$, is taken as the rank of the matrix V$V$.
The different rank statistics are defined by using different weights in the above calculations, for example
 logrank statistic wi = 1${w}_{i}=1$ Wilcoxon rank statistic wi = ni${w}_{i}={n}_{i}$ Tarone–Ware rank statistic wi = sqrt(ni)${w}_{i}=\sqrt{{n}_{i}}$ Peto–Peto rank statistic wi = S̃ (ti) ${w}_{i}=\stackrel{~}{S}\left({t}_{i}\right)$ where S̃ (ti) = ∏ tj ≤ ti   ( nj − dj + 1 )/(nj + 1) $\stackrel{~}{S}\left({t}_{i}\right)=\prod _{{t}_{j}\le {t}_{i}}\phantom{\rule{0.25em}{0ex}}\frac{{n}_{j}-{d}_{j}+1}{{n}_{j}+1}$

### Estimating the Survivor Function and Hazard Plotting

The most common estimate of the survivor function for censored data is the Kaplan–Meier or product-limit estimate,
 i Ŝ(t) = ∏ ((nj − dj)/(nj)),  ti ≤ t < ti + 1 j = 1
$S^(t)=∏j=1i (nj-djnj) , ti≤t
where dj${d}_{j}$ is the number of failures occurring at time tj${t}_{j}$ out of nj${n}_{j}$ surviving to tj${t}_{j}$. This is a step function with steps at each failure time but not at censored times.
As S(t) = eΛ(t)$S\left(t\right)={e}^{-\Lambda \left(t\right)}$ the cumulative hazard rate can be estimated by
 Λ̂(t) = − log(Ŝ(t)). $Λ^(t)=-log(S^(t)).$
A plot of Λ̂(t)$\stackrel{^}{\Lambda }\left(t\right)$ or log(Λ̂(t))$\mathrm{log}\left(\stackrel{^}{\Lambda }\left(t\right)\right)$ against t$t$ or logt$\mathrm{log}t$ is often useful in identifying a suitable parametric model for the survivor times. The following relationships can be used in the identification.
 (a) Exponential distribution: Λ(t) = λt$\Lambda \left(t\right)=\lambda t$. (b) Weibull distribution: log(Λ(t)) = logλ + γlogt$\mathrm{log}\left(\Lambda \left(t\right)\right)=\mathrm{log}\lambda +\gamma \mathrm{log}t$. (c) Gompertz distribution: log(Λ(t)) = logλ + γt$\mathrm{log}\left(\Lambda \left(t\right)\right)=\mathrm{log}\lambda +\gamma t$. (d) Extreme value (smallest) distribution: log(Λ(t)) = λ(t − γ)$\mathrm{log}\left(\Lambda \left(t\right)\right)=\lambda \left(t-\gamma \right)$.

### Proportional Hazard Models

Often in the analysis of survival data the relationship between the hazard function and the number of explanatory variables or covariates is modelled. The covariates may be, for example, group or treatment indicators or measures of the state of the individual at the start of the observational period. There are two types of covariate time independent covariates such as those described above which do not change value during the observational period and time dependent covariates. The latter can be classified as either external covariates, in which case they are not directly involved with the failure mechanism, or as internal covariates which are time dependent measurements taken on the individual.
The most common function relating the covariates to the hazard function is the proportional hazard function
 λ(t,z) = λ0(t)exp(βTz) $λ(t,z)=λ0(t)exp(βTz)$
where λ0(t)${\lambda }_{0}\left(t\right)$ is a baseline hazard function, z$z$ is a vector of covariates and β$\beta$ is a vector of unknown parameters. The assumption is that the covariates have a multiplicative effect on the hazard.
The form of λ0(t)${\lambda }_{0}\left(t\right)$ can be one of the distributions considered above or a nonparametric function. In the case of the exponential, Weibull and extreme value distributions the proportional hazard model can be fitted to censored data using the method described by Aitkin and Clayton (1980) which uses a generalized linear model with Poisson errors. Other possible models are the gamma distribution and the log-normal distribution.

### Cox's Proportional Hazard Model

Rather than using a specified form for the hazard function, Cox (1972) considered the case when λ0(t)${\lambda }_{0}\left(t\right)$ was an unspecified function of time. To fit such a model assuming fixed covariates a marginal likelihood is used. For each of the times at which a failure occurred, ti${t}_{i}$, the set of those who were still in the study is considered this includes any that were censored at ti${t}_{i}$. This set is known as the risk set for time ti${t}_{i}$ and denoted by R(ti)$R\left({t}_{i}\right)$. Given the risk set the probability that out of all possible sets of di${d}_{i}$ subjects that could have failed the actual observed di${d}_{i}$ cases failed can be written as
 (exp(siTβ))/( ∑ exp(zlTβ)) $exp(siTβ) ∑exp(zlTβ)$ (1)
where si${s}_{i}$ is the sum of the covariates of the di${d}_{i}$ individuals observed to fail at ti${t}_{i}$ and the summation is over all distinct sets of ni${n}_{i}$ individuals drawn from R(ti)$R\left({t}_{i}\right)$. This leads to a complex likelihood. If there are no ties in failure times the likelihood reduces to
 nd L = ∏ (exp(ziTβ))/([∑l ∈ R(ti)exp(zlTβ)]) i = 1
$L=∏i=1ndexp(ziTβ) [∑l∈R(ti)exp(zlTβ)]$
(2)
where nd${n}_{d}$ is the number of distinct failure times. For cases where there are ties the following approximation, due to Peto [2], can be used:
 nd L = ∏ (exp(siTβ))/([∑l ∈ R(ti)exp(zlTβ)]di). i = 1
$L=∏i=1ndexp(siTβ) [∑l∈R(ti)exp(zlTβ)]di .$
(3)
Having fitted the model an estimate of the baseline survivor function (derived from λ0(t)${\lambda }_{0}\left(t\right)$ and the residuals) can be computed to examine the suitability of the model, in particular the proportional hazard assumption.

## Recommendations on Choice and Use of Available Functions

The following functions are available.
 nag_surviv_kaplanmeier (g12aa) computes Kaplan–Meier estimates of the survivor function and their standard deviations. nag_surviv_logrank (g12ab) comparison of survival curves using rank statistics. nag_surviv_coxmodel (g12ba) fits the Cox proportional hazards model for fixed covariates. nag_surviv_coxmodel_risksets (g12za) creates the risk sets associated with the Cox proportional hazards model for fixed covariates.
Depending on the rank statistic required, it may be necessary to call nag_surviv_logrank (g12ab) twice, once to calculate the number of failures (di${d}_{i}$) and the total number of observations (ni${n}_{i}$) at time ti${t}_{i}$, to facilitate in the computation of the required weights, and once to calculate the required rank statistics.
The following functions from other chapters may also be useful in the analysis of survival data.
 nag_stat_mills_ratio (g01mb) the reciprocal of Mills' Ratio, that is the hazard rate for the Normal distribution. nag_correg_glm_poisson (g02gc) fits generalized linear model with Poisson errors (see Aitkin and Clayton (1980)). nag_correg_glm_gamma (g02gd) fits generalized linear model with gamma errors. nag_univar_estim_normal (g07bb) fits Normal distribution to censored data. nag_univar_estim_weibull (g07be) fits Weibull distribution to censored data. nag_nonpar_rank_regsn_censored (g08rb) fits linear model using likelihood based on ranks to censored data (see Kalbfleisch and Prentice (1980)). nag_contab_condl_logistic (g11ca) fits a conditional logistic model. When applied to the risk sets generated by nag_surviv_coxmodel_risksets (g12za) the Cox proportional hazards model is fitted by exact marginal likelihood in the presence of tied observations.

## Functionality Index

 Cox's proportional hazard model,
 create the risk sets nag_surviv_coxmodel_risksets (g12za)
 parameter estimates and other statistics nag_surviv_coxmodel (g12ba)
 Survival,
 Rank statistics nag_surviv_logrank (g12ab)
 Survivor function nag_surviv_kaplanmeier (g12aa)

## References

Aitkin M and Clayton D (1980) The fitting of exponential, Weibull and extreme value distributions to complex censored survival data using GLIM Appl. Statist. 29 156–163
Cox D R (1972) Regression models in life tables (with discussion) J. Roy. Statist. Soc. Ser. B 34 187–220
Gross A J and Clark V A (1975) Survival Distributions: Reliability Applications in the Biomedical Sciences Wiley
Kalbfleisch J D and Prentice R L (1980) The Statistical Analysis of Failure Time Data Wiley