\section*{Survival Analysis}

References:  Rupert G.~Miller~Jr., ``Survival Analysis'', Wiley 1981.

John Kalbfleisch and R.L.~Prentice, ``The Statistical Analysis of
Failure Time Data'', Wiley, 1980

Miller, Efron, et. al, ``Biostatistics Casebook'', Wiley, 1980.

Venables and Ripley, Chapter 11.

Cryer, J. D., ``Time Series Analysis'', Duxbury Press, 1986.\footnote{This is the most readable of these texts.}


\subsection*{Aim}

For both of us ({\bf you} and I, to learn something about the problems
and methods of survival analysis.  We will eventually take some data
sets and attempt to do at least some rudimentary, and perhaps somewhat
complicated, analyses on them.


\subsection*{Problem}

We observer a {\bf positive valued} random variable (the survivor time
or failure time) and some co-variates.  Typically the value of the
random variable is the time to failure of some physical component or
the time to death of a biological unit (or the time to relapse of some
disease, etc.).  There are all sorts of examples of waiting times
which can be similarly analyzed---e.g., the length of time required to
learn a new skill.

Early work in this area was some of the first (modern) statistics ever
recorded.  Graunt and Edmund Halley were responsible for the first
useful {\bf life tables} (estimate of the force of mortality) in about
166X.  Medical and insurance applications carried the work steadily
forward.  

The modern era was ushered in by the need to test mechanical
components for reliability.  Rapid industrial expansion in the early
1940's was a major impetus.  

The most novel feature of failure-time data is that notion of
censoring---not all sample units (e.g., patients in a medical trial)
will be observed to failure (e.g., death), yet even the ``partially''
observed data provide useful information on the lifetime
distribution.


\paragraph{Example}

Carcinogenesis.  Two groups of rats were given pre-treatments then
exposed to a carcinogen.  We would like to compare the pre-treatments.
Presumably, we could transform the data and perform a two-sample
t-test, except there is some censoring in the data.  Some of the rats
died of causes unrelated to the cancer, or survived the study.  
\begin{center}
Days to Vaginal Cancer Mortality in Rats
  \begin{tabular}{r l l l l l l l l l l}
\hline
Group 1 & 143& 164& 188& 188& 190& 192& 206& 209& 213& 216\\
        & 220& 227& 230& 234& 246& 265& 304& 216*& 244*\\
Group 2 & 142& 156& 163& 198& 205& 232& 232& 233& 233& 233\\
        & 233& 239& 240& 261& 280& 280& 296& 296& 323& 204*\\
        & 344* \\    \hline
  \end{tabular}
\end{center}
The items with (*) are those that are censored.  This is like a form
of missing data except the ``missing-ness'' of the data contains
information about the questions of interest.  


\paragraph{Example}

Accelerated Life Test.  Many physical components are {\bf designed} to
be very reliable and have long lifetimes. (and then there are the
other components).  It is still of interest to know something about
the expected lifetime of a component or a comparison between
differing designs.  When the lifetime is very long, it is often
difficult to actually test the component.  (Think about a computer
chip.  If you tested the lifetime in the standard way, the chip would
be obsolete before you saw any failures).

In these circumstances and ``accelerated life test'' is often
performed.  This consists of an experiment where the components are
subjected to {\em more\/} stress than would be expected under normal
operating conditions with the aim of extrapolating or comparing.

For a certain type of motor designed to operate ate 130 degrees C, the
following severely censored data were obtained.
\begin{description}
\item[150C] All 10 motorettes without failure at 8064 hours
\item[170C] 1764, 2772, 3444, 3542, 3780, 4860, 5196, 3 motorettes
  without failure at 5448 hours
\item[190C] 408, 408, 1344, 1344, 1440, 5 motorettes without failure
  at 1680 hours
\item[220C] 408, 408, 504, 504, 504, 5 motorettes without failure at
  528 hours
\end{description}

How would you analyze these data?  Suppose you wanted to compare two
types of motor?


\paragraph{African Children}

A Stanford psychiatrist wanted to know the age at which a group of
African children learned to perform a particular task---so he wrote a
grant proposal and went to Africa.  When he arrived in the village
some of the children could already perform the task, some children
learnt while he was there and some hadn't learnt when he left.  In
this case there is both left and right censoring.  For some
individuals we know an age by which they could perform the task (left
censoring); for some we know an age at which they could not perform
the task, and for the rest we know the age at which they learnt the
task.


\subsection*{Survival Functions and Hazard Rates}

Let the random variable $T > 0$ denote survival time.  Let $T$ have
density $f(t)$ and distribution
$$
F(t) = \mbox{Prob}(T < t) $$
(Note that some books use $F(t) = \mbox{Prob}(T > t)$ which makes
sense but is not standard).

The survival function is
$$
S(t) = 1 - F(t) = P(T > t) $$
and the hazard rate (function) is
$$
\lambda(t) = \frac{f(t)}{1 - F(t)}$$
which can be interpreted as
$$ 
\lambda(t) = \frac{\mbox{Prob}( t < T < d + \Delta t | T > t)}{\Delta
  t}
$$
which can be read as
$$
\mbox{Prob}\{\mbox{dying in} (t, t+\Delta t) | \mbox{lived to time }
t\} $$
or the instantaneous force of mortality.

Now, notice that if we integrate $\lambda(t)$ we get
\begin{eqnarray*}
  \int_0^t{\lambda(u) du}& = &\int_0^t{ \frac{f(u)}{1- F(u)} du} \\
  & = & - \log{[ 1- F(u)]_0^t}  \\
  & = & - \log{(S(t))}
\end{eqnarray*}

Thus
$$
S(t) = \exp{ \{ - \int_0^t{ \lambda(u) du}\} } $$

These various forms, definitions and {\bf ideas} are very important.
If you don't understand them now, stop me.  


\subsection*{Types of Censoring}

Let $T_1, \ldots,  T_n$ be iid $F$.
\begin{itemize}
\item Type I Censoring.  Let $t_c$ be some preassigned fixed number,
  the fixed censoring time.  Instead of seeing $T_1, \ldots, T_n$ we
  get to see
  $$ Y_i = \min{(T_i, t_c)} $$ so the density of $Y$ might look like \qquad 
   \psfig{figure=surv.eps,angle=-90,width=1.5in}
 \item Type II Censoring.    Let $r < n$ be fixed.  All we observe are
   $T(1), \ldots, T(r)$.
\end{itemize}

Both types of censoring are common---either we have only a certain
amount of time in which to experiment (I) or we don't know how long we
need to go but know how what percentage of failures we need to wait
for (II).

\begin{itemize}
\item III Random Censoring.  Let $C_1, \ldots, C_n$ be iid $G$ were
  $C_i$ is the censoring time associated with $T_i$.  We observe
  $$ Y_i = \min(T_i, C_i) $$ and
  $$ \delta_i = {\rm I} (Ti \leq C_i) =
  \left\{ \begin{array}{l l} 1 & \mbox{ if $ T_i $ not censored} \\
      0 & \mbox{ if $ T_i $ is censored}
      \end{array}
      \right . $$
\end{itemize}

Many medical examples have people enter the study at different times
by and the exit the study because of
\begin{enumerate}
\item loss of follow up (they move, etc.)
\item drop out (side effects, contrariness)
\item end of study.
\end{enumerate}

The following picture shows what might happen in a study.

\ \\
\centerline{\psfig{figure=surv1.eps,angle=-90,width=6in}}

One might also imaging left censoring, where the start date is only
know to be after a certain date.


\subsection*{Some distributions for Survival Analysis}


\paragraph{Exponential}

This is the basic model with constant risk, $\lambda(t) = \lambda > 0
$.  Thus
$$ \int_0^t{\lambda(u) du} = \lambda \times t $$ 
and so
\begin{eqnarray*}
  S(t) &=& \exp{ \int_0^t{\lambda(u) du}} = e^{- \lambda t} \\
  F(t) &=& 1 - e^{- \lambda t} \\
  f(t) &=& \lambda e^{- \lambda t} 
\end{eqnarray*}
and thus
$ E(T) = 1/\lambda$ and $\mbox{Var}(T) = 1 / {\lambda}^2$.
Here is a plot of the density.  

\psfig{figure=expplot.eps,angle=-90,width=4in}

A reasonable thing to do might be to plot $\log{ (\hat{S}(t))}$
against $t$.  It should be a straight line.


\paragraph{Gamma}

One possible generalization of the exponential.
$$
f(t) = \frac{\lambda^{\alpha}}{\Gamma(\alpha)} t^{\alpha -1}
e^{\lambda t} \quad \alpha>0 \quad \lambda>0$$

\begin{eqnarray*}
  F(t) &=& \mbox{Incomplete Gamma Integral, no closed form} \\
  \lambda(t) &=& 1 - e^{- \lambda t} \\
\end{eqnarray*}

$ E(T) = \alpha/\lambda$ and $\mbox{Var}(T) = \alpha / {\lambda}^2$.
Here is a plot of the density, for $\alpha = 0.5, 1, 2, 3 $

\psfig{figure=gamplot.eps,angle=-90,width=4in}

You really already know a lot about the Gamma distribution because
$$
\chi^2_k = \mbox{Gamma}(k/2, 1/2) $$
that is, $\alpha = k/2 $ and $\lambda = 0.5 $.


\paragraph{Weibull}

Another generalization of the exponential.  

$$ S(t) = e^{ - (\lambda t) ^ \alpha} \quad \alpha >0 \quad \lambda >
0 $$

This can be forced to look a LOT like a gamma distribution.  The
Weibull is very important, but we will not dwell on it.


With complete (i.e., uncensored) data estimation of any of these
distributions is relatively easy.  Just write out the likelihood
equations and differentiate.  For some of the models there are even
closed-form estimates.  We will defer looking at the parametric ML
estimates until we have looked at a non-parametric estimate of the
survivor function.


\subsection*{Estimation of the Survivor Function}
(Without parametric models).

Remember what the survivor function is, just the probability of
surviving to be older than time $t$.  A convenient and familiar object
is the sample cumulative distribution function (cdf), written as
$$
\bar{F}_n(x) = \frac{ \mbox{ \# of obs} \leq x}{n} $$

Estimates of the sample survival function are obviously of the form, $
1 - \bar{F}_n(x) $.  In the analysis of survival data it is often
useful to summarize the survival experience in terms of 
$\bar{F}_n(x)$ or $\bar{S}_n(x)$.  (Why?  Well it might direct us to
an appropriate parametric model, or it might be interesting in its own
right).  When the data is censored we need to think of something other
than the simple cdf.


\paragraph{Life Tables}

The classical method of estimating $S(t)$ is the so-called actuarial
method which depends on the life table.

Let time be partitioned into a fixed sequence of intervals, $I_1,
\ldots, I_k$ mostly of equal length.  The intervals should be open at
the left and closed at the right.  

Consider the line

\centerline{\psfig{figure=surv2.eps,angle=-90,width=6in}}

with
\begin{eqnarray*}
  n_i &=& \mbox{\# alive at beginning of } I_i \\
  d_i &=& \mbox{\# died in } I_i \\
  l_i &=& \mbox{\# lost to follow-up in } I_i \\
  w_i &=& \mbox{\# withdrawn in } I_i \\
  p_i &=& P(\mbox{ surviving through } I_i | \mbox{ alive at start of } I_i
\end{eqnarray*}

As an example, look at the data below, where we would like to estimate
$S(5)$ say.

\begin{center}
  \begin{tabular}{r r r r r}
\hline
Years after diag & Alive at start & Died & Lost to Followup & Withdrawn\\
\hline
0--1 & 126 & 47 & 4 & 15 \\
1--2 & 60 & 5 & 6 & 11 \\
2--3 & 38 & 2 & --- & 15 \\
3--4 & 21 & 2 & 2 & 7 \\
4--5 & 10 & --- & --- & 6 
  \end{tabular}
\end{center}


\paragraph{Reduced Sample Method}

To estimate $S(\tau_k)$ we only use those subjects who are at risk for
the entire period $(0, \tau_k]$.  That is then
$$
  \begin{array}{r c l}
  n &=& n_1 - \sum_{i=1}^{k}{l_i} - \sum_{i=1}^{k}{w_i} \\
  d &=& \sum_{i=1}^{k}{d_i} \\
\mbox{and  } S(\tau_k)& =& 1 - d/n
  \end{array}
\qquad
\begin{array}{r c l}
  n &=& 126 - 12 - 54 = 60 \\
  d &=& 56 \\
    & = & 1 - \frac{56}{60} = 0.078 
\end{array}
$$

What (if anything) is wrong with this?  It throws away information,
all the iformation in $l_i$ and $w_i$, and it is thus an inefficient
(and biased) estimate.  


\paragraph{Actuarial Method}

The survival probability $S(\tau_k)$ can be broken into a product of
probabilities, 
\begin{eqnarray*}
  S(\tau_k) &=& P(T > \tau_k) \\
            &=& P(T > \tau_1) \cdot P(T > \tau_1| T > \tau_2) \cdots
            P(T > \tau_k| T > \tau_{k-1}) \\ 
            &=& p_1 \cdot p_2 \cdot p_3 \cdots p_k \\
\mbox{where} \\
p_i &=& P(T > \tau_i| T > \tau_{i-1})
\end{eqnarray*}

The actuarial method estimates each $p_i$ separately and multiplies to
get $S(\tau_k)$.  To estimate $p_i$ one could use
$$ 1 - \frac{d_i}{n_i}$$, but this ignores losses and withdrawals in
interval $I_i$.  A sensible(?) approach might be to assume that on
average those that were lost or withdrew were at risk for $1/2$ of the
interval.  Thus one might define the effective sample size as
$$
n_i' = n_i - 1/2 \cdot (l_i - w_i) $$
with
\begin{eqnarray*}
 \hat{q}' &=& \frac{d_i}{n_i'} \\
    \hat{p}_i &=& 1 - q_i \\
\mbox{and   }  \hat{S}(\tau_k) &=& \prod_{i=1}^{k}{(\hat{p}_i)}
\end{eqnarray*}

For the above data, this yields,
$$ \hat{S}(\mbox{5 years}) = 0.44 $$
(You work this out!)


\paragraph{Variance of $\hat{S}(\tau_k)$}

By assuming that 
$$
n_i' \hat{p_i} \approx \mbox{Binomial}(n_i', p_i) $$
and some manipulations an estimate of the variance of the estimate can
be found as
$$
\widehat{\mbox{Var}}(\hat{S}(\tau_k)) = \hat{S}^2(\tau_k)
                  \sum_{i=1}^{k}{\frac{d_i}{n_i' (n_i' - d_i)}} $$

This is known as Greenwood's formula (named after Major Greenwood).




