\subsection*{The Kaplan-Meier Product Limit Estimate}

This is very similar to the actuarial method we just saw, except that
the intervals $I_i$ are not fixed in advance.  In fact the intervals
are determined by when an event occurs.  So a revised picture might
be, 

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

Recall that we observe 
$ 
\begin{array}{r c l}
  Y_1, \ldots , Y_n & & \mbox{survival times, and} \\
  \delta_1, \ldots , \delta_n & & \left (
  \begin{array}{r c l}
1 & & \mbox{censored} \\
0 & & \mbox{not censored} \\
  \end{array} \right ) \\
\end{array}
$

Order the $Y_i$ into
$
\begin{array}{r}
  Y_{(1)}, \ldots , Y_{(n)} \\
  \delta_{(1)}, \ldots , \delta_{(n)} \\
\end{array}
$ and carry only the $\delta_i$

Let ${\cal R}(t)$ be the risk set at $t$ and
\begin{eqnarray*}
  n_i &=& \mbox{\# in ${\cal R}(Y_{(i)}) = $ \# alive just before } Y_{(i)}\\
  d_i &=& \mbox{\#died at }  Y_{(i)} \\
  p_i \mbox{ and } q_i &=& \mbox{same as before---conditional
  probability of surviving}
\end{eqnarray*}

Then
$$
\hat{p}_i = \left \{ 
\begin{array}{c c l}
1 - 1/n_i & \mbox{if} & \delta_{(i)} = 1 \\
1 & \mbox{if} & \delta_{(i)} = 0 \\
\end{array}
\right . $$

Thus, without ties, the PL estimate is 
\begin{eqnarray*}
\hat{S}(t) &=& \prod_{y(i) \leq t}{ \hat{p}_i} = \prod{(1 -
    1/n_i)^{\delta(i)}} \\
 &= & \prod{(1 - \frac{1}{(n - i + 1)})^{\delta(i)}} 
\end{eqnarray*}


\paragraph{Notes}
\begin{enumerate}
\item Tied Uncensored observations---do the obvious thing which is
 $ p_i = 1 - \frac{d}{m} $ where $d$ deaths occur.
\item Censored and uncensored tied data---consider the uncensored data
  to be just before the censored data.  (Why in this order?)
\end{enumerate}


\subsection*{Two Views of the PL estimate as an MLE}

The PL estimate would not be nearly so popular if it could not be
interpreted in other ways.

Let $t_1,\ldots,t_k$ be the observed failure times in a sample of size
$n_0$ with survivor function, $S(t)$.  Suppose $d_j$ items fail at
time $t_j$ and $m_j$ items are censored in $[t_j, t_{j+1})$ at known
times $t_{j_l}$.

Let $n_j = (m_j + d_j) + \cdots + (m_k + d_k) $ be the number of items
at risk at time just before $t_j$.  The probability of failure at
$t_j$ is
$$
 S(t_j) - S(t_j + 0) $$
We assume that the contribution to the likelihood of a survival time
 censored at $t_{j_l}$ is
$P(T > t_{j_l}) = S(t_{j_l} + 0) $.

Thus,
$$
{\cal L} = \prod_{j=0}^{k} \left\{
{\left[  S(t_j) - S(t_j + 0) \right]}^{d_j} \cdot
  \prod_{l=1}^{m_j} S(t_{j_l} + 0) \right\} $$

Given the data, this can be viewed as a likelihood function on the
space of all survivor functions, $S(t)$.  There are problems with this
definition of the MLE, but the intuition is good.   

Clearly to maximize ${\cal L}$ one needs and $S$ which is {\em
  discontinuous\/} at the observed failure times---otherwise ${\cal
  L}$ is zero.   

Further, ${\cal L}$ is maximized when 
$S(t_{j_l}+0) = S(t_{j}+0)$

An estimate $\hat{S}(t)$ is then a discrete survivor function with
hazard components $\hat{\lambda}_i$ at $t_i$.

Thus,
$$
\hat{S}(t_j) = \prod_{l=1}^{j-1} (1 - \hat{\lambda}_l)
$$
and
$$
\hat{S}(t_j + 0) = \prod_{l=1}^{j} (1 - \hat{\lambda}_l)
$$
where the $\lambda$ are chosen to maximize the function
\begin{eqnarray*}
 & &\prod_{j=1}^{k}\left\{ \left[ 
  \prod_{l=1}^{j-1} (1 - \lambda_j)^{d_j} \right]
       \lambda_j^{d_j} \left[ 
            \prod_{l=1}^{j} (1 - \lambda_l)^{m_j} \right] \right \} \\  
 & = & \prod_{j=1}^{k} \lambda_j^{d_j} ( 1 - \lambda_j)^{n_j - d_j}
\end{eqnarray*}
because
$$
S(t_j) - S(t_j + 0) = \prod_{l=1}^{j-1} \lambda_j (1 - \lambda_j)
$$

Obviously the maximum is achieved when
$ \hat{\lambda_j} = \frac{d_j}{n_j} $ and the product limit estimate
is
$$
\hat{F}(t) = \prod_{j \atop{t_j < t}} \left( \frac{n_j - d_j}{n_j} \right)
$$


\paragraph{Interpretation}

We are in effect making the conditional probability of failure at each
$t_j$ agree exactly with the observed conditional relative frequency
of failure at $t_j$.  

Problem:  If we assume that $S(t)$ is continuous then $\hat{S}(t)$ is
always bumpy, which might give you cause for concern.  (This can all
be fixed by appropriate technical mumbling).


\subsection*{Example}


\paragraph{No Censoring}

Assume that $n=5$ and $x = \mbox{death}$ and $o = \mbox{censor}$ with
$$
\begin{array}{c c c c c c}
 & x & x & x & x & x \\
p_i & (1 - \frac{1}{5}) &(1 - \frac{1}{4}) &(1 - \frac{1}{3}) &(1 -
\frac{1}{2}) & (1 - \frac{1}{1}) \\
\end{array}
$$

where $p_i$ is the $p(\mbox{survive to age } i | \mbox{alive before
  time} i) $
$$
\begin{array}{c c c c c c c}
  S_i & = & \frac{4}{5} & \frac{4}{5}\cdot\frac{3}{4} &
  \frac{4}{5}\cdot\frac{3}{4}\cdot\frac{2}{3} & \ldots \\
 & = & \frac{4}{5}& \frac{3}{5}& \frac{2}{5}& \frac{1}{5} & 0
 \end{array}
$$

Thus, without censoring the PL estimate is the obvious one


\paragraph{Some censored data}

Let $n=7$ and take the data
$$
\begin{array}{c c c c c c c c}
 & x & x & o &x &o & x  \\
p_i & \frac{6}{7} & \frac{5}{6} & & \frac{3}{4} & & \frac{1}{2} \\
\end{array}
$$

The PL estimate is
$$
\begin{array}{c c c c c c c c c}
  S_i & = & \frac{6}{7} & \frac{6}{7}\cdot\frac{5}{6}
  &\frac{6}{7}\cdot\frac{5}{6} & \frac{6}{7}\cdot\frac{5}{6}\cdot\frac{3}{4} &
  \ldots \\
 & = & \frac{6}{7}& \frac{5}{7}& \frac{5}{7}& \frac{15}{28} 
     & \frac{15}{28} & \frac{15}{56} & \frac{15}{56} \\
 \end{array}
$$


\section*{How does one compare two survival curves?}

That is a very good question.  One answer is to say, just use the
Splus function {\tt survdiff}, which has a description of
\begin{quote}
       Test if there is a difference between two or more survival
       curves  using  the  G-rho family of tests.  Or, test for a
       single curve against a known alternative.
\end{quote}

But, to really understand what this is all about, we need a
digression.


\section*{Digression on some Nonparametric Techniques}


\subsection*{Wilcoxen or Mann-Whitney Test}

Suppose we have samples from two distributions $F$ and $G$ and we wish
to compare the distributions.  Let $X_1, \ldots X_m$ and
$Y_1, \ldots Y_m$ denote the random samples from the two
distributions.  

One test of the difference between the distributions proceeds as
\begin{enumerate}
\item Rank all the data---i.e., mix the two samples.
\item Let $R(X_i)$ denote the ranks of the F sample in the combined
  data.
\item Form $ S = \sum_{i=1}^{m}{R(X_i)} $, the rank sum.
\end{enumerate}

Now note that $S$ will be large if the F-distribution has larger
values than the G-distribution and vice-versa.  It is easy to show
that
$$
\mbox{E}(S) = \frac{m(n+m+1)}{2} $$ and
$$
\mbox{Var}(S) = \frac{m(n+m+1)n}{2} $$ 
(This comes from assuming that $F= G $ in which case the ranks are
just sampled form a discrete uniform distribution.

For small values of $n$ and $m$ there are tables of exact tail
probabilities for $S$ (or some linear transformation of it).  

For large samples, just take a normal approximation
$$
z = \frac{ S - \mbox{E}(S)}{\sqrt{\mbox{Var}(S)}} $$

The test makes minimal assumptions and works quite well in certain
situations?  However, what is the test sensitive to?  When normality
can be assumed, the 2-sample t-test is a competitor and like most
specific purpose tools does better than the Wilcoxon test.

You can find this test as {\tt wilcox.test} in Splus.

\subsection*{Kolmogorov-Smirnov test}

The previous test is not very sensitive to differences in the
distribution shape, as opposed to solely location differences.  A test
which is sometimes useful here is the (Kolmogorov-)Smirnov test.

Let $\hat{F}$ and $\hat{G}$ be estimates of two distribution
functions.  (For censored data one might think about using the
Kaplan-Meier estimates).  The exact form of the test depends on the
alternate hypothesis (i.e., what do you want to test) so I will assume
something like
$ H_1: F \neq G$.

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


Now form $T = \max{|F - G |}$, and look up the results in a table, OR
just use the Splus function {\tt ks.gof}. 

The test is sensitive to all sorts of departures, particularly changes
in tails and different variances.
  \begin{flushright}
{\Large\bf End Digression}    
  \end{flushright}


So how does one compare two Kaplan-Meier estimates?   One of the
things that {\tt survdiff} does is a (version) of Gehan's test, which
itself is just a modification of the Wilcoxon test.

Instead of using $S$, the rank-sum, we need to develop another
equivalent statistic.

Let
$$
U(X_i, Y_j) = \left\{
\begin{array}{c c c }
+1 & \mbox{if} & X_i > Y_j \\
0 & \mbox{if} & X_i = Y_j \\
-1 & \mbox{if} & X_i < Y_j \\
\end{array} \right .
$$

Let $U = \sum_{i=1}^{m}\sum_{j=1}^{n} U_{i,j}$, which can be read as
the average number of times that $X > Y$.  (Would
$\sum_{i=1}^{m}\sum_{j=1}^{n}| U_{i,j}|$ make sense? ).

It isn't hard to see that
$$
S = \frac{m(m+n+1)}{2} = \frac{1}{2} U $$
thus
$\mbox{E}(U) = 0$ and $ \mbox{Var}(U) = \frac{n m (n+m+1)}{3}$.

For censored data, the Gehan test suggests defining $U_{i,j}$ in a
different form.

Our basic data is $(X_i, \delta_i)$ and $(Y_i, \epsilon_i)$ where
$\delta$ and $\epsilon$ are the indicators of censoring.

The null hypothesis is $ F_1 = F_2$ (distribution) and $G_1 = G_2$
(censoring), with alternate interest in $F_1 \neq F_2$.  (I am not
sure what happens when the censoring is different).

Now define
$$
U_{i,j} = \left\{
\begin{array}{c c c }
+1 & \mbox{if}
 & 
 \begin{array}{l}
  \mbox{we know that } t_i > u_j \\
  \mbox{i.e., } (x_i > y_j, \epsilon_j =1) \mbox{ or }
                (x_i = y_j, \delta_i = 0, \epsilon_i=1) \\
 \end{array} \\
0 &  \mbox{otherwise} & \\
-1 & \mbox{if}
 & 
 \begin{array}{l}
  \mbox{we know that } t_i < u_j \\
  \mbox{i.e., } (x_i < y_j, \delta_j =1) \mbox{ or }
                (x_i = y_j, \delta_i = 1, \epsilon_i=0) \\
 \end{array}
\end{array} \right .
$$
and $U = \sum_{i=1}^{m}\sum_{j=1}^{n} U_{i,j}$


It is possible to show (via some permutation theory) that
$\mbox{E}(U) = 0$ and $ \mbox{Var}(U) = \mbox{a big mess}$, but there
is a computational formula, and packages.

\section*{Digression on the Mantel-Haenszel Test}

Suppose there are two populations and each individual in each
population can have one of two characteristics.  For example, two
populations of cancer patients with different treatments who can
either survive beyond a year or die within the year.   We could write
out the table of data as, 
\begin{center}
  \begin{tabular}{l| l | l | l}
   &Dead & Alive \\ \hline
Pop1 &  $a$ & $b$ & $n_1$ \\ \hline
Pop2 &  $c$ & $d$ & $n_2$\\ \hline
     & $m_1$ & $m_2$ &$ n$
  \end{tabular}
\end{center}

Let
\begin{eqnarray*}
  p_1 &=& P(\mbox{Dead} | \mbox{Pop1}) \\
  p_2 &=& P(\mbox{Dead} | \mbox{Pop2}) \\
\end{eqnarray*}

A test of $p_1 = p_2$ can be based on (thing that we already know
about homogeneity of proportions) or the statistic
$$
\chi^2 = {\left[
  \frac{ \hat{p}_1 - \hat{p}_2}{\sqrt{\hat{p}(1-\hat{p})(\frac{1}{n_1}
      + \frac{1}{n_2})}} \right] }^2
 = \frac{n(ad - bc)^2}{n_1 n_2 m_1 m_2}
$$
where
$ \hat{p}_1 = \frac{a}{n_1} $, $ \hat{p}_2 = \frac{c}{n_2} $,  and
$ \hat{p} = \frac{m_1}{n} $.  The $\chi^2$ is approximately
distributed as a $\chi^2_1$ random variable. 

The intuition for this test comes (again) from the 2-sample t-test,
but the test itself is derived from a hypergeometric distribution for
the contingency table.

$$
P(A=a) = \frac{{n_1 \choose a }{n_2 \choose{ m_1 -a}}}{{n \choose m_1}} $$

This test assume that the $n$ and $m$ margins are BOTH fixed in
advance.  

Under the hypergeometric distribution (which is our null)
$$
\mbox{E}(A) = \frac{n_1 m_1}{n} $$
and 
$$
\mbox{Var}(A) = \frac{n_1 n_2 m_1 m_2}{n^2 (n - 1)} $$
whence
$$ ad-bc = n(a - \mbox{E}(A)$$ and
$$
n_1 n_2  m_1  m_2 = n^2 (n-1) \mbox{Var}(A) $$
thus
$$
\chi^2 = {\frac{n}{n-1}\left[
  \frac{a - \mbox{E}(A)}{ \sqrt{ \mbox{Var}(A)}}
   \right]}^2
$$

For the example above, the hypergeometric assumption is suspect.  If I
collected the data by finding $m_1$ deaths and $m_2$ survivors and
then recorded their populations, then it is probably OK.  On the other
hand, if I just start out with two populations and watch for deaths a
binomial (or product binomial) seem more reasonable.  But the
particular form of the $\chi^2$ test is (I think) just the same as the
binomial one would have been, modulo a fudge about $n$ or $n-1$.


{\small\em The real reason that the hypergeometric is used here is
  based on the statistical notion of ancilarity\/}


\subsection*{Continue Digression, a sequence of $2 \times 2$ tables}
The same setup as before, but now there are $K$ hospitals.  Because
there may be differences between hospitals we don't wish to combine all
the $2 \times 2$ tables.  Rather we would like to test
$$
p_1^1 = p_2^1, p_1^2 = p_2^2, \ldots , p_1^k = p_2^k, $$
i.e., equality in each population, but not overall equality.

The data are now
\begin{center}
  \begin{tabular}{l| l | l | l}
   &Dead & Alive \\ \hline
Pop1 &  $a_1$ & $b_1$ & $n_11$ \\ \hline
Pop2 &  $c_1$ & $d_1$ & $n_12$\\ \hline
     & $m_11$ & $m_12$ &$ n_1$
  \end{tabular} \qquad \ldots \qquad
  \begin{tabular}{l| l | l | l}
   &Dead & Alive \\ \hline
Pop1 &  $a_k$ & $b_k$ & $n_k1$ \\ \hline
Pop2 &  $c_k$ & $d_k$ & $n_k2$\\ \hline
     & $m_k1$ & $m_k2$ &$ n_k$
  \end{tabular}
\end{center}

The Mantel-Haenszel statistic for this problem is 
$$
\mbox{MH} = \frac{ \sum_{i=1}^{k} ( a_i - \mbox{E}(A_i))}
                 {\sqrt{\sum_{i=1}^{k} (\mbox{Var}(A_i))}} $$

If the tables are independent then MH tends towards a $N(0,1)$ when
either $k$ is fixed and $n_i \rightarrow \infty$ or $ k \rightarrow
\infty$ and the tables are iid.

Of course we could use our contingency table techniques to answer
these questions---but you end up with slightly different tests.

When used in survival analysis, neither of these conditions for
asymptotic normality holds---but the result still does.  

  \begin{flushright}
{\Large\bf End Digression}    
  \end{flushright}


\section*{Kaplan-Meier in Survival Analysis}

In survival analysis this test is used in the following manner.  Let
$ 
\begin{array}{r }
  X_1, \ldots , X_n \\
  \delta_1, \ldots , \delta_n \\
\end{array}
$
be one population and
$ 
\begin{array}{r }
  Y_1, \ldots , Y_n \\
  \epsilon_1, \ldots , \epsilon_n \\
\end{array}
$
be the other.   Consider $Z_{(1)}, \ldots, Z_{(n+m)}$ to be the
combined ordered populations.  For each $Z_{(i)}$ which corresponds to a
death (in either population) construct a $2 \times 2$ table and
compute the MH test for this sequence of tables. 
(E.G., at time (i) we have
\begin{center}
  \begin{tabular}{l| l | l | l}
   &Dead & Alive \\ \hline
Pop1 &  $a$ & $b$ & $n_1$ \\ \hline
Pop2 &  $c$ & $d$ & $n_2$\\ \hline
  \end{tabular}  
\end{center}
and at time (i$+$1) we have
\begin{center}
  \begin{tabular}{l| l | l | l}
   &Dead & Alive \\ \hline
Pop1 &  $a+1$ & $b-1$ & $n_1$ \\ \hline
Pop2 &  $c$ & $d$ & $n_2$\\ \hline
  \end{tabular}  
\end{center}
 etc.,)

In other words, test $F_1 = F_2$ (which is the same as $S_1 = S_2$) by
using the MH test.   (Independence of tables is violated here, but
asymptotic normality of the test can still be pushed through).



\subsection*{What next}

Before proceeding to other (maybe more complicated) matters, we need
to ask if any of these ideas are useful or used in practice.  The best
I can do is present a section of the book ``Biostatistics Casebook''
by Miller, Efron, et. al. which shows a very simple yet convincing
application.




Next topic, {\bf Spring Break}

% Local Variables: 
% mode: latex
% eval: (LaTeX-math-mode t)
% eval: (my-latex-mode)
% eval: (read-abbrev-file "~/.abbrev_defs")
% abbrev-mode: t
% TeX-master: "class"
% End: 


