Reconstruct Kaplan–Meier Estimator as M-estimator and Its Confidence Band

We assume that the survival times to an event of interest are denoted by $T_{1},\ldots,T_{n}$, which are independently and identically distributed (i.i.d.) under a cumulative distribution function $F_{0}$ and the corresponding survival function $S_{0}=1-F_{0}$. In a similar way, we assume i.i.d. censoring times $C_{1},\ldots,C_{n}$ from a censoring distribution $G_{0}$. The observed time of subject $i$ is $X_{i}=\min\{T_{i},C_{i}\}$ with an indicator $\Delta_{i}=I\{T_{i}<C_{i}\}$ which equals 1 if the event of interest is observed before censoring and 0 otherwise. Often, independence is assumed between event time $T_{i}$ and censoring time $C_{i}$ for $i=1,\ldots,n$. Let $X_{(1)}<\cdots<X_{(K)}$ be the $K$ distinct observed event times. In what follows, we define the M-estimator of the survival function and express the Kaplan–Meier estimator as a special case of the M-estimator.

We start with the case where there is no censoring (i.e., $\Delta_{i}=1$ for all $i$). Consider a known functional $m_{S}:\mathcal{S}\to\mathbb{R}$ where $\mathcal{S}=\{S(x):[0,\infty)\to[0,1];S(x)\text{ is nonincreasing}\}$. A popular method to find the estimator $\hat{S}(x)$ is to maximize a criterion function as follows,

$$\hat{S}(x)=\arg\max_{S(x)\in\mathcal{S}}M_{n}({S})=\arg\max_{S(x)\in\mathcal{S}}\frac{1}{n}\sum_{i=1}^{n}m_{S}(X_{i}).$$

One special case of the M-function is the $L_{2}$ functional norm (or a quadratic norm) such that

$$m_{S}(X)=\int_{0}^{\infty}\big{[}-I\{X>x\}^{2}+2S(x)I\{X>x\}-S(x)^{2}\big{]}d\mu(x),$$

where $\mu(x)$ is a cumulative probability function.

Let $\#\{i:\text{ Condition }\}$ be the number of observations that meet the condition. It is clear that when the $L_{2}$ functional norm is used, the empirical M-function is

	$\displaystyle M_{n}(S)$	$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\int_{0}^{\infty}\big{[}-I\{X_{i}>x\}^{2}+2S(x)I\{X_{i}>x\}-S(x)^{2}\big{]}d\mu(x)$		(1)
		$\displaystyle=\int_{0}^{\infty}\Big{[}-\frac{\#\{i:X_{i}>x\}}{n}+2S(x)\frac{\#\{i:X_{i}>x\}}{n}-S(x)^{2}\Big{]}d\mu(x),$

and the Kaplan–Meier estimator

	$\displaystyle\hat{S}(x)$	$\displaystyle=\prod_{k:X_{(k)}\leq x}\Bigg{(}1-\frac{\#\{i:X_{i}=X_{(k)};\Delta_{i}=1\}}{\#\{i:X_{i}\geq X_{(k)}\}}\Bigg{)}$
		$\displaystyle=\prod_{k:X_{(k)}\leq x}\Bigg{(}1-\frac{\#\{i:X_{i}=X_{(k)}\}}{\#\{i:X_{i}\geq X_{(k)}\}}\Bigg{)}$
		$\displaystyle=\prod_{k:X_{(k)}\leq x}\Bigg{(}\frac{\#\{i:X_{i}>X_{(k)}\}}{\#\{i:X_{i}\geq X_{(k)}\}}\Bigg{)}$
		$\displaystyle=\prod_{k:X_{(k)}\leq x}\Bigg{(}\frac{\#\{i:X_{i}\geq X_{(k+1)}\}}{\#\{i:X_{i}\geq X_{(k)}\}}\Bigg{)}$
		$\displaystyle=\frac{\#\{i:X_{i}>x\}}{n}$

is the maximizer of $M_{n}(S)$ in (1).

When there are censored observations in the data, the empirical M-function of the observed data is

$\displaystyle\widetilde{M}_{n}(S)$

$\displaystyle=\frac{1}{n}\sum_{i=1}^{n}\widetilde{m}_{S}(X_{i},\Delta_{i}),$

(2)

where

$$\widetilde{m}_{S}(X,\Delta)=\begin{cases}m_{S}(X),&\Delta=1,\\ \displaystyle\int_{0}^{X}\big{[}-I\{X>x\}^{2}+2S(x)I\{X>x\}-S(x)^{2}\big{]}d\mu(x),&\Delta=0.\\ \end{cases}$$

To obtain the optimizer,

$$\hat{S}(x)=\arg\max_{S(x)\in\mathcal{S}}\widetilde{M}_{n}({S}),$$

(3)

we can apply the EM algorithm as follows:

• •

  E-step: Given the $g$th step estimator $\hat{S}^{(g)}(x)$, compute the expectation of the empirical M-function,

  	$\displaystyle E[M_{n}(S)|\hat{S}^{(g)}]=$	$\displaystyle\frac{1}{n}\sum_{\Delta_{i}=1}\int_{0}^{\infty}\big{[}-I\{X_{i}>x\}^{2}+2S(x)I\{X_{i}>x\}-S(x)^{2}\big{]}d\mu(x)$		(4)
  		$\displaystyle+\frac{1}{n}\sum_{\Delta_{i}=0}\int_{0}^{\infty}\big{[}-\frac{\hat{S}^{(g)}(\max\{x,X_{i}\})}{\hat{S}^{(g)}(X_{i})}+2S(x)\frac{\hat{S}^{(g)}(\max\{x,X_{i}\})}{\hat{S}^{(g)}(X_{i})}-S(x)^{2}\big{]}d\mu(x).$

• •

  M-step: Compute

  $$\hat{S}^{(g+1)}(x)=\arg\max_{S(x)\in\mathcal{S}}E[M_{n}(S)|\hat{S}^{(g)}(x)].$$

  (5)

The validity of this EM algorithm is guaranteed by Theorem 1.

Based on Theorem 1, we conclude that

	$\displaystyle\widetilde{M}_{n}({\hat{S}^{(g+1)}})$	$\displaystyle=E[M_{n}({\hat{S}^{(g+1)}})|\hat{S}^{(g)}]-\Big{(}E[M_{n}(\hat{S}^{(g+1)})|\hat{S}^{(g)}]-\widetilde{M}_{n}({\hat{S}^{(g+1)}})\Big{)}$
		$\displaystyle\geq E[M_{n}({\hat{S}^{(g)}})|\hat{S}^{(g)}]-\Big{(}E[M_{n}(\hat{S}^{(g+1)})|\hat{S}^{(g)}]-\widetilde{M}_{n}({\hat{S}^{(g+1)}})\Big{)}$
		$\displaystyle\geq E[M_{n}({\hat{S}^{(g)}})|\hat{S}^{(g)}]-\Big{(}E[M_{n}(\hat{S}^{(g)})|\hat{S}^{(g)}]-\widetilde{M}_{n}({\hat{S}^{(g+1)}})\Big{)}$
		$\displaystyle=\widetilde{M}_{n}({\hat{S}^{(g)}})$

and thus the M-estimator would be obtained as the convergent point of the EM algorithm. Through a proof by induction, it can be shown that the EM algorithm converges to the KM estimator, implying that the KM estimator is an M-estimator (3).

The proofs of Theorems 1 and 2 are provided in the Appendix.

Define

$$\kappa_{x}(X,\Delta)=\begin{cases}I\{X>x\},&\Delta=1,\\ \displaystyle\frac{S_{0}(\max\{x,X\})}{S_{0}(X)},&\Delta=0,\end{cases}$$

and we have

$$E[\widetilde{m}_{S}(X,\Delta)|S_{0}]=\int_{0}^{\infty}\big{[}-\kappa_{x}(X,\Delta)+2S(x)\kappa_{x}(X,\Delta)-S(x)^{2}\big{]}d\mu(x)$$

and

$$\frac{\partial E[\widetilde{m}_{S}(X,\Delta)|S_{0}]}{\partial S(x)}=2\kappa_{x}(X,\Delta)-2S(x).$$

Given that $X=\min\{T,C\}$ and indicator $\Delta=I\{T<C\}$ where $T\sim F_{0}=1-S_{0}$ and $C\sim G_{0}$, it is straightforward that $E\big{\{}\kappa_{x}(X,\Delta)\big{\}}=S_{0}(x)$. As a result,

$\displaystyle E\Bigg{\{}\frac{\partial E[\widetilde{m}_{S}(X,\Delta)|S_{0}]}{\partial S(x)}\Bigg{\}}=2S_{0}(x)-2S(x),$

which equals 0 when $S(x)=S_{0}(x)$. By the law of large numbers,

$$\frac{\partial E[\widetilde{M}_{n}(S)|S_{0}]}{\partial S(x)}=\frac{1}{n}\sum_{i=1}^{n}\frac{\partial E[\widetilde{m}_{S}(X_{i},\Delta_{i})|S_{0}]}{\partial S(x)}$$

converges to $2S_{0}(x)-2S(x)$ in probability and thus $\hat{S}(x)$ converges to $S_{0}(x)$ in probability.

If there exists an $X_{F}$ where $F_{0}(X_{F})<1$, we have that when $0<x_{1}\leq x_{2}<X_{F}$,

$\displaystyle E\big{\{}\kappa_{x_{1}}(X,\Delta)\kappa_{x_{2}}(X,\Delta)\big{\}}$

$\displaystyle=S_{0}(x_{2})(1-G_{0}(x_{1}))+S_{0}(x_{1})S_{0}(x_{2})\int_{0}^{x_{1}}\frac{dG_{0}(u)}{S_{0}(u)}$

by the conditional expectation formula. As a result, when $S=S_{0}$,

		$\displaystyle E\Bigg{\{}\frac{\partial E[\widetilde{m}_{S}(X,\Delta)|S_{0}]}{\partial S(x_{1})}\frac{\partial E[\widetilde{m}_{S}(X,\Delta)|S_{0}]}{\partial S(x_{2})}\Bigg{\}}$
	$\displaystyle=$	$\displaystyle 4E\{\big{(}\kappa_{x_{1}}(X,\Delta)-S(x_{1})\big{)}\big{(}\kappa_{x_{2}}(X,\Delta)-S(x_{2})\big{)}\}$
	$\displaystyle=$	$\displaystyle 4E\{\kappa_{x_{1}}(X,\Delta)\kappa_{x_{2}}(X,\Delta)\}-4S(x_{1})S(x_{2})$
	$\displaystyle=$	$\displaystyle 4S_{0}(x_{1})S_{0}(x_{2})\Biggl{\{}\frac{1-G_{0}(x_{1})}{S_{0}(x_{1})}+\int_{0}^{x_{1}}\frac{dG_{0}(u)}{S_{0}(u)}-1\Biggr{\}}$
	$\displaystyle=$	$\displaystyle 4S_{0}(x_{1})S_{0}(x_{2})\int_{0}^{x_{1}}\frac{1}{S_{0}^{2}(1-G_{0})}d(1-S_{0}),$

which leads to the joint asymptotic distribution as

$$\sqrt{n}\begin{pmatrix}\hat{S}(x_{1})-S_{0}(x_{1})\\ \hat{S}(x_{2})-S_{0}(x_{2})\\ \end{pmatrix}\mathop{\to}^{\mathcal{D}}N\Biggl{\{}\begin{pmatrix}0\\ 0\end{pmatrix},\begin{pmatrix}\displaystyle S_{0}(x_{1})^{2}\int_{0}^{x_{1}}\frac{d(1-S_{0})}{S_{0}^{2}(1-G_{0})}&\displaystyle S_{0}(x_{1})S_{0}(x_{2})\int_{0}^{x_{1}}\frac{d(1-S_{0})}{S_{0}^{2}(1-G_{0})}\\ \displaystyle S_{0}(x_{1})S_{0}(x_{2})\int_{0}^{x_{1}}\frac{d(1-S_{0})}{S_{0}^{2}(1-G_{0})}&\displaystyle S_{0}(x_{2})^{2}\int_{0}^{x_{2}}\frac{d(1-S_{0})}{S_{0}^{2}(1-G_{0})}\end{pmatrix}\Biggr{\}}.$$

(7)

Under the log-transformation, the pointwise asymptotic distribution of $\log\hat{S}(x)$ is

$$\sqrt{n}(\log\hat{S}(x)-\log S_{0}(x))\mathop{\to}^{\mathcal{D}}N\Big{\{}0,\int_{0}^{x}\frac{1}{S_{0}^{2}(1-G_{0})}d(1-S_{0})\Big{\}},$$

where the variance can be estimated by

$$\sum_{k:X_{(k)}\leq x}\frac{n\Delta_{(k)}}{\big{(}\sum_{l=k}^{K}n_{(l)}\big{)}\big{(}\sum_{l=k+1}^{K}n_{(l)}\big{)}}.$$

Given the asymptotic distribution (7), the asymptotic distribution of the process $\hat{Z}(x)=\sqrt{n}\{\hat{S}(x)-S_{0}(x)\}$ ($0<x<X_{F}$) is provided in Theorem 3.

To deduce the confidence band by Theorem 3, let $A(x)=\int_{0}^{x}{S_{0}(u)^{-2}(1-G_{0}(u))^{-1}}d(1-S_{0}(u))$, $H(x)=A(x)/\{1+A(x)\}$. It is clear that

$$\lim\limits_{x\to\infty}A(x)\geq\lim\limits_{x\to\infty}\int_{0}^{x}S_{0}(u)^{-2}d(1-S_{0}(u))=\infty$$

and $\lim\limits_{x\to\infty}H(x)=1$. For convenience, let $H(x)=1$ for $x\geq X_{F}$. Thus, the covariance function of the zero-mean Gaussian process ${Z}(x)$ in Theorem 3 can be rewritten as

$${\rm cov}\big{\{}{Z}(x_{1}),{Z}(x_{2})\big{\}}=S_{0}(x_{1})S_{0}(x_{2})\frac{H(x_{1})}{1-{H}(x_{1})}=S_{0}(x_{1})S_{0}(x_{2})\frac{H(x_{1})(1-{H}(x_{2}))}{(1-{H}(x_{1}))(1-{H}(x_{2}))}.$$

Under the log-transformation, the process $\hat{Z}^{*}(x)=\sqrt{n}\{\log\hat{S}(x)-\log S_{0}(x)\}$ converges weakly to a zero-mean Gaussian process

$$\frac{B^{0}(H(x))}{1-{H}(x)},\quad 0<x<X_{F},$$

where $B^{0}(\cdot)$ is a standard Brownian bridge process on $[0,1]$. With the constant

$$c_{\alpha}(a,b)=\inf\Biggl{\{}c:\mathbb{P}\Bigg{(}\sup_{[a,b]}|B^{0}(\cdot)|\leq c\Bigg{)}\geq 1-\alpha\Biggr{\}},\quad 0<\alpha<1,0<a<b<1,$$

the asymptotic $(1-\alpha)$ confidence band of the survival function in the interval $[x_{1},x_{2}]\subset[0,X_{F}]$ is

$$\Bigg{[}\hat{S}(x)\exp\Bigg{(}-\frac{c_{\alpha}(\hat{H}(x_{1}),\hat{H}(x_{2}))}{\sqrt{n}(1-\hat{H}(x))}\Bigg{)},\hat{S}(x)\exp\Bigg{(}\frac{c_{\alpha}(\hat{H}(x_{1}),\hat{H}(x_{2}))}{\sqrt{n}(1-\hat{H}(x))}\Bigg{)}\Bigg{]},\quad x_{1}<x<x_{2},$$

where

$$\hat{H}(x)=1-\Biggl{\{}1+\sum_{k:X_{(k)}\leq x}\frac{n\Delta_{(k)}}{\big{(}\sum_{l=k}^{K}n_{(l)}\big{)}\big{(}\sum_{l=k+1}^{K}n_{(l)}\big{)}}\Biggr{\}}^{-1}.$$

We conduct several simulation studies to compare the Kaplan–Meier estimator obtained by optimizing the M-function (2) with the existing maximum likelihood approach (Kaplan and Meier,, 1958). We refer to the Kaplan–Meier estimator obtained by the proposed method as “M-KM” and the existing maximum likelihood approach as “KM”. Both the estimation and inference, including the coverage probability of confidence intervals and confidence bands, are studied. We consider two examples with $100$ replications for each.

Example 1. For each sample of Example 1, $n=200$ observations of survival data are generated, where $F_{0}$ and $G_{0}$ are exponential distributions with rates 1/3 and 1/6 respectively.

Example 2. For each sample of Example 2, $n=500$ observations of survival data are generated, where $F_{0}$ and $G_{0}$ are Weibull$(1,1)$ and Weibull$(1,2)$ respectively.

Figure 1 displays the KM estimators, pointwise 95% confidence intervals and simultaneous 95% confidence bands computed by the proposed method and existing maximum likelihood approach under one sample from Examples 1 and 2. It can be seen that although with different target functions, the KM estimator by the proposed method is the same as the traditional one. The 95% confidence intervals are exactly the same as that by Greenwood’s formula. In addition, the 95% confidence band coincides with the one proposed by Hall and Wellner, (1980).

Tables 1 and 2 show the coverage probability and the average length of the 95% confidence interval at several time points over 100 repetitions for Examples 1 and 2. We choose seven time points $1,2,\ldots,7$ in Example 1 and eight time points $0.1,0.3,\ldots,1.5$ in Example 2. Given that the KM estimator obtained by the proposed method equals the one computed via the maximum likelihood approach, the results of M-KM is the same as the existing one. The coverage probability of the 95% confidence interval in Example 1 is around 0.95 at all seven time points. However, in Example 2, the coverage probability is less than 0.95 for larger time points. We also compute the coverage probability for the confidence bands. In Example 1, the coverage probability of the confidence band is 0.97, while the value is 0.94 in Example 2. This means that the properties of the KM estimator can be successively recovered with the M-estimation theory, implying that KM estimator can be interpreted as an M-estimator.

We apply the proposed method to two datasets to evaluate its performance. The first dataset is from a diabetic study, including $n=394$ observations and 8 variables. Among all observations, 197 patients are labeled “high-risk” diabetic retinopathy according to the diabetic retinopathy study (DRS). For each patient, one eye was randomly chosen to receive the laser treatment and the other one received no treatment. The event of interest is that the visual acuity dropped below 5/200. Some records were censored due to death, dropout, or the end of the study. The second dataset is from a lung cancer study involving the North Central Cancer Treatment Group of 228 patients of advance lung cancer. The survival time was measured from initiation of the study to the time when the patient died or was censored.

The KM curves and the corresponding 95% confidence intervals and 95% confidence band computed by the proposed method and existing maximum likelihood approach for both datasets are shown in Figure 2. Again, the estimation and inference result of the proposed method are the same as the existing ones, suggesting that the proposed method provides a new interpretation of the KM estimator as an M-estimator.