Actuarial (lifetable) estimator of the survival

Revision: Kaplan-Meier estimator

The Kaplan-Meier [KM] estimator of the survival function $$\widehat{S}(t) = \prod_{\{i: t_i\leq t\}}\biggl[1-\frac{\Delta\overline{N}(t_i)}{\overline{Y}(t_i)}\biggr] = \widehat{S}(t-) \cdot \biggl[1-\frac{\Delta\overline{N}(t)}{\overline{Y}(t)}\biggr]$$ is motivated by decomposition of the survival probability by chain rule. If the survival distribution is continuous the increment process $\Delta \overline{N}$ attains only values 0 (no death at time $t$) or 1 (exactly one death at time $t$). In practice, however, time is not measured continuously or with necessary precision. Then ties begin to occur in the data and more death at the same time can be observed. If the ties come “accidentally”, for example, by rounding the true value, and its presence is rare, then [KM] works still fine. However, when the grouping of time (now rather categorical than continuous) comes from the design of the study itself, the resulting estimate could possibly result in unsatisfactory outputs.

Several methods for adjusting this fact have been proposed, we will introduce the Lifetable or Actuarial estimator.

Lifetable or actuarial estimator

Applies when the event times are grouped

Our goal is still to estimate the survival function, hazard, and density function, but this is complicated by the fact that we don’t know exactly when during each time interval an event occurs.

Notation:

• the $j$-th time interval is $[t_{j-1},t_j)$,
• $c_j$ - the number of censorings in the $j$-th interval,
• $d_j$ - the number of failures (deaths) in the $j$-th interval,
• $r_j$ - the number of units entering the interval.

We could apply the K-M formula directly. However, it treats the problem as though it were in discrete time, with events happening only at endpoints of intervals. But the true events event times lie somewhere within the preceding interval. In fact, what we are trying to calculate here is the conditional probability of dying within the interval, given survival to the beginning of it.

We now make the assumption that the censoring process is such that the censored survival times occur uniformly throughout the $j$-th interval. So, we can assume that censorings occur on average halfway through the interval: $$r_j'=r_j-c_j/2.$$ I.e., average number of individuals who are at risk during this interval. The assumption yields the Actuarial Estimator. It is appropriate if censorings occur uniformly throughout the interval. This assumption is sometimes known as the actuarial assumption.

The event that an individual survives beyond $[t_{j-1},t_j)$ is equivalent to the event that an individual survives beyond $[t_{0},t_1)$ and … and that an individual survives beyond $[t_{j-2},t_{j-1})$ and that an individual survives beyond $[t_{j-1},t_j)$. Let us define the following quantities:

• $P_j=P\{\mbox{an individual survives beyond }[t_{j-1},t_j)\}$,
• $p_j=P\{\mbox{an individual survives beyond }[t_{j-1},t_j) |\mbox{ he survives beyond }[t_{j-2},t_{j-1})\}, \quad j>1$,
• $p_1=P_1.$

Then by chain rule, $$P_j=p_1 \cdot \; \cdots \; \cdot p_j.$$ The conditional probability of an death in $[t_{j-1},t_j)$ given that the individual survives beyond $[t_{j-2},t_{j-1})$ is estimated by $$d_j/r_j'.$$ Thus, $p_j$ is estimated by $$\frac{r_j'-d_j}{r_j'}.$$ So, the actuarial (lifetime) estimate of the survival function is $$\widehat{S}(t)=\prod_{j=1}^{k-1}\frac{r_j'-d_j}{r_j'} \qquad \text{for } t_{k-1} \leq t < t_{k}.$$

Remarks:

• Because the intervals are defined as $[t_{j-1}, t_j)$, the first interval typically starts with $t_0 = 0$.
• The implication in R is that $\widehat{S}(t_0)=1$.

The form of the estimated survival function obtained using this method is sensitive to the choice of the intervals used in its construction. On the other hand, the lifetable estimate is particularly well suited to situations in which the actual death times are unknown, and only available information is the number of deaths and the number of censored observations which occur in a series of consecutive time intervals.

Another way to deal with such type of data is to use Interval-censoring approaches that better fit this situation. Traditional non-parametric estimator of survival function is then Turnbull estimator (very similar to [KM]) which accounts for jumps not at a given time but during whole intervals. More can be found in Survival Analysis with Interval-Censored Data (2018) by Bogaerts, Komárek and Lesaffre.

library(KMsurv)
lifetab(tis,    #[K+1] a vector of end points of time intervals (including t_0)
        ninit,  #[1]   the number of subjects initially entering the study
        nlost,  #[K]   a vector of the number of individuals lost follow or withdrawn alive for whatever reason
        nevent) #[K]   a vector of the number of individuals who experienced the event

Duration of nursing home stay data

The National Center for Health Services Research studied 36 for-profit nursing homes to assess the effects of different financial incentives on length of stay. ‘Treated’ nursing homes received higher per diems for Medicaid patients, and bonuses for improving a patient’s health and sending them home.

Study included 1601 patients admitted between May 1, 1981 and April 30, 1982.

Variables include:

• los - Length of stay of a resident (in days)
• age - Age of a resident
• rx - Nursing home assignment (1:bonuses, 0:no bonuses)
• gender - Gender (1:male, 0:female)
• married - (1:married, 0:not married)
• health - Health status (2:second best, 5:worst)
• censor - Censoring indicator (1:censored, 0:discharged) (Beware!)

nurshome = read.table("http://www.karlin.mff.cuni.cz/~vavraj/cda/data/NursingHome.dat",
  head=F, skip=14, col.names = c("los", "age", "rx", "gender", "married", "health", "censor"))

In this dataset, length of stay (los) can be considered continuous enough so that classical Kaplan-Meier estimator looks fairly smooth:

We wish to demonstrate the use of the actuarial method. To do so, consider the treated nursing home patients, with length of stay (los) grouped into 100 day intervals:

nurshome$int = floor(nurshome$los/100)*100    # group in 100 day intervals

Comparison of different approaches of estimation:

We can see that for low values of time there is merely any difference since $r_j'$ and $\overline{Y}(t_j)$ are large values with almost no impact on ratios used for estimation. However, there seems to emerge some difference for high values of time, since the negligible differences have accumulated and the actuarial update of ratios begins to notably lower the estimating ratios.

library(muhaz)
haz = muhaz(nurshome$los, nurshome$censor==0)
plot(haz)

Mortality in the Czech Republic in 2019

Download dataset mort.RData. There are three dataframes mort.m, mort.f and mort containing mortality data from the Czech Republic in the year 2019 for males, females and both genders combined. Each dataframe has three columns:

• i is an age group (going from the age $i$ years to $i+1$ years),
• Di is the number of deaths in the $i$-th age group in 2019,
• Yi is the number of people alive when they reached the age $i$ (all at risk in $i$-th age group).

Think properly about the structure of this dataset. Does it follow the same structure as supposed in previous Exercise and the one supposed above? What does belonging to a certain age group mean?

The ratio (Yi-Di)/Yi estimates the probability of survival of another year in $i$-th age group. The resulting estimator is then:

Inspired by the actuarial approach we base the estimate on (Yi/2-Di)/(Yi/2) and obtain the following estimate, which inappropriately shortens the life span: If we blindly use lifetab function we obtain: