Large Deviations of the Estimated Cumulative Hazard Rate

TL;DR

This paper investigates the large deviation properties of the Nelson-Aalen estimator for the cumulative hazard rate, revealing asymmetries that cause overestimation especially in small samples and rare event scenarios.

Contribution

It derives large deviation principles for the Nelson-Aalen estimator, providing new insights into its small sample behavior and bias tendencies.

Findings

Rate functions are asymmetric, favoring overestimation.

Overestimation bias is strongest with small samples.

Implications for risk assessment of rare events.

Abstract

Survivorship analysis allows to statistically analyze situations that can be modeled as waiting times to an event. These waiting times are characterized by the cumulative hazard rate, which can be estimated by the Nelson-Aalen estimator or diverse confidence estimators based on asymptotic statistics. To better understand the small sample properties of these estimators, the speed of convergence of the estimate to the exact value is examined. This is done by deriving large deviation principles and their rate functions for the estimators and examining their properties. It is shown that these rate functions are asymmetric, leading to a tendency of the estimated cumulative hazard rate to overestimate the true cumulative hazard rate. This tendency is strongest in the cases of (1) small sample sizes and (2) low tail probabilities. Taking this tendency into account can improve risk assessments of rare events and of cases where only little data is available.

Equations

$$H(t):=-\ln (S(t)).$$

$$H_n(t):=-\ln (S_n(t)),$$

$$P\left(\frac{1}{n}\sum _{i=1}^n{X}_i\in C\right)\approx \exp (-n\inf I(C))$$

$$Y_i=\{\begin{array}{ll}1& \text{ if }{X}_i>\tilde{t}\\ 0& \text{ if }{X}_i\le \tilde{t}\end{array},$$

$$H_n(\tilde{t})=-\ln \left(\frac{1}{n}\sum _{i=1}^n{Y}_i\right).$$

$$(P_n{)}_{n\in \mathbb{N}}={\left(\mathcal{L}\left(\frac{1}{n}\sum _{i=1}^n{Y}_i\right)\right)}_{n\in \mathbb{N}}$$

$${\tilde{I}}_p(x)=x\ln \left(\frac{x}{p}\right)+(1-x)\ln \left(\frac{1-x}{1-p}\right)$$

$$I_p(y)$$

$$=\exp (-y)\ln \left(\frac{\exp (-y)}{p}\right)+(1-\exp (-y))\ln \left(\frac{1-\exp (-y)}{1-p}\right)$$

$$J_p(z)$$

$$=-zp\exp (-z)+(1-p\exp (-z))\ln \left(\frac{1-p\exp (-z)}{1-p}\right)$$

$$\frac{d}{dp}{J}_p(z)=-z\exp (-z)+\frac{1-\exp (-z)}{1-p}-\exp (-z)\ln \left(\frac{1-p\exp (-z)}{1-p}\right)$$

$$\frac{d^2}{d{p}^2}{J}_p(z)=\frac{\exp (-z)(\exp (z)-1{)}^2}{(1-p{)}^2(\exp (z)-p)}.$$

$$\frac{d}{dp}{J}_0(z)=\exp (-z)(-z-1)+1>0\text{ for }z\ne 0.$$

$$(1-x)\ln (1-x)=-x+\sum _{k=2}^{\infty }\frac{x^k}{k^2-k}\text{ for }|x|<1$$

$$J_p(z)=p\exp (-z)\left[-z+\ln \left(\frac{1-p}{e}\right)+\sum _{k=2}^{\infty }\frac{(p\exp (-z){)}^{k-1}}{k^2-k}\right]-\ln (1-p)$$

$$2\sinh (z)=\exp (z)-\exp (-z)\text{ and }2\cosh (z)=\exp (z)+\exp (-z)$$

$$\displaystyle\big{|}J_{p}(-z)-J_{p}(z)\big{|}$$

$$\displaystyle\quad\quad+\sum_{k=2}^{\infty}\frac{2p^{k}\sinh(kz)}{k^{2}-k}\;\Big{|}\;.$$

$$\big{|}J_{p}(-z)-J_{p}(z)\big{|}\approx 2p\left(z\cosh(z)+\sinh(z)\ln(\tfrac{1-p}{e})\right)+p^{2}\sinh(2z)$$

$$I(y,t)$$

$$=\exp (-y)\ln \left(\frac{\exp (-y)}{S(t)}\right)+(1-\exp (-y))\ln \left(\frac{1-\exp (-y)}{F(t)}\right).$$

Full text

Keywords: Survival Analysis; Large Deviations; Nelson-Aalen estimator

1 Introduction

Survival analysis is the standard framework to statistically analyze waiting times until an event occurs (Miličič, 2008). Classical applications of survival analysis arise in epidemiology, where the event can be the recovery or the death of a patient, and the results of the statistical analysis can decide on the admission or nonadmission of a new medical procedure (Sasieni and Brentnall, 2014; Kantoff et al., 2010).

A fundamental concept in characterizing the waiting time up to an event is the cumulative hazard rate (short CH)(Aalen, Borgan, and Gjessing, 2008). It is estimated using the Nelson-Aalen estimator and can be complemented by different confidence intervals and bands to compensate for uncertainties of the results (Aalen, Borgan, and Gjessing, 2008; Aalen, 1978; Bie, Borgan, and Liestøl, 1986; Nelson, 1969, 1972).

However these confidence area estimators are based on the asymptotic properties of the Nelson-Aalen estimator(Bie, Borgan, and Liestøl, 1986), making it impossible to derive analytical statements regarding their performance for small sample sizes. As a result, no reliable assessments of the uncertainties for small samples are available, which is problematic for studies where larger sample sizes cannot be achieved, be it because of ethical considerations, financial reasons, or rarity.

The aim of this paper is to derive analytical expressions regarding the small sample properties of the estimated CH to better assess its deviations from the true CH in these cases.

For this, the pointwise speed of convergence of the estimated CH to the true CH is examined using the theory of large deviations.

2 Outline

The conventions $\ln(0):=\lim_{x\downarrow 0}\ln(x)=-\infty$ and $0\ln(0):=0$ will be used throughout the paper.

Let $X,X_{1},X_{2},\dots$ be i.i.d. positive random variables, whose values model the waiting times for an event to occur. With the survival function $S(t):=P(X>t)$(Kleinbaum and Klein, 2010, p. 9), the cumulative hazard rate (short CH) can be written as(Kleinbaum and Klein, 2010, p. 294)

[TABLE]

It can be estimated using the $n$-th empirical CH

[TABLE]

where $S_{n}(t)=\tfrac{1}{n}\sum_{i=1}^{n}\mathbf{1}_{(t,+\infty)}(X_{i})$ is the $n$-th empirical survival function. The aim of this paper is to examine the pointwise speed of convergence of $H_{n}$ towards $H$.

From the theory of large deviations, it is known that for integrable i.i.d. $(X_{i})_{i\in\mathbb{N}}$, the relation

[TABLE]

holds for a large class of sets $C$(Varadhan, 2016). The function $I$ is called a rate function and determines the speed of convergence of the averaged random variables towards their expectation value.

In this paper, a rate function for the $H_{n}$ is derived (section 3). The properties of this rate function are then examined to draw conclusions about the convergence of $H_{n}$ to $H$ (section (4).

For this, fix any $\tilde{t}\in\mathbb{R}^{+}$, and define the tail probabilities as $p:=p(\tilde{t})=P(X>\tilde{t})$. Assume that $p\in(0,1)$, since the cases $p=0$ and $p=1$ are trivial. Define the i.i.d. random variables

[TABLE]

so $Y_{i}\sim\operatorname{Ber}(p)$. The behavior of $H_{n}$ at $\tilde{t}$ is uniquely determined by the $Y_{i}$, since

[TABLE]

3 Establishing the Rate Functions

First, a large deviation principle (LDP) for the $Y_{i}$ is derived. By either Cramér’s theorem (Cramér, 1938)(Klenke, 2008, p. 508) or by Sanov’s theorem (Sanov, 1958)(Klenke, 2008, p. 518), the series of probability measures

[TABLE]

satisfies a large deviation principle with rate $n$ and rate function

[TABLE]

for $x\in[0,1]$. This is the relative entropy of two Bernoulli distributions, one with success probability $p$ and one with success probability $x$(Klenke, 2008, p. 515).

Applying the contraction principle (Klenke, 2008, p. 518) to this LDP and the function $f(x)=-\ln(x)$ shows that the series $\left(\mathcal{L}\left(-\ln\left(\tfrac{1}{n}\sum_{k=1}^{n}Y_{k}\right)\right)\right)_{n\in\mathbb{N}}$ satisfies a LDP with rate $n$ and rate function

[TABLE]

for $y\in[0,+\infty]$. This is the rate function of $H_{n}(\tilde{t})$, and is displayed in fig. 1 for four different values of $p$. Next, substitute $y=z-\ln(p)$ and define the centered rate function as

[TABLE]

for $z\in[\ln(p),+\infty]$ and $p\in(0,1)$. It is defined since the main interest of this examination is to compare the behavior of the rate functions $I_{p}$ for different $p$ close to the pointwise limit value $-\ln(p)$ of the empirical cumulative hazard rate. This value is shifted to the origin in the centered rate function and therefore allows to compare the behavior of the empirical cumulative hazard rate at fixed distances from the respective limit values for different $p$.

4 Properties of the Rate Functions

4.1 Monotonicity in p

In this section, it is shown that $J_{p}$, taken as a function of $p$, is strictly increasing. This implies that the speed of convergence is decreasing as $\tilde{t}$ increases and correspondingly the tail probabilities decrease. Without loss of generality, it is assumed that $X$ takes on every tail probability, so $J_{p}$ is well defined for all $p\in(0,1)$.

First, the function $J_{p}$ can not be expected to be strictly increasing in $p$ for $z=0$, since by definition $J_{p}(0)=0$ for all $p\in(0,1)$. Therefore the case $z=0$ will be excluded.

The first derivative of $J_{p}$ with respect to $p$ is given by

[TABLE]

and the second derivative of $J_{p}$ with respect to $p$ by

[TABLE]

Since the inequality $\exp(z)-p>0$ holds for all $p\in(0,1)$ and all $z$ in the domain of $J_{p}$, termwise analysis of eq. (13) shows that the second derivative is positive for all feasible $p$, and zero only when $z=0$, which was excluded above. This makes the second derivative strictly positive, so $J_{p}$ is strictly convex in $p$.

If $J_{p}$ as defined in eq. (11) is taken as a function of $p$ for $p<1$, it is well-defined and its first derivative, evaluated at $p=0$, yields

[TABLE]

So $J_{p}$ is strictly convex in $p$ and its gradient at $p=0$ is strictly positive, therefore $J_{p}$ is strictly increasing in $p$ for $p\in(0,1)$. This shows that the rate of convergence $J_{p}(z)$ is decreasing for all $z$ as the tail probabilities $p$ decrease.

4.2 Asymmetry in z

First, it is shown that $J_{p}$ is not axis symmetric with respect to the ordinate axis. As an aid, the identity

[TABLE]

is used. Splitting the fraction in the logarithm in $J_{p}$ from equation (11) and then applying the power series from eq. (15) with $x=p\exp(-z)$ yields

[TABLE]

for $p\exp(-z)<1$. Therefore $J_{p}$ is asymmetric by the asymmetry of $\exp(-z)$.

Next, the symmetry defect of $J_{p}$, given by $|J_{p}(-z)-J_{p}(z)|$, is examined. For this, the identities

[TABLE]

are used. With the representation of $J_{p}$ in eq. (16), directly subtracting $J_{p}(z)$ and $J_{p}(-z)$ yields

[TABLE]

Using only the second order term of the sum gives the approximation for the symmetry defect

[TABLE]

for $p\in(0,1)$ and $z\in(0,-\ln(p))$. These two statements show that the rate of convergence is not symmetric around the true cumulative hazard rate $H$, and will always be slower from above the cumulative hazard rate than from below.

5 Example

Let $X$ be a positive random variable with distribution function $F$ and survival function $S$. Taking the tail probabilities $p$ as functions of $t$, meaning $p=p(t)=S(t)$, the rate function $I_{p}$ in equation (9) is a function in the variables $t$ and $y$:

[TABLE]

For the case where $X$ is exponentially distributed with mean $\tfrac{1}{\lambda}=2$, the contour plot of $I(y,t)$ is displayed in fig. 2, alongside the cumulative hazard rate. The asymmetry of the rate functions is clearly visible.

6 Discussion

The main result of this paper is that the estimate of the cumulative hazard rate (short CH) has a tendency to overestimate the true CH, which is strongest in the cases of (1) small sample sizes or (2) small tail probabilities. This is a direct conclusion from the asymmetry of the rate function governing the convergence of the estimated cumulative hazard rate in combination with the fundamental relation of the theory of large deviations given in equation (3). The degree of this effect is determined by the tail probabilities and therefore unique for every distribution.

It is notable that the asymmetry observed is not generated by the distribution of the random variables, but is rooted in the asymmetry of the rate function of the Bernoulli distribution. It is amplified by the logarithm in the definition of the cumulative hazard rate, and unavoidable in the sense that it naturally arises from observing whether a random variable takes on a value over or under a given threshold.

What remains unclear is how these results translate into the standard framework of survival analysis, i.e. censoring and using point processes instead of i.i.d. random variables. Especially censoring might have a strong influence on the case with low tail probabilities, since early censoring can render long survival times (that are commonly associated with low tail probabilities) irrelevant. Extending the results presented in this papers to this more general setting is relevant for applications, but requires further work.