\documentclass{article}
\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}

%%%%
%
% Defs for the footer with the date
%
%%%%

\makeatletter
\newcount\m \newcount\n
\def\twodigits#1{\ifnum #1<10 0\fi \number#1}
\def\hours{\n=\time \divide\n 60
  \m=-\n \multiply\m 60 \advance\m \time
  \twodigits\n:\twodigits\m}
%
\def\draftfooter{\jobname{} \today{} \hours}
\def\footertext{\draftfooter}
\def\ps@headings{%
  \def\@oddfoot{\hfil \footertext \hfil}
  \let\@evenfoot\@oddfoot
  \def\@oddhead{{\slshape\rightmark}\hfil\thepage}%
  \let\@mkboth\markboth
  \def\chaptermark##1{%
    \markright {\MakeUppercase{%
      \ifnum \c@secnumdepth >\m@ne
        \if@mainmatter
          \@chapapp\ \thechapter. \ %
        \fi
      \fi
      ##1}}}}
\ps@headings
\makeatother
%
\def\shorttoday{%
  \number\year%
  \ifcase\month\or jan\or feb\or mar\or apr\or may\or
     jun\or jul\or aug\or sep\or oct\or nov\or dec\fi
  \twodigits\day}
\def\edrxnotesfooter#1{{\tiny #1 (typeset \shorttoday{} \hours)
  --- working draft}}
\def\edrxnotes#1{\gdef\footertext{\edrxnotesfooter{#1}}}
\def\noedrxfooter{\gdef\footertext{}}
\edrxnotes{Zhang/Asgharian/Wolfson}

%%%%
%
% All the other defs
%
%%%%

\usepackage{graphicx}
\def\epsscale{0.8}
\def\eps#1{\includegraphics[scale=\epsscale]{zawfigs/#1.eps}}

% These are used in some tables, but they need fixing
%
\def\twotextlines#1#2{\vtop{\hbox{#1}\hbox{#2}}}
\def\twotextlinesc#1#2{\vtop{\hbox{\hfil#1\hfil}\hbox{\hfil#2\hfil}}}

% Of the following two blocks, the one that comes last will override
% the defs in the previous one - the last one is the one that counts.

% Defs for the ``serious'' version
%
\def\notepbreak#1{}
\def\myeqref#1{\ref{#1}}
\def\myfigure#1#2{#2}	% \myfigure{text}{\eps{...}}
%
% (eev "cd ~/LATEX/; latex wolfson.tex; cp wolfson.dvi zaw.dvi")

% Defs for the draft version
%
\def\notepbreak#1{${\blacktriangleleft}#1{\blacktriangleright}$}
\def\myeqref#1{\ref{#1}{\scriptsize(#1)}}
\def\myfigure#1#2{#1}	% \myfigure{text}{\eps{...}}
%
% (eev "cd ~/LATEX/; latex wolfson.tex && pdflatex wolfson.tex && cp -v wolfson.dvi zawdraft.dvi && mv -v wolfson.pdf zawdraft.pdf")



%%%%
%
% The document itself
%
%%%%

\begin{document}

% Started 2002may09 16:54

% Cover page
% (find-node "(latex)\\author")

\title{Estimating survival \\
  from partially observed data}
\author{Xun Zhang \\
        Division of Clinical Epidemiology, Montreal General Hospital
  \and  Masoud Asgharian \\
        Department of Mathematics and Statistics, McGill University, Montreal
  \and  David B.\ Wolfson \\
        Department of Mathematics and Statistics, McGill University, Montreal}
\date{Working draft, \today}
\maketitle

\section*{Abstract:}

Often, in cross-sectional-follow-up studies, survival data are
obtained from prevalent cases only. This sampling mechanism introduces
lenght-bias. An added difficulty is that in some cases the times of
the onset cannot be ascertained or are recorded with great
uncertainty. Such was the situation in the Canadian Study of Health
and Aging (CSHA), an ongoing nation-wide study of dementia conducted
by Health Canada. This paper proposes methods to estimate the survival
function nonparametrically, when the data are length-biased and only
partially observed. By using the ``forward recurrence times'' only, we
suggest how one can overcome the difficulty caused by missing onset
times, while by using the ``backward recurrence times'' only, one can
avoid the cost and effort of follow-up. We illustrate our methods
through an application to data derived from the CSHA.


%%%%
%
% Section 1: Introduction
%
%%%%

%
% Page 1:
\notepbreak{1:}%
%

\section{Introduction}

The Canadian Study of Health and Aging (CSHA) (CSHA 1994) is one of
the largest epidemiological studies of dementia conducted so far. In
the study, over 10,000 subjects were screened for cognitive
impairment. A total of 1,132 people were identified with dementia,
including Alzheimer's disease, as CSHA-1, the initial data collection
process, which took place between February 1991 and May 1992. During
this initial phase, dates
%
% There's a "*" in the margin here that I don't understand
%
of onset of symptoms were ascertained from caregivers. Those subjects
with dementia were then followed for a further five-year period.
CSHA-2 data collection started in January 1996 and ended by May 1997.

One might term the CSHA a cross-sectional-follow-up study. The
observed survival times do not constitute a random sample from the
true survival distribution. Rather, they are {\sl length-biased} (Cox
1969); only those who have survived to the sampling time point (CSHA-1
in our example) have a chance of being selected into the study.
%
% A sentence deleted together with its figure:
% "This is illustrated in Figure 1.1.

Consider a cohort of subjects who experience two major events: the
initiating event and the terminating event. Let $X$ denote the time
interval between these two events, which is of primary interest.
Suppose that the random variable $X$ is independent of the
%
% page 1 ends, page 2 begins
\notepbreak{1|2}%
%
initiating event. This is the case in many practical simulations. For
example, in the CSHA, the initiating event corresponds to the onset of
dementia, the terminating event to the death of the individual, and
$X$ to the survival time of an individual from onset. For the majority
of dementias, survival times have remained relatively independent of
dates of onset. One focus of CSHA was to estimate the survival
distibution of subjects with an umbrella diagnosis of dementia, as
well as the survival distributions of the subgroups of those diagnosed
with possible Alzheimer's disease, probable Alzheimer's disease, and
vascular dementia.

% Moved to the end of this section:
% "While much of the ensuing discussion is motivated by the problem of
% estimating survival with dementia from a cross-sectional follow-up,
% the methods will be seen to have much broader applicability."

% A deleted block containing eq. 1.1.

Wolfson et al (2001) analyzed the CSHA data in the presence of right
censoring, and estimated the survival function of patients with
dementia adjusting for length-bias. There are several difficulties
that often arise when cross-sectional data are used fro survival
analysis, some of which arose in the analysis of the CSHA data.

1. The survival times and the censoring times are not independent; in
   fact censoring is informative, which means that most of the 
%
% page 2 ends, page 3 begins
\notepbreak{2|3}%
%
  standard results of surviavl analysis, established under the
  assumption of non-informative censoring, are not automatically
  valid. This drawback can be avoided when estimation of the survival
  function is of interest, by using a so-called ``conditional''
  estimator (Wang 1991). An added benefit of conditioning is that the
  estimator obtained is robust against non-stationarity of the onset
  times (Wang 1991). Alternatively, under stationarity (our situation)
  Asgharian et al (2001) derived the nonparametric maximum likelihood
  estimator and its asymptotic.

% I'm not absolutely sure of the order of the first two sentences here

2. There may either be complete or partial ignorance of the initiation
   date. In fact, among the 1,132 CSHA subjects, 185 of them, i.e.,
   more than 10\% of the data, have missing recorded dates of onset.
   If the full (possibly censored) observed survival times are to be
   used in the analysis, these 185 subjects must be excluded. Apart
   from a fairly large loss of information, their exclusion could be
   related to their survival. This would lead to an aditional bias,
   apart from length-bias. Further, since the onset of dementia is
   insidious, the dates of onset recorded cannot be precise.

%
% page 3 ends, page 4 begins
\notepbreak{3|4}%
%

We shall show how the difficulties laid out in 1 and 2 above may be
overcome when basing our analysis on partially observed prevalence
data. Recently, Helmer et al (2000, 2001) have proposed a method that
accounts for onset date uncertainty, in incident data.

The main objective in this paper is to estimate the survival function,
nonparametrically, in studies involving cross-sectional sampling, and
in which part of the data could be missing or unreliable. We shall
consider two ways in which our data may be incomplete: Type I data:
The initiation dates are unkown or only known with great uncertainty.
Type II data: Alternatively, assuming the initiation dates are known
it may be desired to estimate the survival function based only on the
current survival times collected as part of a cross-sectional study
without follow-up. The latter type of study, of course, avoids the
cost and effort of follow-up, and its analysis is discussed here
because of its similarity with the analysis of Type I data.

For data of type I we are required to estimate the survival function
using only the (possibly censored) ``forward recurrence times'', while
for data of type II we must use only the ``backward recurrence
times''. It may sometimes be assumed in length-biased sampling that
the (unseen) initiating events follow a stationary Poisson process on
a certain interval (Blumenthal 1967, Cox 1969). We shall refer to
analyses based on this assumption as ``unconditional methods''; this
paper focuses on ``unconditional methods'', in contrast to
``conditional methods,'' for fully observed survival times, that allow
the initiating points to be arbitrary (Wang 1991). We shall exploit
the key relationship between the forward recurrence time density and
the unbiased survival function,
%
% eq (1.2):
%
\begin{equation} \label{1.2}
  f_{FT}(x) = \frac{1-F(x)}{\mu} = \frac{S(x)}{\mu}
\end{equation}
%
where $F(x)$ is the distribution function of $X$ and $S(x) = 1 - F(x)$
is the survival function. Equation (\myeqref{1.2}) is also
well-known in the theory of renewal processes as the stationary
%
% page 4 ends, page 5 begins
\notepbreak{4|5}%
%
forward and backward recurrence time densities (Resnick 1992, chapter
3). It is convenient that the ``backward recurrence times'' have the
same distribution as the ``forward recurrence times'', owing to the
stationary assumption (Cox 1969), which means that (\myeqref{1.2}) holds when we
replace $f_{FT}(x)$ by $f_{BT}(x)$, the backward recurrence time
density.

% I have many doubts here:

It is often reasonable to assume that we have random right censoring
of the forward recurrence times. Thus the Kaplan-Meier estimator of
the forward recurrence time survivor function retains its usual
properties when based on the randomly right-censored [forward
recurrence] times that arise from type I data.

A kernel estimator that averages this Kaplan-Meier estimator may be
used to estimate $f_{FT}(x)$ and hence, finally, we may define
%
% this is labeled (1.3)
%
\begin{equation} \label{eqOneThree}
  \hat S(x) = \frac{\hat f_{FT}(x)}{\hat f_{FT}(0)},
\end{equation}
%
using (\myeqref{1.2}). In similar fashion one may estimate $S(x)$ by using only
the uncensored backward recurrence times and a kernel estimator of
$f_{BT}(x)$ based on the empirical distribution function of the
backward recurrence times.

In section 2 we introduce the terminology of this paper. Section 3
considers kernel density estimators. It is a well-known phenomenon
that if a probability density function has bounded support then kernel
density estimators are often biased at and near the boundary points of
of the support. Since we need to estimate $f_{FT}(0)$, where 0 is the
boundary point we must address this issue. A method of boundary
correction for kernel density estimation, the reflection method (Zhang
et al 1999), will be discussed in section 3. Section 4 outlines some
algorithms that we use for estimation of the survival function,
%
% page 5 ends, page 6 begins
\notepbreak{5|6}%
%
$S(x)$ as well as its mean and median. A bootstrap scheme is also
given to estimate the variances of these estimators. A simulation
study is performed in section 5, which demonstrates the plausibility
of our method. Finally in section 6 we apply our methods to the CSHA
data to estimate the survival function of patients with dementia, and
compare our results with those found by using the full survival times
after excluding those with missing dates of onset.

% This block came from page 2:

While much of the ensuing discussion is motivated by the problem of
estimating survival with dementia from a cross-sectional follow-up,
the methods will be seen to have much broader applicability.



%%%%
%
% Section 2 (pages 7..11)
%
%%%%

% Title removed?

% 2. Length-biased Sampling and the Forward Recurrence Time

% An "e"-blob:

% Length-biased sampling was systematically studied by Blumenthal
% (1967) and Cox (1969), in estimating electron tube life and lengths
% of textile fibers, respectively. The phenomenon occurs in many
% follow-up studies. In epidemiologic studies, the lifetime is defined
% as the time from initiation, such as onset of disease, to certain
% failure, such as death. Left truncation, a term by which length-bias
% is often known, arises when those individuals who failed before the
% time of recruitment are excluded from the study.

% Stopped at 18:30; worked for 1h35mins

% Restarted at 22:54

\notepbreak{7}

%%%%
%
% Section 2: Terminology and a proposed estimator
%
%%%%

\section{Terminology and a proposed estimator}

Let $U$ and $V$ be the times of meaningful initiating and terminating
events, respectively, and let
%
% eq (2.1):
%
\begin{equation} \label{2.1}
  X = V - U
\end{equation}
%
denote the ``lifetime''. Suppose that the random variable $X$ has the
cumulative distribution function $F$ with the probability density
function $f$. Let $\mu$ be the (finite) mean of $X$, which is unknown.

In a cross-sectional study, suppose that for a random variable $T$,
only an individual whose time of initiation $U' \le T$ and whose
failure time $V'$ satisfies $V' > T$, will be observed. Let
%
% eq (2.2):
%
\begin{equation} \label{2.2}
  X' = U' - V'.
\end{equation}
%
% page 7 ends, page 8 begins
\notepbreak{7|8}%
%
The random variable $X'$ is left truncated, or length-biased, and $T$
is called the left truncation time of $X$.

Let
%
% eq (2.3):
%
\begin{equation} \label{2.3}
  T_F = V' - T
\end{equation}
%
and let
%
% eq (2.4):
%
\begin{equation} \label{2.4}
  T_B = T - U'.
\end{equation}
%

Then $T_F$ and $T_B$ are called, respectively, the {\sl forward} and
{\sl backward recurrence times} in analogy with their counterparts in
renewal theory (Resnick 1992). Notice that
%
% eq (2.5):
%
\begin{equation} \label{2.5}
  X' = T_B + T_F.
\end{equation}

In cross-sectional sampling, under the stationary assumption (Wang,
1991) it is well-known that $X$ had the p.d.f., $X'$,
%
% eq (2.6):
%
\begin{equation} \label{2.6}
  f_{LB}(x) = \frac{xf(x)}{\mu},
\end{equation}
%
where $f$ is the density of $X$ (see (\myeqref{2.1})) (Cox 1969).

% Several "e"-ed things.

%
% page 8 ends, page 9 begins
\notepbreak{8|9}%
%

Let $C'$ be the censoring variable which measures the time interval
from initiation $U'$ to the time of censoring or of potential
censoring, and let
%
% eq (2.7):
%
\begin{equation} \label{2.7}
  Y' = \min(X', C').
\end{equation}
%
Further, write
%
% eq (2.8):
%
\begin{equation} \label{2.8}
  C' = T_B + D',
\end{equation}
%
where $D'$ is the forward recurrence censoring time measured from the
recruitment date. Since $X'$ (see (\myeqref{2.5})) and $C'$ have $T_B$ in
common, they are not, in general, independent. The censoring mechanism
of $X'$ is, therefore, informative, and, consequently, the Meier
estimator is not the non-parametric maximum likelihood estimator for
the survival function of the length-biased survival time (Vardi 1985).

% Two "e"-ed blocks.

% "e"-ed subsection title: "2.2 Forward Recurrence Times"

Under the assumption that the initiating events follow a stationary
Poisson process, which we shall term ``the stationarity assumption'',
we shall exploit the relationship, (\myeqref{1.2}), between the forward
recurrence time density and the unbiased survival function to estimate
latter.

%
% page 9 ends (with an "e"-ed sentence), page 10 begins
% lots of "e"-ed stuff.
\notepbreak{9|10}%
%

Since
%
% eq (2.10) (but DW crossed out the number): 
%
\begin{equation} \label{2.10}
  S(z) = \mu f_{FT}(z),   \qquad   z \le 0,
\end{equation}
%
where $S(0)=1$,
%
% eq (2.11):
%
\begin{equation} \label{2.11}
  \mu = 1 / f_{FT}(0).
\end{equation}
%
It follows that $f_{FT}(z)$ is a decreasing function.

%
% page 10 ends, 11 begins
\notepbreak{10|11}%
%

Let
%
$$
  Z = \min(T_F, D'),
$$
%
the censored forward recurrence time, and suppose that $D'$ is
independent of $T_F$, so called random censoring. Let $\delta$ be the
censoring indicator, which takes the value 1 if $T_F \le D'$, and 0
otherwise.

Suppose that $T_{F1}, T_{F2}, \ldots T_{Fn}$ are $n$ i.i.d.\ of
forward recurrence times with p.d.f. (\myeqref{1.2}) and with corresponding
censoring times $D'_1, D'_2, \ldots, D'_n$. We actually observe
%
% eq (2.12):
%
\begin{equation} \label{2.12}
  \{ (Z_i, \delta_i); i=1, 2, \ldots, n \},
\end{equation}
%
where $Z_i = \min(T_{Fi}, D'_i)$ and $\delta_i$ are the censoring
indicators. From (\myeqref{2.10}) and (\myeqref{2.11}) a natural estimator for $S(z)$ is
then
%
% eq (2.13):
%
\begin{equation} \label{2.13}
  \hat S(z) = \hat f_{FT}(z) / \hat f_{FT}(0),
\end{equation}
%
although it should be pointed out that this is not the nonparametric
maximum likelihood estimator of $S(z)$. Our task is, therefore, to
estimate the density function $f_{FT}(z)$ from the observed data,
$\{(Z_i, \delta_i); i=1,2,\ldots,n\}$, for all $z \le 0$, noting that
%
% eq (2.14):
%
\begin{equation} \label{2.14}
  \hat \mu = 1 / \hat f_{FT}(0)
\end{equation}
%
would be an estimator of the mean $\mu$ of the unbiased lifetime $X$.

% "e"-ed things.

In section 3 we discuss point estimation of $f_{FT}$ and later address
the issue of how to estimate the variance of $\hat S(z)$.

%
% page 11 ends, page 12 begins
\notepbreak{11|12}%
%



%%%%
%
% Section 3: Kernel density estimator (pages 12..17)
%
%%%%

\section{Kernel density estimator}

% the rest of page 12, page 13, and most of p. 14 are "e"-ed.
% Page 14:
\notepbreak{12|13|14}%

It $\tilde G$ is the empirical distribution function of sample $\{X_i;
i=1,2, \ldots, n\}$, which is the nonparametric maximum likelihood
estimator of $G$ (NPLME), then the kernel density estimator of the
p.d.f.\ $g$, of $G$, is given by,
%
% eq (3.4):
% 
%
\begin{equation} \label{3.4}
  \tilde g_n(x) = K_{h_n} * d \tilde G(x)
\end{equation}
$$ = \frac{1}{h_n} \int K(\frac{x-y}{h_n}) d \tilde G(y), $$
%
for given kernel $K$ (Devroye, 1984).

%
% page 14 ends, page 15 begins
\notepbreak{14|15}%
%

% Suggest to take out the "extra" commas, below.

In the presence of right censoring, however, the expression (\myeqref{3.4})
suggests the use of the Kaplan-Meier estimator of $G$ [,] to replace
the empirical distribution function [,] $\tilde G$. Let $\tilde G$
denote the Kaplan-Meier estimator of the distribution function $G$
from the censored data $\{(Y_i, \delta_i); i=1,2,\ldots,n\}$. That is,
$\tilde G = 1 - \tilde S$, where $\tilde S$ [is] the Kaplan-Meier
estimator of the survival function.

Indeed, it can be shown that the modified kernel estimator
%
% eq (3.5) (copy of (3.4)):
%
\begin{equation} \label{3.5}
  \tilde g_n(x) = K_{h_n} * d \tilde G(x)
\end{equation}
$$ = \frac{1}{h_n} \int K(\frac{x-y}{h_n}) d \tilde G(y) $$
%
is consistent for $g$ (Padgett and McNichols 1984).

Equation (\myeqref{3.5}) is the key in the sequel.

% Check if he didn't mean [0,\infty) instead.

The density of the forward recurrence times $f_{FT}$, defined by
(1.2), has support $[0, \infty]$ with 0 as the only boundary point.
Further, in order to use (\myeqref{2.13}) to estimate $S(t)$, estimation of $\mu
= f_{FT}(0)$ is crucial. We shall, without loss of generality,
therefore, focus on this special case, assuming that $g$ has support
$[0, \infty]$.

When the support of the p.d.f., $g$, has boundary points, the kernel
estimator $\tilde g_n$ defined by (\myeqref{3.4}) is often seriously biased at
and near the boundary points (Schuster 1985).
%
% page 15 ends, page 16 begins
\notepbreak{15|16}%
%
Various methods have been developed to adjust for this bias at the
expense of increasing the variance of the estimator. Among these, the
most notable one, perhaps, is the reflection method (Schuster 1985,
Zhang et al 1999): this consists of reflecting the density $g$ to the
other side of the boundary point 0 and modifying the estimator (\myeqref{3.4})
as
%
% eq (3.6):
%
\begin{equation} \label{3.6}
  \tilde g_n^R = \frac{1}{h_n} \int \left( K(\frac{x-y}{h_n}) +
   K(\frac{x+y}{h_n}) \right) \, d \tilde G(y).
\end{equation}
%
This estimator works well when $g'(0)=0$, as an estimator of $g$.
Zhang et al (1999) proposed an improved estimator which minimizes the
Mean Square Error (MSE) at and near the boundary point. Their
estimator is
%
% eq (3.7):
%
\begin{equation} \label{3.7}
  \tilde g_n^t = \frac{1}{h_n} \int \left( K(\frac{x-y}{h_n}) +
   K(\frac{x+t(y)}{h_n}) \right) \, d \tilde G(y).
\end{equation}
%
where $t$ is the transform
%
% eq (3.8):
%
\begin{equation} \label{3.8}
  t(y) = y+dy^2 + Ad^2y^3
\end{equation}
%
with $d = g'(0)/g(0)$, provided that $g$ is differentiable at $x=0$,
that $g(0)>0$, and $A > \frac{1}{3}$. Unfortunately, the choice of $A$
requires the knowledge of $g''(0)$, the second order derivative of $g$
at zero. If we have no information about $g''(0)$, we could take $A=1$
to avoid serious bias. Zhang et al (1999) suggested that although $d$
is unknown it can be estimated by
%
% eq (3.9):
%
\begin{equation} \label{3.9}
  d_n = \frac{\log \tilde g_n(h_n) - \log \tilde g_n^0 (0)}{h_n},
\end{equation}
%
% page 16 ends, page 17 begins
\notepbreak{16|17}%
%
where $\tilde g_n(h_n)$ is the usual kernel estimator of $g$ at $h_n$,
defined by (\myeqref{3.4}), and
%
% eq (3.10):
%
\begin{equation} \label{3.10}
  \tilde g_n^0(0) = \frac{1}{nbh_0}
                    \sum_{i=1}^{n} K_{(0)}(\frac{-X_i}{bh_n}).
\end{equation}
%
Here $K_{(0)}$ is a so-called endpoint kernel supported on $[-1,0]$
satisfying
%
$$
  \int_{-1}^0 K_{(0)}(t)\,dt = 1, \quad
  \int_{-1}^0 t K_{(0)}(t)\,dt = 0, \quad
  \text{and} \quad
  \int_{-1}^0 t^2 K_{(0)}(t)\,dt \neq 0,
$$
%
and
%
% eq (3.11):
%
\begin{equation} \label{3.11}
  b = \left( \frac{
     \left( \int t^2 K(t)\,dt \right)^2 \int (K_{(0)}(t))^2\,dt
   }{
     \left( \int t^2 K_{(0)}(t)\,dt \right)^2 \int (K(t))^2\,dt
   }
  \right)^{1/5} .
\end{equation}
%
The function $\tilde g_n^0(0)$ is an estimator of $g$ proposed by
Zhang and Karunamuni (1998).

Thus, in the presence of censoring, a proposed kernel density
estimator is
%
% eq (3.12):
%
\begin{equation} \label{3.12}
  \hat g_n^t(x) = \frac{1}{h_n} \int \left(
      K(\frac{x-y}{h_n}) + K(\frac{x+t_n(y)}{h_n})
   \right) \, d \hat G(y) ,
\end{equation}
%
where
%
$$
  t_n(y) = y + d_n y^2 + A d_n^2 y^3
$$
%
with the problematic $A$ an unresolved difficulty.
%
% eq (3.13):
%
\begin{equation} \label{3.13}
  d_n = \frac{
    \log \hat g_n(h_n) - \log \hat g_n^0(0)
   }{
     h_n
   } .
\end{equation}

In particular, if $g'(0)$ is known to be 0, we propose,
%
% eq (3.14):
%
\begin{equation} \label{3.14}
  \hat g_n^t(x) = \frac{1}{h_n} \int \left(
     K(\frac{x-y}{h_n}) + K(\frac{x+y}{h_n})
   \right) \, d \hat G(y) .
\end{equation}

%
% page 17 ends, page 18 begins
\notepbreak{17|18}%
%

%%%%
%
% Section 4: An algorithm and a bootstrap scheme (pages 18..21)
%
%%%%

\section{An algorithm and a bootstrap scheme}

In this section we present an algorithm for estimating the survival
function, $S(x)$, of the unbiased survival time $X$, using the
estimator $\hat S(x)$ defined by (\myeqref{2.13}). We assume that the
derivative $f'_{FT}(0)$ exists.

%%%%
%
% Subsection 4.1: Algorithm
%
%%%%

\subsection{Algorithm}

1. Given the observed data (\myeqref{2.12}) of censored forward
   recurrence times
%
$$
  \{(Z_i, \delta_i); i=1,2,\ldots,n\},
$$
%
where $Z_i = \min(T_{Fi}, D_i)$ and $\delta_i$ are the censoring
indicators, calculate the Kaplan-Meier estimator $\hat G$ of the
distribution function of the forward recurrence times $T_F$. For
convenience we assume $Z_1 \le Z_2 \le \ldots \le Z_n$. Then
%
$$
  \hat G(x) = 1 - \prod_{i=1}^{k_x} (\frac{n-i}{n-i+1})^{\delta_i} .
$$
%
Here $k_x$ is the value of $k$ such that $x \in [Z_k, Z_{k+1}].$

2. Calculate the variance $\hat \sigma^2$ of $Z_1, Z_2, \ldots, Z_n$,
   and let the bandwidth
%
$$
  h_n = \hat\sigma( \frac{15e\sqrt{2\pi}}{8n} )^{\frac{1}{5}}
       = 1.6644 \hat\sigma n^{-\frac{1}{5}} ,
$$
%
(see (Devroy 1984) for a discussion of bandwidth choice).
%
% page 18 ends, page 19 begins
\notepbreak{18|19}%
%
Although our choice of bandwidths is probably not ``optimal'' because
of the censoring, our experience with large data sets is that the
optimality of bandwidths is not important.

3. Choosing an appropriate kernel $K$ and a suitable endpoint kernel
   $K_{(0)}$, calculate the crude kernel density estimator (\myeqref{3.5})
%
$$
  \tilde g_n(x) = \frac{1}{h_n} \int K(\frac{x-y}{h_n}) d \tilde G(y)
$$
%
for $x = h_n$ and the endpoint kernel density estimate (\myeqref{3.10}) at
$x=0$:
%
$$
  \tilde g_n^0(0) = \frac{1}{nbh_0}
                    \sum_{i=1}^{n} K_{(0)}(\frac{-X_i}{bh_n}).
$$

4. Estimate the derivative of $\log(f_{FT}(x))$ at $x=0$, i.e.,
   ${f_{FT}}'(0)/f_{FT}(0)$, by (\myeqref{3.13}):
%
$$
  d_n = \frac{
    \log \hat g_n(h_n) - \log \hat g_n^0(0)
   }{
     h_n
   } .
$$

5. Estimate the density function for $x \le 0$ by (\myeqref{3.12}):
%
$$
  \hat f_{FT}(x) = \frac{1}{h_n} \int \left(
      K(\frac{x-y}{h_n}) + K(\frac{x+t_n(y)}{h_n})
   \right) \, d \hat G(y) ,
$$
%
where $t_n(y) = y + d_n y^2 + d_n^2 y^3$. Let $g_i$ denote the jump of
$\hat G$ at $Z_i$. Since $g_i = 0$ if $T_{Fi}$ is censored (i.e., when
$T_{Fi} \neq Z_i$),
%
$$ \hat f_{FT}(x)  = \frac{1}{h_n} \sum_{i=1}^n \left(
      K(\frac{x-T_{Fi}}{h_n}) + K(\frac{x+t_n(T_{Fi})}{h_n})
   \right) \, g_i
$$
%
$$ = \frac{1}{h_n} \sum_{i=1}^n \left(
      K(\frac{x-Z_i}{h_n}) + K(\frac{x+t_n(Z_i)}{h_n})
   \right) \, g_i .
$$

%
% page 19 ends, page 20 begins.
\notepbreak{19|20}%
%

6. Estimate the mean $\mu$ of the lifetime $X$ by (\myeqref{2.13}):
%
$$
  \hat\mu = 1/ \hat f_{FT}(0) .
$$

7. Combine steps 5 \& 6 to finally obtain an estimator of the survival
   function $\hat S(x)$ by (\myeqref{2.14}):
%
$$
  \hat S(x) = \hat \mu \hat f_{FT}(x) = \hat f_{FT} / \hat f_{FT}(0),
  \qquad x \le 0.
$$

Remark: $d = {f_{FT}}'(0)/f_{FT}(0) = S'(0)$. In cases where $S'(0)$
is known, or can be approximated from prior knowledge of the
distribution $X$, steps 3 \& 4 can be omitted.

%%%%
%
% Subsection 4.2: Bootstrap scheme
%
%%%%

\subsection{Bootstrap scheme}

In order to find confidence intervals of the estimators derived above,
we use the bootstrap (Davidson and Hinkley 1997). There are several
bootstrap resampling schemes with regard to survival data. The most
straightforward method is simply to resample from the observed pairs
$\{(Z_1, \delta_i); i=1,2,\ldots,n\}$ (Efron 1981). The procedure is
as follows: let $\hat H$ be the empirical distribution function on
$\mathbb{R} \times \{0,1\}$, of the observed $n$ pairs; the
distribution puts mass $1/n$ at each pair $(Z_i, \delta_i)$.

1. Draw a bootstrap sample $\{(Z_i^*, \delta_i^*); i=1,2,\ldots,n\}$
   by independently sampling $n$ times with replacement from the set
   $\{(Z_1, \delta_i); i=1,2,\ldots,n\}$. This is equivalent to
   drawing a random sample from $\hat H$.

%
% page 20 ends, page 21 begins.
\notepbreak{20|21}%
%

2. Applying steps 1--7 listed in the previous section to these
   artificial bootstrapped data, calculate accordingly $\hat \mu^*$,
   $\hat S^*(x)$, and $\hat m^*$, the estimated mean, survival
   function, and median, respectively.

3. Repeat independently steps 1 \& 2 $N$ times, obtaining $\{\hat
   \mu_k; k=1,2,\ldots,N\}$, $\{\hat S_k^*(x); k=1,2,\ldots,N\}$, and
   $\{\hat m_k^*; k=1,2,\ldots,N\}$.

4. Calculate the bootstrap variances of $\hat\mu$, $\hat S(x)$, and
   $\hat m$, i.e., the sample variances respectively from those data
   obtained in step 3. Confidence intervals can be constructed based
   on normal theory.

5. Various other types of bootstrap confidence intervals, such as
   percentile, basic bootstrap, and studentized confidence intervals,
   can also be constructed based on the data obtained in setp 3
   (Davidson and Hinkley 1997).

%
% page 21 ends, page 22 begins.
\notepbreak{21|22}%
%

%%%%
%
% Section 5: Simulations
%
%%%%

\section{Simulations}

Consider the one-parameter Gamma family
%
% eq (5.1):
%
\begin{equation} \label{5.1}
  f_a(x) = \frac{x^{a-1} e^{-x}}{\Gamma(\alpha)}, \quad x \ge 0
\end{equation}
%
with parameter $\alpha>0$. Let $X$ be a random variable with density
(5.1), mean
%
% eq (5.2):
%
\begin{equation} \label{5.2}
  \mu = \mathrm{mean}(X) = \alpha
\end{equation}
%
and survival function 
%
% eq (5.3):
%
\begin{equation} \label{5.3}
  S_\alpha(x) = \frac{1}{\Gamma(\alpha)} \int_x^\infty t^{a-1} e^{-t} \, dt,
   \quad x \ge 0,
\end{equation}
%
which cannot, in general, be expressed in closed form. The density of
the forward recurrence times induced by length-biased sampling, is, by
(\myeqref{1.2}) and (\myeqref{5.2}),
%
% eq (5.4):
%
\begin{equation} \label{5.4}
  g_a(x) = \frac{S_a(x)}{\mu} = \frac{S_a(x)}{\alpha}, \quad x \le 0.
\end{equation}

The goal of the following simulation study was to evaluate the
performance of the nonparametric estimator, $\hat S(x)$, computed from
the forward recurrence times only. If $\alpha<1$, the density $g$ has
no derivative at 0, and so our method does not apply; if $\alpha = 1$,
then (\myeqref{5.1}), (\myeqref{5.3}), and (\myeqref{5.4}) are all identical to the p.d.f.\ of the
exponential distribution. Hence, we assumed $\alpha>1$. In this case
$g'(0) = 0$, and formula (\myeqref{3.14}) can be used for the estimation of
(\myeqref{5.4}).

We began by examining the behaviour of the estimator $\hat
S_\alpha(x)$ when there is no censoring.

Let $\{Y_i; i=1,2,\ldots,n\}$ be a random sample from probability
density (\myeqref{5.4}), let 
%
$$
  h_n = \hat\sigma( \frac{15e\sqrt{2\pi}}{8n} )^{\frac{1}{5}}
       = 1.6644 \hat\sigma n^{-\frac{1}{5}} ,
$$
%
% page 22 ends, page 23 begins.
\notepbreak{22|23}%
%
the bandwidth, where ${\hat \sigma_n}^2$ is the sample variance of
$\{Y_i; i=1,2,\ldots,n\}$, and let $K$ be a kernel. Without censoring,
the kernel density estimate of (\myeqref{5.4}) is
%
% eq (5.5):
%
\begin{equation} \label{5.5}
  \hat g_\alpha(x) = \frac{1}{h_n}
      \sum_{i=1}^n \left(
          K(\frac{x-Y_i}{h_n}) + K(\frac{x+Y_i}{h_n})
        \right), \quad x \ge 0 .
\end{equation}

Our simulations consisted of three parts. In all simulations, samples
of size 500 were used. To generate samples from density (\myeqref{5.4}), Von
Neumann's rejection method (Devroye 1985, chapter 8) was employed.
Three typical kernels were examined. They are the Gaussian kernel
%
$$
  K(x) = \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}} ,
$$
%
the Epanechnikov kernel
%
$$
  K(x) = \frac{3}{4} (1-x^2), \quad -1 \le x \le 1,
$$
%
and the box kernel
%
$$
  K(x) = \frac{1}{2} \quad -1 \le x \le 1.
$$

As pointed out earlier, since the estimation of the density function
at the boundary point $x=0$, is crucial, the first part of the
simulation was to investigate the estimation of the mean $\mu$ for
different values of $\alpha$:
%
$$
  \hat \mu_\alpha = \frac{1}{\hat g_\alpha(0)},
$$
%
which also estimates $\alpha$ since $\alpha = \mu$. The summary
results are reported in Tables 5.1--3 for the three different kernels
given above. The value of $\overline{\hat \mu_\alpha}$ represents the
average over 400
%
% page 23 ends, page 24 begins.
\notepbreak{23|24}%
%
repetitions. The bias and the variance of $\hat \mu_\alpha$ were also
estimated based on the 400 repetitions.

\medskip

\centerline{Table 5.1. Estimation of $\mu=\alpha$ with the Gaussian
kernel}

$$\begin{array}{|c|c|c|c|} \hline
  \alpha & \overline{\hat \mu_\alpha} & \text{Bias} & \text{Variance} \\ \hline
  1.5 &  1.5399 & 0.0399 & 0.0239 \\ \hline
   2  &  2.0374 & 0.0374 & 0.0510 \\ \hline
   3  &  3.0543 & 0.0543 & 0.1704 \\ \hline
   5  &  5.1752 & 0.1752 & 0.8811 \\ \hline
  10  & 10.6078 & 0.6078 & 9.8478 \\ \hline
  \end{array}
$$

\medskip

\centerline{Table 5.2. Estimation of $\mu=\alpha$ with the
Epanechnikov kernel}

$$\begin{array}{|c|c|c|c|} \hline
  \alpha & \overline{\hat \mu_\alpha} & \text{Bias} & \text{Variance} \\ \hline
  1.5 &  1.5499 & 0.0499 & 0.0503 \\ \hline
   2  &  2.0573 & 0.0573 & 0.1153 \\ \hline
   3  &  3.1134 & 0.1134 & 0.4265 \\ \hline
   5  &  5.4065 & 0.4065 & 2.9459 \\ \hline
  10  & 10.9903 & 0.9903 & 21.4525 \\ \hline
  \end{array}
$$

\centerline{Table 5.3. Estimation of $\mu=\alpha$ with the box kernel}

$$\begin{array}{|c|c|c|c|} \hline
  \alpha & \overline{\hat \mu_\alpha} & \text{Bias} & \text{Variance} \\ \hline
  1.5 &  1.5567 & 0.0567 & 0.0407 \\ \hline
   2  &  2.0466 & 0.0466 & 0.0999 \\ \hline
   3  &  3.1176 & 0.1176 & 0.4067 \\ \hline
   5  &  5.3701 & 0.3701 & 2.0475 \\ \hline
  10  & 10.9541 & 0.9541 & 20.0773 \\ \hline
  \end{array}
$$

\medskip

%
% page 24 ends, page 25 begins.
\notepbreak{24|25}%
%

The second part of the simulation examined the behavior of the
nonparametric estimator for the survival function $S_\alpha(x)$:
%
$$
  \hat S_\alpha(x) = \frac{\hat g_\alpha(x)}{\hat g_\alpha(0)},
  \quad x \ge 0 .
$$
%
We plot 10 typical realizations of the estimator in Figures 5.1--3
with the three given kernels.

%
% page 25 ends, page 26 begins.
\notepbreak{25|26}%
%

\myfigure{
  \centerline{Figure 5.1. Estimates of $S_\alpha(x)$ with the Gaussian kernel}
}{
  \eps{figure51}
}

%
% page 26 ends, page 27 begins.
\notepbreak{26|27}%
%

\myfigure{
  \centerline{Figure 5.2. Estimates of $S_\alpha(x)$ with the
  Epanechnikov kernel}
}{
  \eps{figure52}
}

%
% page 27 ends, page 28 begins.
\notepbreak{27|28}%
%

\myfigure{
  \centerline{Figure 5.3. Estimates of $S_\alpha(x)$ with the box kernel}
}{
  \eps{figure53}
}

%
% page 28 ends, page 29 begins.
\notepbreak{28|29}%
%

Remarks:

1. The modified kernel density estimator works very well for densities
   with substantial mass near the boundary point (Zhang et al 1999).
   This is the case for the densities $g_\alpha$ defined by (\myeqref{5.4}) with
   small values of $\alpha$, as illustrated through Tables 5.1--3 and
   Figures 5.1--3. For larger $\alpha$, there is much less mass near
   $x=0$. We demonstrate this by plotting the densities $g_\alpha$ for
   $\alpha=2$ and $\alpha=5$ in Figure 5.4. The mean and the variance
   of $g_\alpha$ can be easily calculated: $\mathrm{E}[g_\alpha(X)] =
   (\alpha+1)/2$ and $\mathrm{Var}[g_\alpha(X)] =
   (\alpha+1)(\alpha+5)/12$. With the value of $\alpha$ increasing,
   the mean of the sample is further away from 0 and the variance
   becomes larger. As a result, the bias and the variability of our
   estimators increase as $\alpha$ increases.

\myfigure{
  \centerline{Figure 5.4. Densities $g_\alpha$ for $\alpha = 2$ \& 5}
}{
  \eps{figure54}
}

2. Among the three kernels, the Gaussian kernel performs best. There
   is little difference between the Epanechnikov kernel and the box
   kernel.

%
% page 29 ends, page 30 begins.
\notepbreak{29|30}%
%

3. For large values of $\alpha$, our estimators still perform
   reasonably well when the sample size is large. This is illustrated
   in Table 5.4 and Figure 5.5, in which samples of size 1000 and the
   Gaussian kernel was used. The values in Table 5.4 were calculated
   over 400 repetitions.

\medskip

\centerline{Table 5.4. Estimation of $\mu=\alpha$ for large $\alpha$}

$$\begin{array}{|c|c|c|c|} \hline
  \alpha & \overline{\hat \mu_\alpha} & \text{Bias} & \text{Variance} \\ \hline
   10 & 10.4434 & 0.4434 &  3.5558 \\ \hline
   15 & 15.6831 & 0.6831 & 13.0304 \\ \hline
  \end{array}
$$

\medskip

\myfigure{
  \centerline{Figure 5.5. Estimates of $S_{10}(x)$ from samples of size 1000}
}{
  \eps{figure55}
}

%
% page 30 ends, page 31 begins.
\notepbreak{30|31}%
%

% Stopped at 3:29; worked for 4h35 mins (but with a 20 mins break).
% Total until now, 1h35 + 4h35 - 20mins = 5h70-20mins = 5h50.


% Restarted 2002may13, 18:52 (panda time):

In the last part of the simulation, we investigated the eefect on the
estimator due to censoring. Samples from the density $g_\alpha$ (see
(\myeqref{5.4})), were randomly censored with an exponential random variable
with mean $\lambda$ and with censoring proportion set at about 15\%.
For instance, with $\alpha=3$ and $\lambda=12$, the censoring
probability is about 0.15. In this case, formula (\myeqref{5.5}) is no longer
valid. Instead, to estimate $g_\alpha$, the density in (\myeqref{5.4}), we use
the kernel density estimator
%
$$
  \hat g_\alpha(x) = \frac{1}{h_n n} \sum_{i=1}^n
    \left( K(\frac{x-z_i}{h_n}) + K(\frac{x+z_i}{h_n}) \right) g_i,
  \quad
  x \le 0,
$$
%
where $Z_i, i=1,2,\ldots,n$, are the censored data, and $g_i$ is the
jump of the Kaplan-Meier estimator at $Z_i$ (see (\myeqref{3.14})).

Again, samples of size 500 were taken, and estimation was based on 400
repetitions. Summary results of the estimator $\hat \mu_\alpha$ are
reported in Table 5.5, and 10 realizations of the estimator $\hat
S_\alpha$ are plotted in Figure 5.6, using the Gaussian kernel. We did
this for moderate values of $\alpha$. With large $\alpha$, the sample
size needs to be increased in order to achieve reasonble accuracy, as
noted in Remark 1 above.

\medskip

Table 5.5 Estimation of $\mu=\alpha$ in the presence of censoring

$$\begin{array}{|c|c|c|c|} \hline
  \alpha & \overline{\hat \mu_\alpha} & \text{Bias} & \text{Variance} \\ \hline
   3  &  3.0601 & 0.0601 & 0.0346 \\ \hline
   5  &  5.0863 & 0.0863 & 0.0946 \\ \hline
  \end{array}
$$


\medskip

%
% page 31 ends, page 32 begins.
\notepbreak{31|32}%
%

\myfigure{
  \centerline{Figure 5.6 Estimates of $S_\alpha(x)$ in the presence of
  censoring}
}{
  \eps{figure56}
}

%
% page 32 ends, page 33 begins.
\notepbreak{32|33}%
%


%%%%
%
% Section 6: Application to survival with dementia
%
%%%%

\section{Application to survival with dementia}

In the Canadian study of Health and Aging (see Introduction), there
were questions raised about the reliability of the recorded dates of
onset of patients, and concerns over a substantial proportion of
missing onset times. These concerns can be avoided by using the
methods described in Sections 2 and 3. In the CSHA, the forward
recurrence times are the time intervals from CSHA-1, the first phase
of the study during the year of 1991, to failure; when there is
censoring, they are the time intervals from CSHA-1 to loss of
follow-up or to CSHA-2, when the study ended in 1996. Among the 1,332
patients diagnosed with dementia, a total number of 901 patients were
classified into the categories of ``probable Alzheimer's disease'',
``possible Alzheimer's disease'', and ``vascular dementia'', with 433,
277, and 201 subjects in these categories, respectively. One of the
goals of the study was the estimation of the survival distribution
from onset, of individuals with dementia of these three types.

We first apply the methods described in Sections 2 and 3 by assuming
that the derivative of the survival function at time zero is zero.
This is essentially equivalent to the assumption that subjects cannot
fail in some very small interval after onset, which seems reasonable
by the nature of the disease.

We computed the estimator, $\hat S(t)$, of the unbiased survival
distribution, $S(t)$, from these forward recurrence times, following
the procedure proposed in 4.1, with the 
%
% page 33 ends, page 34 begins.
\notepbreak{33|34}%
%
Gaussian kernel $K = \frac{1}{\sqrt{2 \pi}} e^{-\frac{x^2}{2}}$. Then
the bootstrap scheme described in 4.2, repeated 1,000 times, was used
to calculate a 95\% pointwise confidence band for $S(t)$. The results
are summarized in Figure 6.1.

We further computed the estimated survival distributions for the three
different diagnostic categories, respectively. They are plotted in
Figure 6.2. There is little difference between the survival
distributions of the groups ``possible Alzheimer's disease'' and
``vascular dementia''. Although the estimated survival probability of
the group ``probable Alzheimer's disease'' is greater than those of
the two other groups, the difference is not statistically significant.

Table 6.1 displays the estimated mean and median survival times of
these groups, along with their respective 95\% confidence intervals,
which were computed over 1,000 bootstrap replications. Notice that all
the confidence intervals of the three groups overlap.

\centerline{Table 6.1. Estimated Mean and Median Survival Times by
Diagnosis Category}

$$\begin{array}{|c|c|c|c|c|} \hline
  \text{Diagnosis} &
    \text{Total group} &
    \twotextlines{Probable}{Alzheimer's} &
    \twotextlines{Possible}{Alzheimer's} &
    \twotextlines{Vascular}{Dementia} \\ \hline
  %
  \twotextlines{Median survival}{time (months)} &
    53 & 56 & 51 & 47 \\ \hline
  %
  \twotextlines{95\% Confidence}{interval} &
    (47,59) & (49,63) & (41,61) & (37,57) \\ \hline
  %
  \twotextlines{Mean survival}{time (months)} &
    61 & 67 & 62 & 56 \\ \hline
  %
  \twotextlines{95\% Confidence}{interval} &
    (55,67) & (58,76) & (52,72) & (46,66) \\ \hline
  \end{array}
$$

\medskip

%
% page 34 ends, page 35 begins.
\notepbreak{34|35}%
%

\myfigure{
  \centerline{Figure 6.1. Estimated Survival Curve with 95\% Confidence Band}
}{
  \eps{figure61}
}

\myfigure{
  \centerline{Figure 6.2. Comparison of Survival Curves --- By
  Diagnostic Categories}
}{
  \eps{figure62}
}

% 2002may13 19:47 (panda time): took a break
% 2002may13 21:53 (panda time): restarted

%
% page 35 ends, page 36 begins.
\notepbreak{35|36}%
%

Although the survival probability at any time between 1 to 5 years
post onset of dementia can be roughly approximated from Figures 6.1
and 6.2, we list the estimated survival probabilities at 1--5 years
post onset, with their 95\% confidence intervals, in Table 6.2.

\centerline{Table 6.2. Estimated Survival Probabilities at 1--5 Years
Post Onset of Dementia}

$$\begin{array}{|l|c|c|c|c|c|} \hline
  \text{Diagnosis} & \text{1 Year} & \text{2 Year} & \text{3 Year} &
    \text{4 Year} & \text{5 Year} \\ \hline
  \text{Total group} &
    \twotextlinesc{0.95}{(0.90, 1.00)} &
    \twotextlinesc{0.84}{(0.73, 0.94)} &
    \twotextlinesc{0.70}{(0.60, 0.80)} &
    \twotextlinesc{0.58}{(0.48, 0.68)} &
    \twotextlinesc{0.39}{(0.31, 0.48)} \\ \hline
  \twotextlines{Probable}{Alzheimer's} &
    \twotextlinesc{0.99}{(0.93, 1.00)} &
    \twotextlinesc{0.94}{(0.78, 1.00)} &
    \twotextlinesc{0.81}{(0.63, 0.98)} &
    \twotextlinesc{0.65}{(0.50, 0.81)} &
    \twotextlinesc{0.43}{(0.31, 0.56)} \\ \hline
  \twotextlines{Possible}{Alzheimer's} &
    \twotextlinesc{0.94}{(0.89, 0.99)} &
    \twotextlinesc{0.79}{(0.66, 0.91)} &
    \twotextlinesc{0.65}{(0.49, 0.81)} &
    \twotextlinesc{0.54}{(0.39, 0.69)} &
    \twotextlinesc{0.38}{(0.24, 0.51)} \\ \hline
  \twotextlines{Vascular}{Dementia} &
    \twotextlinesc{0.92}{(0.87, 0.97)} &
    \twotextlinesc{0.77}{(0.63, 0.90)} &
    \twotextlinesc{0.64}{(0.48, 0.80)} &
    \twotextlinesc{0.49}{(0.34, 0.64)} &
    \twotextlinesc{0.30}{(0.18, 0.42)} \\ \hline
  \end{array}
$$

One of the drawbacks of using the forward recurrence times is that we
are unable to estimate the distribution beyond the follow-up period.
In our example, we cannot estimate the survival distribution beyond
the 5-year period, the follow-up time period of the CSHA. We do not
have this obstacle while working with the backward recurrence times.

To illustrate our method, we also used the backward recurrence times
to estimate the survival distribution of patients with dementia. As
mentioned in the introduction, among the subjects diagnosed in the
CSHA, about 10\% of them had missing dates of onset. Out of the total
number of 901 subjects in the three categories of ``probable
Alzheimer's disease'', ``possible Alzheimer's disease'', and
``vascular dementia'', 81 had no dates of
%
% page 36 ends, page 37 begins.
\notepbreak{36|37}%
%
onset recorded. Thus there are 820 backwards recurrence times, with
396, 251, and 173 in the three categories, respectively. Table 6.3
gives the summary results of the mean and median survival times,
Figure 6.3 plots the survival distribution for the total group, along
with its pointwise 95\% confidence band, and Figure 6.4 plots the
survival distributions for the three categories.

\centerline{Table 6.3. Estimated Median Survival Times by Diagnosis
Category}
\centerline{--- Using Backwards Recurrence Times}

$$\begin{array}{|c|c|c|c|c|} \hline
  \text{Diagnosis} &
    \text{Total group} &
    \twotextlines{Probable}{Alzheimer's} &
    \twotextlines{Possible}{Alzheimer's} &
    \twotextlines{Vascular}{Dementia} \\ \hline
  %
  \twotextlines{Median survival}{time (months)} &
    55 & 58 & 68 & 57 \\ \hline
  %
  \twotextlines{95\% Confidence}{interval} &
    (50, 60) & (51, 65) & (59, 77) & (48, 66) \\ \hline
  %
  \twotextlines{Mean survival}{time (months)} &
    66 & 68 & 82 & 67 \\ \hline
  %
  \twotextlines{95\% Confidence}{interval} &
    (61,71) & (61, 75) & (73, 91) & (57, 77) \\ \hline
  \end{array}
$$

%
% page 37 ends, page 38 begins.
\notepbreak{37|38}%
%

\myfigure{
  \centerline{Figure 6.3. Estimated Survival Curve with 95\% Confidence Band}
  \centerline{--- Using Backward Recurrence Times}
}{
  \eps{figure63}
}

\myfigure{
  \centerline{Figure 6.4. Comparison of Estimated Survival Curves}
  \centerline{--- By Diagnostic Categories}
}{
  \eps{figure64}
}

%
% page 38 ends, page 39 begins.
\notepbreak{38|39}%
%

There is considerable discrepancy between Figure 6.2 and Figure 6.4.
From Figure 6.2 we see that patients in the category ``probable
Alzheimer's disease'' have better survival than that of the patients
in the other two categories; while from Figure 6.4 we find that the
patients in the category of ``possible Alzheimer's disease'' have
better survival. It is important to point out that it is not the
methodology that causes the discrepancy but the two figures are
plotted based on {\sl different} data sets. This shows that the
missing dates of onset were not {\sl randomly} missing and we conclude
that patients with missing onset may be different from the others. It
is, therefore, dangerous to simply exclude those patients with missing
onset in the study. We further illustrate this by overlaying Figure
6.1 and Figure 6.3 together in Figure 6.5, the two survival curves
estimated based on forward and backward recurrence times,
respectively, of the two different data sets. They are not comparable
methodologically, but the comparison of the survival of two slightly
different groups of patients.

If we wish to compare the results of the two methods, we should use
the same set of subjects. Figure 6.6 overlays the two estimated
survival curves by using forward and backward recurrence times,
respectively, based on the same group of patients in the CSHA --- a
total of 820 subjects in the three categories ``probable Alzheimer's
disease'', ``possible Alzheimer's disease'', and ``vascular
dementia'', excluding those who had missing dates of onset. The close
correspondence of the two curves shown in Figure 6.6 also supports our
``stationary assumption'' imposed at the very beginning, the asumption
that the initiating events follow a stationary Poisson process. Under
this assumption, the distributions of the forward and backward
recurrence times are very similar.

%
% page 39 ends, page 40 begins.
\notepbreak{39|40}%
%

\myfigure{
  \centerline{Figure 6.5. Comparison of Survival Curves}
  \centerline{--- Forward and Backward Recurrence Times}
}{
  \eps{figure65}
}

\medskip

\myfigure{
  \centerline{Figure 6.6. Comparison of Survival Curves}
  \centerline{--- Forward and Backward Recurrence Times, Same Subjects}
}{
  \eps{figure66}
}

%
% page 40 ends, page 41 begins.
\notepbreak{40|41}%
%

Notice that in Figure 6.6 there is a moderate separation of the two
estimated survival curves beyond 30 months. Censoring of the forward
recurrence times might account for the discrepancy.

Comparing our results, especially Table 6.1, with those found by
Wolfson et al (2000) by using the full survival times after excluding
those with missing dates of onset, our estimated median survival times
are much longer, though all the corresponding 95\% confidence
intervals overlap. Their estimated median survival time for the total
group was 40 months with a 95\% confidence interval of (32, 48), and
their estimated medians for the three categories, ``probable
Alzheimer's disease'' ``possible Alzheimer's disease'', and ``vascular
dementia'', were 38 (18, 58) months, 42 (28, 54) months, and 40 (28,
52) months respectively.

One possible explanation for this discrepancy could be that the
asumption we made at the beginning of this section, that the
derivative of the survival function at time zero is zero, is not
accurate. To avoid this assumption, we have to estimate the derivative
of the logarithm of the density of the forward recurrence times at
zero, $f'_{FT}(0)/f_{FT}(0)$, following steps 3 and 4 in 4.1, In doing
so we used the Gaussian kernel and the following endpoint kernel
(recall (\myeqref{3.10})) proposed by Zhang and Karunamuni (1998):
%
$$
  K_{(0)}(t) = 12 (1+t) (0.5+t) I_{[-1,0]}(t),
$$
%
where $I_{[-1,0]}(t)$ denotes the indicator function of the interval
$[-1, 0]$.

Summary results are shown in Table 6.4, which are very close to those
estimates calculated by Wolfson et al (2000), although our methods are
different from theirs since they used the full survival data. We also
plot the estimated survival curve for the total
%
% page 41 ends, page 42 begins.
\notepbreak{41|42}%
%
group along with a 95\% confidence band in Figure 6.7, as well as the
estimated survival curves for the three categories, in Figure 6.8.

% 2002may13 22:43 (panda time): took another break :)
% 2002may14  0:27 (panda time): restarted

\centerline{Table 6.4. Estimated Mean and Survival Times by Diagnosis Category}

$$\begin{array}{|c|c|c|c|c|} \hline
  \text{Diagnosis} &
    \text{Total group} &
    \twotextlines{Probable}{Alzheimer's} &
    \twotextlines{Possible}{Alzheimer's} &
    \twotextlines{Vascular}{Dementia} \\ \hline
  %
  \twotextlines{Median survival}{time (months)} &
    44 & 51 & 41 & 41 \\ \hline
  %
  \twotextlines{95\% Confidence}{interval} &
    (36, 52) & (39, 63) & (29, 53) & (31, 51) \\ \hline
  %
  \twotextlines{Mean survival}{time (months)} &
    50 & 55 & 50 & 48 \\ \hline
  %
  \twotextlines{95\% Confidence}{interval} &
    (43, 57) & (37, 73) & (40, 60) & (40, 56) \\ \hline
  \end{array}
$$

\myfigure{
  \centerline{Figure 6.7. Estimated Survival Curve with 95\% Confidence Band}
}{
  \eps{figure67}
}

%
% page 42 ends, page 43 begins.
\notepbreak{42|43}%
%

\myfigure{
  \centerline{Figure 6.8. Comparison of Survival Curves}
  \centerline{--- By Diagnostic Categories}
}{
  \eps{figure68}
}

\medskip

Observe that in Figue 6.8, the survival of patients in the category of
``probable Alzheimer's disease'' is much better than the other two
groups, compared to what we have seen in Figure 6.2 and 6.4. There is
no contradiction here since we get much wider 95\% confidence
intervals this time This is mainly due to the bias introduced by the
estimation of $f'_{FT}(0)/f_{FT}(0)$ and the lack of knowledge of
$f''_{FT}(0)$ (see (\myeqref{3.8})). This is the price we pay for not using the
full information (full survival times); sometimes we do not have a
choice. Notice that all the confidence intervals calculated in Table
6.4 overlap the corresponding confidence intervals presented in Table
6.1. This also suggests that the assumption, that the derivative of
the survival function at time zero is zero, we made at the beginning,
may not be unreasonable.

%
% page 43 ends, page 44 begins.
\notepbreak{43|44}%
%

The same procedures can be also applied to the backward recurrence
times. The estimated survival curves are plotted in Figures 6.9.
Confidence intervals can be obtained by bootstrapping.

\myfigure{
  \centerline{Figure 6.9. Estimated Survival Curves}
  \centerline{--- By Diagnostic Categories}
}{
  \eps{figure69}
}

%%%%
%
% Section 7: Conclusion
%
%%%%

%
% page 44 ends, page 45 begins.
\notepbreak{44|45}%
%

\section{Conclusion}

The methods proposed in this thesis represent a first attempt to
estimate a full survival distribution nonparametrically based only on
forward or backward recurrence times. They allow us to estimate the
full survival distribution without following up subjects. This would
surely save cost and time. Moreover, backward recurrence times cannot
be censored. Alternatively, by using only the forward recurrence times
our methods allow us to estimate the survival distribution without
knowing the onset times. For many situations this is the case. Finally
our methods may be useful even if we had full data. For we could
estimate the survival function of those with missing onset times and
compare it with the survival function of those with observed onset
times. Differences between these two estimated survival functions may
indicate that those with missing onset times are not missing at
random.

Our methods work very well with large data sets with substantial mass
near the origin. The estimator (\myeqref{2.13}), $\hat S(t) = \hat f_{FT}(t) /
\hat f_{FT}(0)$ (or $\hat S(t) = \hat f_{BT}(t) / \hat f_{BT}(0)$), is
very sensitive to the estimation of the value of the density of the
forward or backward recurrence times at the origin. It is always
difficult to estimate a density function at a boundary point. When
there are few observations near the origin, serious bias could occur.
In particular, without the knowledge of the derivative of the survival
function at time zero, we have to, in addition, estimate
$f'_{FT}(0)/f_{FT}(0)$ (or its counterpart with the backward
recurrence times). This could
%
% page 45 ends, page 46 begins.
\notepbreak{45|46}%
%
introduce further serious bias. Therefore, our methods should be used
with caution when there are not enough data points near the origin.

%
% page 46 ends, page 47 begins.
\notepbreak{46|47}%
%

%%%%
%
% References
%
%%%%

\section*{References}

{\bf References}

Asgharian, M., M'Lan, C.\ E., and Wolfson, D.\ B.\ (2000).
Length-biased sampling with right censoring: an unconditional
approach. To appear.

Blumenthal, S.\ (1967). Proportional sampling in life length studies.
{\sl Technometrics} {\bf 9}: 205--218.

Colett, D.\ (1994). {\sl Modelling Survival Data in Medical Research.}
Chapman \& Hall.

Cox, D.\ R.\ (1969). Some Sampling Problems in Technology. In {\sl New
Developments in Survey Sampling.} Edited by Johnson \& Smith. John
Wiley \& Sons.

Cox, D.\ R.\ and Lewis, P.\ A.\ W.\ (1966). {\sl The Statistical
Analysis of Series of Events.} John Wiley \& Sons.

CSHA (1994). Canadian Study of Health and Aging: Study methods and
prevalence of dementia. {\sl Canadian Medical Association Journals}
1994; {\bf 150}: 899--913.

Davison, A.\ C.\ and Hinkley, D.\ V.\ (1997). {\sl Bootstrap Methods
and Their Applications.} Cambridge University Press.

Devroye, L.\ and Györgi, L.\ (1985). {\sl Nonparametric Density
Estimation.} John Wiley \& Sons.

Efron, B.\ (1981). Censored data and the bootstrap. {\sl Journal of
the American Statistical Association} {\bf 76}: 312--319.

Folland, G.\ B.\ (1984). {\sl Real Analysis.} John Wiley \& Sons.

%
% page 47 ends, page 48 begins.
\notepbreak{47|48}%
%

Kalbfleisch, J.\ D.\ and Prentice, R.\ L.\ (1980). {\sl The
Statistical Analysis of Failure Time Data.} John Wiley and Sons.

McFadden, J.\ A.\ (1962). On the lengths of intervals in a stationary
point process. {\sl Journal of the Royal Statistical Society. Series
B} {\bf 24}: 364--382.

Padgett, W.\ T.\ and McNichols, D.\ T.\ (1984). Nonparametric density
estimation from censored data. {\sl Communications in Statistics:
Theory and Methods} {\bf 13}: 1581--1661.

Parzen, E.\ (1962). On estimation of a probability density function
and mode. {\sl Annals of Mathematical Statistics} {\bf 33}:
1065--1076.

Resnick, S.\ I.\ (1992). {\sl Adventures in Stochastic Processes.}
Birkhauser.

Rosenblatt, M.\ (1956). Remarks on some nonparametric estimates of a
density function. {\sl Annals of Mathematical Statistics} {\bf 27}:
832--837.

Schuster, E.\ F.\ (1985). Incorporating support constraints into
nonparametric estimators of densities. {\sl Communications in
Statistics: Theory and Methods} {\bf 14}: 1123--1136.

Stein, E.\ M.\ (1970). {\sl Singular Integrals and Differentiability
Properties of Functions.} Princeton University Press.

Vardi, Y.\ (1985). Empirical distributions in selection bias model.
{\sl Annals of Statistics} {\bf 13}: 178--205.

Wang M.-C.\ (1991). Nonparametric estimation from cross-sectional
survival data. {\sl Journal of the American Statistical Association}
{\bf 86}: 130--143.

Wolfson, C., Wolfson, D.\ B., et al (2000). A reevaluation of the
duration of survival after the onset of dementia. {\sl The New England
Journal of Medicine} {\bf 344}: 1111-1116.

Zhang, S., Karunamuni, R.\ J.\ (1998). On kernel density estimation
near endpoints. {\sl Journal of Statistical Planning and Inference}
{\bf 70}: 301--316.

%
% page 48 ends, page 49 begins.
\notepbreak{48|49}%
%

Zhang, S., Karunamuni, R.\ J., and Jones, M.\ C.\ (1999). An improved
estimator of the density function at the boundary. {\sl Journal of the
American Statistical Association} {\bf 94}: 1231--1241.

\end{document}


%%%%
%
% Other stuff
%
%%%%

% Time logs and totals

% Started 2002may09 16:54
% Stopped at 18:30; worked for 1h35mins
% Restarted at 22:54
% Stopped at 3:29; worked for 4h35 mins (but with a 20 mins break).
% Total until now, 1h35 + 4h35 - 20mins = 5h70-20mins = 5h50.

% Restarted 2002may13, 18:52 (panda time):
% 2002may13 19:47 (panda time): took a break
% 2002may13 21:53 (panda time): restarted
% 2002may13 22:43 (panda time): took another break :)
% 2002may14  0:27 (panda time): restarted
% stopped 2002may14 05:23, but took a 40-min break

(defun hm (h m) (+ (* h 60) m))
(defun hmdif (h1 m1 h2 m2) (- (hm h2 m2) (hm h1 m1)))

(+ (+ (hmdif  16 54  18 30)		; \
      (hmdif  22 54  27 29)		; | first night (may 09 -> may 10)
      -20)				; / 5h51
   (+ (hmdif  18 52  19 47)		; \
      (hmdif  21 53  22 43)		; | second night (may 13 -> may 14)
      (hmdif   0 27   5 23)		; | 6h01
      -40)				; /
   )

% (fooi-re "\\(\\\\myeqref{[0-9.]*}\\){[0-9.]*}" "\\1")
% (find-lsrcfile "base/ltmath.dtx" "\\begin{environment}{equation}")
% (fooi-re "(\\([0-9]\\.[0-9]+\\))" "\\myeqref{\\1}{\\1}")
% (fooi-re "\\(\\\\myeqref{[0-9.]*}{[0-9.]*}\\)" "(\\1)")
% (fooi-re "\\$\\$\\([^$]*\\)\\$\\$" "\\\\begin{equation} \\\\label{??}\\1\\\\end{equation}")
% (find-lsrcfile "base/classes.dtx" "\\begin{environment}{equation}")


# (find-fline "~/LATEX/wolfson.tex" "``serious'' version" "eev")
#
cd ~/LATEX/
dvips -o zawdraft.ps zawdraft.dvi
dvips -o zaw.ps zaw.dvi
gzip -c6 < zaw.ps > zaw.ps.gz
laf wolfson* zaw*.*

#
cd ~/LATEX/
zip zawtexdvi.zip wolfson.tex $(dvi+eps zaw.dvi)
laf wolfson* zaw*.*

#
cd ~/LATEX/
laf wolfson* zaw*.*

#
cd ~/LATEX/
latex wolfson.tex
rexdvi wolfson.dvi
dvips -o wolfson.ps wolfson.dvi

#
cd ~/LATEX/
gv wolfson.ps &

#
cd ~/LATEX/
#pdflatex wolfson.tex		;# doesn't support ".eps"s!
laf      wolfson.pdf
#gv       wolfson.pdf

#
% (find-fline "/debian/usr/share/texmf/source/latex/guide/essential.tex" "|p{1in}|")
#
xdvi /usr/share/texmf/doc/latex/general/essential.dvi &
# (find-node "(latex)eqnarray")
#
cd ~/LATEX/
latex wolfson.tex
pdflatex wolfson.tex
#
cd ~/LATEX/
# anonftp angg.twu.net 'cd public_html/tmp' \
anonftp angg.twu.net 'cd slow_html/tmp' \
  'put zaw.ps.gz' \
  'put wolfson.tex' \
  'put zawdraft.dvi' \
  'put zawdraft.pdf' \
  'put zawdraft.ps' \
  'put zawtexdvi.zip' \
  dir quit
#
echo ... > /tmp/index.html
cd /tmp/
anonftp angg.twu.net 'cd public_html/tmp' \
  'put index.html' \
  dir quit
#
cd ~/LATEX/zawfigs/
for j in 1 2 3 4 5 6;       do cp -v figure5.$j.eps figure5$j.eps; done
for j in 1 2 3 4 5 6 7 8 9; do cp -v figure6.$j.eps figure6$j.eps; done
#
