Survival Analysis Note

Introduce to Survival Analysis

Survival Function and Hazard Rate

Survival Function

Definition

Let $T\ge 0$ be a non-negative Random Variable representing survival time with PDF $f(t)$ and CDF $F(t)$. The survival function of the random variable is defined as $S(t)=1-F(t)=P(T>t)$

The survival function is a function that gives the probability that a patient, device, or other object of interest will survive past a certain time.

Link to original

Hazard Function

Definition

Hazard Function

Let $T\ge 0$ be a non-negative Random Variable representing survival time with PDF $f(t)$ and CDF $F(t)$. The hazard function $0<\lambda (t)<\infty$ of the random variable $T$ is defined as $\lambda (t)=\frac{f(t)}{S(t)}={\lim }_{\Delta t\to 0}\frac{P(T<t+\Delta t|T>t)}{\Delta t}$ where $S(t)$ is the Survival Function.

The hazard function refers to the rate of occurring event at a given time $t$.

Cumulative Hazard Function

The hazard function can alternatively be represented in terms of the cumulative hazard function, defined as $\Lambda (t)={\int }_0^t\lambda (u)du$

Facts

The Survival Function $S(t)$, the Cumulative Hazard Function, the density (PDF) $f(t)$, the Hazard Function, and the distribution function (CDF) of survival time $F(t)$ are related through $S(t)=\exp [-\Lambda (t)]=\frac{f(t)}{\lambda (t)}=1-F(t)$

Link to original

Types of Censoring

Right Censoring

Kinds

Type 1 Censoring

Definition

Type 1 Censoring

when $c_1=\cdots =c_n=t_c$ ($c_i$’s are constant)

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator. We observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$.

In a type 1 censoring setting, the censored times $C_i$ are fixed constants, not random variables.

The PDF of the observation $Y_i$ is derived as $\begin{array}{rl}P(Y_i=C_i,{\delta }_i=0)& =P(T_i>C_i)=S(y_i)\\ P(Y_i=T_i,{\delta }_i=1)& =f(y_i)\end{array}$

Likelihood of Type 1 Censoring Data

The Likelihood Function of type 1 censoring data is defined as $L(\theta )=\prod _{i=1}^{n}f(y_i{)}^{{\delta }_i}S(y_i{)}^{1-{\delta }_i}=\prod _{i=1}^n\lambda (y_i{)}^{{\delta }_i}S(y_i)=\prod _{i=1}^n\lambda (y_i{)}^{{\delta }_i}\exp [-\Lambda (y_i)]$

Link to original

Type 2 Censoring

Definition

Type 2 Censoring

when $r=3$

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator. We observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$.

In a type 2 censoring setting, we observe first $r$ out of $n$ experiment $T_1,\dots ,T_r$. In other words, for the order statistics $T_{(i)}$ of $T_i$, we only observe $T_{(1)}\le T_{(2)}\le \cdots \le T_{(r)}$. Where $C_1=C_2=\cdots =C_n$ are not constants, but random variables.

Likelihood of Type 2 Censoring Data

The likelihood function of type 2 censored data can be computed using the same equation used for type 1 censored data but computing the joint-PDF of the Order Statistic $T_{(1)}\le \cdots \le T_{(r)}$ is easier.

The joint PDF of $T_{(1)}\le \cdots \le T_{(r)}$ is derived as $L(\theta )=\frac{n!}{(n-r)!}\left[\prod _{i=1}^{r}f(t_{(i)})\right]S(t_{(r)}{)}^{n-r}$

Link to original

Random Censoring

Definition

Random Censoring

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times, which can be both random variable or constant. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator.

In a random censoring setting, we only observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$, and the censoring times are $C_1,\dots ,C_n$ are i.i.d. Random Variable follows PDF $g_{\gamma }$ and CDF $G_{\gamma }$.

The PDF of the observation $Y_i$ is derived as $\begin{array}{rl}P(Y_i=t,{\delta }_i=0)& =P(C_i=t,T_i>C_i)=g_{\gamma }(t)S_{\theta }(t)\\ P(Y_i=t,{\delta }_i=1)& =P(T_i=t,T_i\le C_i)=f_{\theta }(t)(1-G_{\gamma }(t))\end{array}$ where $\theta$ is the parameter of interest, and $\gamma$ is the: nuisance parameter

Likelihood of Random Censoring Data

The Likelihood Function of random censored data is defined as $L(\theta )=\prod _{i=1}^n[f_{\theta }(t)(1-G_{\gamma }(t)){]}^{{\delta }_i}[S_{\theta }(t)g_{\gamma }(t){]}^{1-{\delta }_i}=c\cdot \prod _{i=1}^n{f}_{\theta }(t{)}^{{\delta }_i}{S}_{\theta }(y_i{)}^{1-{\delta }_i}$ where $c$ is a constant.

Link to original

Link to original

Left Censoring

Definition

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times, and $C_1,C_2,\dots ,C_n$ are the censoring times, which can be both random variable or constant. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\ge C_i\text{: uncensored}\\ 0& T_i<C_i\text{: censored}\end{array}$ be the censoring indicator.

In a left censoring setting, we observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\max (T_i,C_i)$. In other words the event has already occurred before it becomes the subject of observed.

Link to original

Interval-Censored Data

Definition

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times. Interval censored data is given as an interval, not an exact point of time. We only observe the interval $(L_i,R_i]$ that includes $T_i$ Interval-censored data is divided into four cases.

Case 1 Interval-Censored Data (Current Status Data)

Data is given the form of $[0,C_i]$ or $(C_i,\infty )$, where $C_i$ is a fixed time point.

Case 2 Interval-Censored Data

Data is given the form of $(L_i,R_i]$ where $0<L_i<R_i<\infty$.

Double Censored Data

Data is given the form of $(L_i,R_i]$ where $0\le L_i\le R_i\le \infty$.

If $R_i=\infty$, then it is right-censored data, if $L_i=0$, then it is left-censored data, and if $L_i=R_i$, then it is uncensored data.

Panel Data

Observations are made at discrete time points. The period between these observations can be viewed as an interval.

Link to original

Mean Imputation Method

Definition

Mean imputation method substitute the given interval of case 2 interval-censored data $(L_i,R_i]$ with the mean of the interval $Y_i=\frac{(L_i+R_i)}{2}$ if the interval is finite. If $R_i=\infty$, the data $(L_i,\infty )$ is substituted with the left value $Y_i:=L_i$ and treated as right-censored data.

Link to original

Parametric Models

Distributions

Survival Models based on Distributions

Kinds

Exponential Distribution

Exponential Distribution assumes that the Hazard Function is a constant regardless of time. The distribution is memoryless $P(T>t+s|T>s)=P(T>t),0<s<t<\infty$.

The PDF of survival time follows Exponential Distribution $f(x)=\lambda e^{-\lambda t}$ The Hazard Function $\lambda (t)=\lambda$ The Survival Function $S(t)=e^{-\lambda t}$

Gamma Distribution

The PDF of survival time follows Gamma Distribution $f(t)=\frac{{\lambda }^{\alpha }}{\Gamma (\alpha )}{t}^{\alpha -1}{e}^{-\lambda t}$ The Hazard Function and Survival Function don’t have closed form expression.

Weibull Distribution

Weibull Distribution uses additional parameter $\alpha$ than Exponential Distribution to control the shape of the Hazard Function. If $\alpha 1$ then the hazard function is increasing, if $\alpha =1$ then it is constant, and if $0<\alpha <1$ then it is decreasing.

The PDF of survival time follows Weibull Distribution $f(t)=\alpha \lambda (\lambda t{)}^{\alpha -1}\exp [-(\lambda t{)}^{\alpha }]$ where $\lambda$ is a scale parameter and $\alpha$ is a shape parameter.

The Hazard Function $\lambda (t)=\alpha \lambda (\lambda t{)}^{\alpha -1}$ The Survival Function $S(t)=\exp [-(\lambda t{)}^{\alpha }]$

Rayleigh Distribution

The PDF of survival time follows Rayleigh distribution $f(x)=({\lambda }_0+{\lambda }_{1}t)\exp (-{\lambda }_{0}t-\frac{1}{2}{\lambda }_1{t}^2)$ The Hazard Function $\lambda (t)={\lambda }_0+{\lambda }_{1}t$ The Survival Function $S(t)=\exp (-{\lambda }_{0}t-\frac{1}{2}{\lambda }_1{t}^2)$

Log-Normal Distribution

We assume that the survival time follows a Log-Normal Distribution $T\sim Log-normal(\mu ,{\sigma }^2)$ or $\log T\sim N(\mu ,{\sigma }^2)$. The Hazard Function of the distribution is hump-shaped.

The PDF of survival time follows Log-Normal Distribution $f(t)=\frac{1}{\sqrt{2\pi \sigma }}\frac{1}{t}\exp \left(-\frac{(\log t-\mu {)}^2}{2{\sigma }^2}\right)$ The Survival Function $S(t)=1-\Phi \left(\frac{\log t-\mu }{\sigma }\right)$ where $\Phi$ is the CDF of standard normal distribution

The Hazard Function doesn’t have closed form expression.

Gompertz Distribution

The Hazard Function $\ln \lambda (t)={\lambda }_0+{\lambda }_{1}t\iff \lambda (t)=\lambda e^{\gamma t}$ where $\lambda :=e^{{\lambda }_0}$ and $\gamma ={\lambda }_1$

Gompertz-Makeham Distribution

The Hazard Function $\lambda (t)=\alpha +\lambda (t)=\lambda e^{\gamma t}$

Link to original

Survival Models based on Log-Lifetime

Definition

Assume that the log survival time $T$ can be modeled with Log-Linear Model $Y=\ln T=\sigma W+\alpha$ where $\sigma$ is the Scale Parameter, $\alpha$ is the Location Parameter, and $W$ is some well-known distribution.

Kinds

Standard Gumbel Distribution

If $W$ follows the standard Gumbel distribution (extreme value distribution), then $f(y)=\frac{1}{\sigma }\exp \left[(\frac{y-\alpha }{\sigma })-\exp (\frac{y-\alpha }{\sigma })\right]$ $f(t)=\frac{1}{\sigma t}\exp \left[(\frac{\ln t-\alpha }{\sigma })-\exp (\frac{\ln t-\alpha }{\sigma })\right]$ If $\alpha =-\ln \lambda$ and $\sigma =\frac{1}{\rho }$, then $f(t)=\rho \lambda (\lambda t{)}^{\rho -1}\exp [-(\lambda t{)}^{\rho }]$ follows a Weibull Distribution.

Normal Distribution

If $W$ follows the standard normal distribution, then $T$ follows log-normal distribution $f(t)=\frac{1}{\sqrt{2\pi }}\frac{\rho }{t}\exp \left[-\frac{{\rho }^2}{2}(\ln \lambda t{)}^2\right]$ where $\alpha =-\log \lambda$ and $\sigma =\frac{1}{\rho }$

Logistic Distribution

If $W$ follows a Logistic Distribution, then $T$ follows Log-Logistic Distribution $f(y)=\frac{1}{\sigma }\frac{\exp (\frac{y-\alpha }{\sigma })}{[1+\exp (\frac{y-\alpha }{\sigma }){]}^2}$ $f(t)=\lambda \rho \frac{(\lambda t{)}^{\rho -1}}{[1+(\lambda t{)}^{\rho }{]}^2}$ where $\alpha =-\ln \lambda$ and $\sigma =\frac{1}{\rho }$.

Gumbel Distribution

If $W$ follows a Gumbel Distribution (generalized extreme value distribution), then $f(t)=\lambda \rho (\lambda t{)}^{\rho k-1}\frac{\exp [-(\lambda t{)}^{\rho }]}{\Gamma (k)}$ where $\alpha =-\ln \lambda$ and $\sigma =\frac{1}{\rho }$.

Special Cases

If $\rho =1$ then $T\sim \Gamma (k,\lambda )$ If $k=1$ then $T\sim \mathrm{Weibull}(\lambda ,k)$ If $\rho =k=1$ then $T\sim \mathrm{Exp}(\lambda )$ If $k\to \infty$ then $T\sim \mathrm{Lognormal}(\mu ,{\sigma }^2)$

Exponential-F Distribution

If $e^W$ follows a F-Distribution $F(2m,2n)$, then $T$ follows a generalized F-distribution.

Special Cases

If $m=n=1$ then $T$ follows Log-Logistic Distribution If $n\to \infty$ then $T$ follows generalized gamma distribution If $m,n\to \infty$ then $T$ follows Log-Normal Distribution

Link to original

Survival Models with Surviving Fractions

Definition

Consider an event such as death from a specific cause on the incidence of a particular disease.
Define a binary variable $z$ where $z=1$ indicates that an individual will experience the event eventually, and $z=0$ indicates that the individual will never experience the event.
Let $T$ denote the time of occurrence of the event with PDF $f_{\theta }$ and CDF $F_{\theta }$, $T$ is defined only when $z=1$, and $T$ and $z$ are independent. We observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i{C}_i\text{: censored}\end{array}$.

Suppose that the probability of $z$ is $p$, $p=P(z=1)$.

No Covariate Case

$\begin{array}{rl}P({\delta }_i=1,T_i=Y_i,z_{i=1})& =p{f}_{\theta }(y_i)\\ P({\delta }_i=0,T_i{C}_i,z_{i=1})& =(1-p)+p{S}_{\theta }(y_i)\end{array}$

The likelihood is defined as $L(\theta )=\prod _{i=1}^n[p{f}_{\theta }(y_i){]}^{{\delta }_i}[(1-p)+p{S}_{\theta }(y_i){]}^{1-{\delta }_i}$

With Covariates Case

The likelihood is defined as $L(\theta )=\prod _{i=1}^n[p({\mathbf{x}}_i)f_{\theta }(y_i|{\mathbf{x}}_i){]}^{{\delta }_i}[(1-p({\mathbf{x}}_i))+p({\mathbf{x}}_i)S_{\theta }(y_i|{\mathbf{x}}_i){]}^{1-{\delta }_i}$

Link to original

Nonparametric Methods: One Sample

Life Tables

Empirical Survival Function

Definition

Let $T_i,i=1,2,\dots ,n$ be i.i.d. random variables with the Survival Function $S(t)$, then the empirical survival function is defined as $\hat{S}(t)=\frac{1}{n}\sum _{i=1}^{n}I(T_i>t)$ Since $I(T_i>t)$ is a Bernoulli trial, $n\hat{S}(t)\sim B(n,S(t))$

Link to original

Reduced Sample Estimator

Definition

Notations

Assume intervals have equal length and let the notations:

• $I_i$: the $i$-th interval $({\tau }_{i-1},{\tau }_i]$ where ${\tau }_0=0$
• $n_i$: the number of individuals alive at the beginning of $I_i$
• $d_i$: the number of deaths during $I_i$
• $l_i$: the number of individuals censored during $I_i$
• $P_i$: $P(T>{\tau }_i|T>{\tau }_{i-1})=P(\text{surviving through}{I}_i|\text{alive at beginning of}{I}_i)$ where $T_i$ is the survival time.
• $q_i=1-P_i$

Estimation of Reduced Sample Estimator

Reduced sample method estimates the Survival Function as $\hat{S}({\tau }_k)=1-\frac{d_k^{\ast }}{n_k^{\ast }}$ where $d_k^{\ast }=\sum _{i=1}^k{d}_i$ is the number of deaths in the interval $(0,{\tau }_k)$, and $n_k^{\ast }=n_1-\sum _{i=1}^k{l}_i$ is the number of uncensored data in $(0,{\tau }_k)$

The drawback of the reduced sampling method is that it ignore the information contained in censored observations, therefore it usually a biased (under) estimator of the Survival Function $S(t)$.

Link to original

Life Table Estimator

Definition

Notations

Assume intervals have equal length and let the notations:

• $I_i$: the $i$-th interval $({\tau }_{i-1},{\tau }_i]$ where ${\tau }_0=0$
• $n_i$: the number of individuals alive at the beginning of $I_i$
• $d_i$: the number of deaths during $I_i$
• $l_i$: the number of individuals censored during $I_i$
• $P_i$: $P(T>{\tau }_i|T>{\tau }_{i-1})=P(\text{surviving through}{I}_i|\text{alive at beginning of}{I}_i)$ where $T_i$ is the survival time.
• $q_i=1-P_i$

Estimation of Life Table Estimator

The life table estimator is derived from the expression $S(t)=\prod _{i=1}^k{P}_i$ where $P_i=P(T>{\tau }_i|T>{\tau }_{i-1})$

The life table estimator estimates the Survival Function as $\hat{S}({\tau }_k)=\prod _{i=1}^k{\hat{P}}_i$ where ${\hat{P}}_i=1-\frac{d_i}{n_i^{\prime }}$ and $n_i^{\prime }=n_i-\frac{l_i}{2}$ is called the effective sample size.

We assume that, on average, those individuals who became censored during $I_i$ were at risk for half the interval

Variance of Life Table Estimator

For a given $n_i^{\prime }$, assume that $n_i^{\prime }{\hat{P}}_i\approx B(n_i^{\prime },P_i)$, where $P_i=\frac{S({\tau }_i)}{S({\tau }_{i-1})}$. Since $\ln \hat{S}({\tau }_k)=\sum _{i=1}^k\ln {\hat{P}}_i$, under the assumption of the independence of $\ln {\hat{P}}_i$‘s, the variance of the life table estimator is approximated as $\mathrm{Var}[\ln \hat{S}({\tau }_k)]\approx \sum _{i=1}^k\mathrm{Var}[\ln {\hat{P}}_i]$ Then, by the Delta Method, Greenwood’s formula is derived as $\widehat{\mathrm{Var}}(\hat{S}({\tau }_k))=\hat{S}({\tau }_k{)}^2\sum _{i=1}^k\frac{d_i}{n_i^{\prime }(n_i^{\prime }-d_i)}$

Confidence Interval for Life Table Estimator

Under the asymptotic normality of the estimator, we can use $\hat{S}(t)\pm z_{\alpha /2}\sqrt{\widehat{\mathrm{Var}}(\hat{S}(t))}$ as a $100(1-\alpha )\%$ confidence interval for $S(t)$. However, this region could take values outside $(0,1)$. To avoid this kind of problem, the log-log-transformation is used.

Log-log Transformation

To guarantee for the CI of $S(t)$ to be within $(0,1)$, use log-log transformation. Let $V(t)=\ln (-\ln S(t))$. Then, by delta-method $\widehat{\mathrm{Var}}[\hat{V}(t)]\approx \frac{\widehat{\mathrm{Var}}(\hat{S}(t))}{[\ln \hat{S}(t){]}^2\hat{S}(t{)}^2}$ i.e. CI for $V(t)$ is $\hat{V}(t)\pm z_{\alpha /2}\sqrt{\widehat{\mathrm{Var}}(\hat{V}(t))}$ and CI for $S(t)$ is calculated as $\exp \left[-\exp \left(\hat{V}(t)\pm z_{\alpha /2}\sqrt{\widehat{\mathrm{Var}}(\hat{V}(t))}\right)\right]\in (0,1)$

Examples

Consider a life table

$I_i$	$n_i$	$d_i$	$l_i$	$n_i^{\ast }=n_1-\sum _{j=1}^i{l}_j$	$d_i^{\ast }=\sum _{j=1}^i{d}_j$	${\hat{S}}_R=1-\frac{d_i^{\ast }}{n_i^{\ast }}$	$n_i^{\prime }=n_i-\frac{l_i}{2}$	${\hat{P}}_i=1-\frac{d_i}{n_i^{\prime }}$	${\hat{S}}_L=\prod _{j=1}^i{\hat{P}}_j$
$(0,1]$	$126$	$47$	$19$	$107$	$47$	$0.56$	$116.5$	$0.60$	$0.60$
$(1,2]$	$60$	$5$	$17$	$90$	$52$	$0.42$	$51.5$	$0.90$	$0.54$
$(2,3]$	$38$	$2$	$15$	$75$	$54$	$0.28$	$30.5$	$0.93$	$0.50$
$(3,4]$	$21$	$2$	$9$	$66$	$56$	$0.15$	$16.5$	$0.88$	$0.44$
$(4,5]$	$10$	$0$	$6$	$60$	$56$	$0.07$	$7$	$1.00$	$0.44$
where ${\hat{S}}_R$ is the Reduced Sample Estimator and ${\hat{S}}_R$ is the life table estimator.

The survival functions are estimated as

Link to original

Kaplan-Meier Estimator

Kaplan-Meier Estimator

Definition

Kaplan-Meier Estimator

Consider a Random Censoring case $Y_i=\min (T_i,C_i)$. Assume that $y_1\le y_2\le \cdots \le y_n$, and the distinct failure times are ${\tau }_1<{\tau }_2<\cdots <{\tau }_k$ where $k\le n$. Let $d_j:=\sum _{i=1}^{n}I(y_i={\tau }_j,{\delta }_i=1)$ be the number of deaths at ${\tau }_j$, and $n_j:=\sum _{i=1}^{n}I(y_i\ge {\tau }_j)$ be the number of alive at ${\tau }_j$ where the set$\{y_i|y_i\ge {\tau }_j\}$ is called the risk set at ${\tau }_j$. We only observe $(Y_i,{\delta }_i),i=1,\dots ,n$

The Kaplan-Meier estimator is derived from the expression $S(t)=\prod _{i=1}^k{P}_j,P_j=P(T>{\tau }_j|T>{\tau }_{j-1})$ where ${\tau }_0=0$

General Case

As an estimator of $P_j$, consider ${\hat{P}}_j=\hat{P}(T>{\tau }_j|T>{\tau }_{j-1})=\frac{\hat{P}(T>{\tau }_j)}{\hat{P}(T>{\tau }_{j-1})}=\frac{n_{j+1}/n_1}{n_j/n_1}=\frac{n_j-d_j}{n_j}=1-\frac{d_j}{n_j}$

The Kaplan-Meier estimator is defined with the estimated ${\hat{P}}_j$‘s $\hat{S}(t)=\prod _{j=1}^k{\hat{P}}_j=\prod _{\{j|{\tau }_j\le t\}}\left(1-\frac{d_j}{n_j}\right)$ The cumulative hazard function is estimated by Nelson-Aalen Estimator in the same logic $\hat{\Lambda }(t)=\sum _{\{j|{\tau }_j\le t\}}\frac{d_j}{n_j}$

No Ties Case

When there’s no tie in the observation, $y_1<y_2<\cdots <y_n$, then the failure times are equal to the observation ${\tau }_j=y_i$, death is equal to the censoring indicator $d_j={\delta }_j$, and $n_j=n-j+1$. Thus, the Kaplan-Meier estimator is defined as $\hat{S}(t)=\prod _{\{j|y_j\le t\}}(1-\frac{{\delta }_j}{n-j+1})=\prod _{\{j|{\tau }_j\le t\}}(1-\frac{1}{n-j+1}{)}^{{\delta }_j}$

Properties

Self-Consistency

An estimator $\hat{S}(t)$ is self-consistent if $\begin{array}{rl}\hat{S}(t)& =\frac{1}{n}\sum _{i=1}^n\left[1I(Y_i>t)+0I(Y_i\le t,{\delta }_i=1)+\frac{\hat{S}(t)}{\hat{S}(Y_i)}I(Y_i<t,{\delta }_i=0)\right]\\ & =\frac{1}{n}\left[N_y(t)+\sum _{i:y_i\le t}(1-{\delta }_i)\frac{\hat{S}(t)}{\hat{S}(Y_i)}\right]\end{array}$ where $N_y(t)=\sum _{i=1}^{n}I(Y_i>t)=\#(Y_i>t)$

The Kaplan-Meier estimator is the unique self-consistent estimator for $t<Y_{(n)}$ where $Y_{(n)}$ is the largest observation.

Generalized MLE

The Kaplan-Meier estimator gives the Generalized Maximum Likelihood Estimation of the Survival Function $S$.

Strong Consistency

The Kaplan-Meier estimator uniformly Almost Surely converges to $S(t)$ $\hat{S}(t)\stackrel{a.s.}{\to }S(t),\forall t\in {\mathbb{R}}^{+}$

Proof

Consider a function $S^{\ast }(t)=P(Y>t)$ and decompose it to the sum of the subsurvival functions $S_u^{\ast }(t)$ and $S_c^{\ast }(t)$. $S^{\ast }(t)=S_u^{\ast }(t)+S_c^{\ast }(t)$ where $S_u^{\ast }(t)=P(Y>t,\delta =1)$ is the uncensored case and $S_c^{\ast }(t)=P(Y>t,\delta =0)$ is the censored case.

Then, the survival function $S(t)=P(T>t)$ can be expressed as a function of the subsurvival functions. $S(t)=\Psi (S_u^{\ast },S_c^{\ast },t)$

Define the empirical subsurvival functions ${\hat{S}}_u^{\ast }(t)=\frac{1}{n}\sum _{i=1}^{n}I(Y_i>t,{\delta }_i=1)$ and ${\hat{S}}_c^{\ast }(t)=\frac{1}{n}\sum _{i=1}^{n}I(Y_i>t,{\delta }_i=0)$. The Kaplan-Meier estimator also can be expressed as a function of the empirical subsurvival functions. $\hat{S}(t)=\Psi ({\hat{S}}_u^{\ast },{\hat{S}}_c^{\ast },t)$

By Glivenko-Cantelli theorem, ${\hat{S}}_u^{\ast }(t)\stackrel{a.s.}{\to }{S}_u^{\ast }(t)$ and ${\hat{S}}_c^{\ast }(t)\stackrel{a.s.}{\to }{S}_c^{\ast }(t)$ for all $t\in {\mathbb{R}}^{+}$. Since $\Psi$ is a continuous function of $S_u^{\ast }(t)$ and $S_c^{\ast }(t)$, $\hat{S}(t)=\Psi ({\hat{S}}_u^{\ast },{\hat{S}}_c^{\ast },t)\stackrel{a.s.}{\to }\Psi (S_u^{\ast },S_c^{\ast },t)=S(t),\forall t\in {\mathbb{R}}^{+}$

Asymptotic Normality

Kaplan-Meier estimator has asymptotic normality. $\sqrt{n}(\hat{S}(t)-S(t))\stackrel{D}{\to }N\left(0,S(t{)}^2{\int }_0^t\frac{d{F}_u(X)}{(1-H(x){)}^2}\right)$ where $F_u(t)=P(T\le t,\delta =1)={\int }_0^t(1-G(x))dF(x)$, $C\sim G$, and $Y\sim H$.

The variance of the estimator is estimated by Greenwood’s formula ${\sigma }_S^2(t)\approx \hat{S}(t{)}^2\sum _{\{j|{\tau }_j\le t\}}\frac{d_j}{n_j(n_j-d_j)}$ For the no ties case, the formula is ${\sigma }_S^2(t)\approx \hat{S}(t{)}^2\sum _{\{j|{\tau }_j\le t\}}\frac{{\delta }_j}{(n-j)(n-j+1)}$

Examples

$n=8$ case where $\{\begin{array}{ll}\times & {\delta }_i=1\\ \circ & {\delta }_i=0\end{array}$

Facts

Kaplan-Meier estimator has Self-Consistency and Asymptotic Normality, and it is generalized MLE

If no censoring, Kaplan-Meier estimator is just the Empirical Survival Function.

Link to original

Hazard Function Estimators

Nelson-Aalen Estimator

Definition

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator. In a Random Censoring setting, we only observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$.

The Nelson-Aalen estimator estimates cumulative hazard function as $\hat{\Lambda }(t)=\sum _{i:y_i\le t}\frac{{\delta }_i}{(n-i+1)}$ It is derived from Kaplan-Meier Estimator

Link to original

Peterson Estimator

Definition

Suppose that $T_1,T_2,\dots ,T_n$ are i.i.d. random variables represent survival times with CDF $F$, and $C_1,C_2,\dots ,C_n$ are the censoring times. And let ${\delta }_i=\{\begin{array}{ll}1& T_i\le C_i\text{: uncensored}\\ 0& T_i>C_i\text{: censored}\end{array}$ be the censoring indicator. In a Random Censoring setting, we only observe $(Y_i,{\delta }_i),i=1,\dots ,n$ where $Y_i=\min (T_i,C_i)$.

The Peterson estimator estimates cumulative hazard function as $\hat{\Lambda }(t)=\sum _{i:y_i\le t}-\ln \left(1-\frac{{\delta }_i}{n-i+1}\right)$ It is derived from Kaplan-Meier Estimator

Link to original

Robust Estimators

Estimators for Survival Function

Definition

Mean of Survival Time

Without Censoring

$\hat{\theta }={\int }_0^{\infty }xd{F}_n(x)=\bar{x}$ where $F_n$ is the empirical CDF

With Censoring

$\hat{\theta }={\int }_0^{\infty }xd\hat{F}(x)={\int }_0^{\infty }\hat{S}(t)dt=\sum _{i=1}^n{S}_i{Y}_i$ where $\hat{F}$ is calculated by the Kaplan-Meier Estimator and $S_i$ is the jump size at $Y_i$.

The asymptotic variance of the parameter is obtained by $\mathrm{AVar}(\hat{\theta })=\frac{1}{n}{\int }_0^{\infty }\frac{1}{(1-H(S){)}^2}{\left({\int }_0^{\infty }S(t)dt\right)}^{2}d{F}_u(S)$ And the asymptotic variance is estimated by $\widehat{\mathrm{AVar}}(\hat{\theta })=\sum _{i=1}^n{\left({\int }_{y_i}^{\infty }S(t)dt\right)}^2\frac{{\delta }_i}{(n-i)(n-i+1)}$ where $(1-\widehat{H(S)}{)}^2=\frac{n-i+1}{n}$ and $d{\hat{F}}_u(y_i)=\frac{{\delta }_i}{n}$

Median of Survival Time

A reasonable estimator for $\theta$ is $\hat{\theta }={\hat{S}}^{-1}(0.5)$ where $\hat{S}$ is the Kaplan-Meier Estimator

If ${\hat{S}}^{-1}$ does not have a unique solution, then $\hat{\theta }$ is defined as the midpoint of the interval constituting of the solutions.

However, the estimator $\hat{S}(t)$ over-estimate the true parameter, so linear smooth of the estimator $\hat{\hat{S}}(t)$ is used to estimate the $\theta$ $\hat{\hat{\theta }}={\hat{\hat{S}}}^{-1}(0.5)$

The asymptotic variance of the estimated parameter is obtained by $\mathrm{AVar}(\hat{\theta })=\frac{\mathrm{AVar}(\hat{S}(\theta ))}{f^2(\theta )}$ where $\mathrm{AVar}(\hat{S}(\theta ))$ is estimated by Greenwood’s formula, and $f$ may estimated by Kernel Estimation.

Link to original

Bayes Estimator for Survival Function

Definition

The Bayes estimator of $S(t)$ with the Squared Error Loss and a Dirichlet process prior ${\mathcal{P}}_{\alpha }$ is given by ${\hat{S}}_{\alpha }(t)=\frac{\alpha (t,\infty )+N_y(t)}{\alpha (0,\infty )+n}{\prod }_{i:y_i\le t}{\left[\frac{\alpha [y_i,\infty )+(n-i+1)}{\alpha [y_i,\infty )+(n-i)}\right]}^{1-{\delta }_i}$ where the squared loss is $L(\hat{S},S)={\int }_0^{\infty }(\hat{S}(t)-S(t){)}^{2}dw(t)$ where $w$ is any non-negative non-decreasing function, the parameter $\alpha$ which is a finite non-negative measure on $(0,\infty )$, and $N_y(t)=\sum _{i=1}^{n}I(Y_i>t)=\#(Y_i>t)$.

Examples

If $\frac{\alpha (t,\infty )}{\alpha (0,\infty )}=\exp [-{\lambda }_{0}t]$, then In many cases, $\mathrm{MLE}[{\hat{S}}_{\alpha }(t)]\le \mathrm{MLE}[\hat{S}(t)]$

Link to original

Nonparametric Density Estimation

Kernel Estimation for Survival Analysis

Definition

Kaplan-Meier estimator is a step function. So it is difficult to calculate its quantile function and Density Function. The Kernel Density Estimation is used to make it smooth function.

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time with a distribution $F$, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time with a distribution $G$. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$ where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator.

Without Censoring

For the complete data, the Kernel Density Estimation is defined as $\hat{f}(t)=\frac{1}{n}\sum _{i=1}^n{K}_{\lambda }(t-Y_i)=\frac{1}{n\lambda }\sum _{i=1}^{n}K\left(\frac{t-Y_i}{\lambda }\right)$ where $K\text{s.t.}{\int }_{-\infty }^{\infty }K(t)dt=1,{\int }_{-\infty }^{\infty }{K}^2(t)dt<\infty$ is the kernel, $K_{\lambda }(x):=\frac{1}{h}K(\frac{x}{h})$ is the scaled kernel, and $\lambda$ is a smoothing parameter.

The kernel estimator for the Distribution Function is defined as $\hat{F}(t)={\int }_{-\infty }^t\hat{f}(u)du=\frac{1}{n}\sum _{i=1}^{n}W\left(\frac{t-Y_i}{\lambda }\right)$ where $W(t)={\int }_{-\infty }^{t}K(u)du$

With Censoring

For the censored data, the weights for each observation is defined as a jump size in Kaplan-Meier Estimator. $\hat{f}(t)=\sum _{i=1}^n{S}_i{K}_{\lambda }(t-Y_i)$ $\hat{F}(t)=\sum _{i=1}^n{S}_{i}W\left(\frac{t-Y_i}{h}\right)$ where $S_i$ is the jump size at $Y_i$ in Kaplan-Meier Estimator.

Thus, the Kernel Density Estimation for the Survival Function is $1-\sum _{i=1}^n{S}_{i}W\left(\frac{t-Y_i}{h}\right)$

Link to original

Nonparametric Methods: Two Samples

Gehan Test

Definition

For the first sample, let $T_1,T_2,\dots ,T_m$ be i.i.d. survival time with a distribution $F_1$, and $C_1,C_2,\dots ,C_m$ be i.i.d. censoring time with a distribution $G_1$. We can observe $(X_i,{\delta }_i),i=1,2,\dots ,m$ where $X_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator.

For the second sample, let $U_1,U_2,\dots ,U_m$ be i.i.d. survival time with a distribution $F_2$, and $D_1,D_2,\dots ,D_n$ be i.i.d. censoring time with a distribution $G_2$. We can observe $(Y_i,{ϵ}_i),i=1,2,\dots ,n$ where $Y_i=\min (U_i,D_i)$ and ${ϵ}_i=I(U_i\le D_i)$ is censoring indicator.

Gehan’s test is an extension of Signed-Rank Wilcoxon Test. The test statistic of Gehen’s test is defined as $U=\sum _{i=1}^m\sum _{j=1}^n{U}_{ij}$ where $U_{ij}=\{\begin{array}{ll}1& t_i>u_j:(x_i>y_j,{ϵ}_j=1)\text{or}(x_i=y_j,{\delta }_i=0,{ϵ}_j=1)\\ 0& \text{otherwise}\\ -1& t_i<u_j:(x_i<y_j,{\delta }_i=1)\text{or}(x_i=y_j,{\delta }_i=1,{ϵ}_j=0)\end{array}$.

Under the null hypothesis $H_0:F_1=F_2,G_1=G_2$, let $(z_i,{\zeta }_i),i=1,2,\dots ,n+m$ be the combined samples and define $U_{kl}$, $U_k^{\ast }$, and $U$ $U_{kl}=\{\begin{array}{ll}1& (z_k>z_l,{\zeta }_l=1)\text{or}(z_k=z_l,{\zeta }_k=0,{\zeta }_l=1)\\ 0& \text{otherwise}\\ -1& (z_k<z_l,{\zeta }_k=1)\text{or}(z_k=z_l,{\zeta }_k=1,{\zeta }_l=0)\end{array}$ $U_k^{\ast }=\sum _{l=1\ne k}^{m+n}{U}_{kl}$ $U=\sum _{k=1}^{m+n}{U}_k^{\ast }I(k\in I_1)$ where $I_1$ is the indicator set for sample 1.

Then, the variance of $U$ is calculated as $\mathrm{Var}(U)=\frac{mn}{(m+n)(m+n-1)}\sum _{k=1}^{m+n}(U_k^{\ast }{)}^2$ and $\frac{U}{\sqrt{\mathrm{Var}(U)}}\stackrel{D}{\to }N(0,1)$

Link to original

Hypothesis Test for a 2 by 2 Contingency Table

Definition

The hypothesis test for a $2\times 2$ contingency table is used to determine if there’s a significant association between two categorical variables

Consider a $2\times 2$ Contingency Table

	Dead	Alive
Population 1	$a$	$b$	$n_1$
Population 2	$c$	$d$	$n_2$
	$m_1$	$m_2$	$n$
i.e. $a\sim B(n_1,p_1)$ and $c\sim B(n_2,p_2)$

We want to test $H_0:p_1=p_2$. Let ${\hat{p}}_1=\frac{a}{n_1}$, ${\hat{p}}_2=\frac{c}{n_2}$, $\hat{p}=\frac{m_1}{n}$, and $B(r,n,p)=\left(\genfrac{}{}{0ex}{}{n}{r}\right)p^r(1-p{)}^{n-r}$

Uncorrelated Chi-squared Test

The test statistic is defined as $T=\frac{n(ad-bc{)}^2}{n_1{n}_2{m}_1{m}_2}\sim {\chi }^2(1)$

We reject $H_0$ if $T\ge {\chi }_{\alpha }^2(1)$

Yates’ Corrected Chi-squared Test

The test statistic is defined as $T=\frac{n(|ad-bc|-\frac{n}{2}{)}^2}{n_1{n}_2{m}_1{m}_2}\sim {\chi }^2(1)$

We reject $H_0$ if $T\ge {\chi }_{\alpha }^2(1)$

Fisher’s Exact Test

The test statistic is defined as $T=\sum _{j=x_1}^{x_2}\frac{\left(\genfrac{}{}{0ex}{}{n_1}{j}\right)\left(\genfrac{}{}{0ex}{}{n_2}{m_1-j}\right)}{\left(\genfrac{}{}{0ex}{}{n}{m_1}\right)}$ where $x_1=\max (0,m_1-n_2)$, $x_2=\min (n_1,m_1)$

We reject $H_0$ if $T\le \alpha$

Corrected Chi-squared Test

The test statistic is defined as $T=\frac{n(|ad-bc|-\frac{n}{4}{)}^2}{n_1{n}_2{m}_1{m}_2}\sim {\chi }^2(1)$

We reject $H_0$ if $T\ge {\chi }_{\alpha }^2(1)$

Liddell’s Exact Test

The test statistic is defined as $T=\sum _{x=0}^{n_1}\sum _{y=0}^{n_2}B(x,n_1,\hat{p})B(y,n_2,\hat{p})I(\left|\frac{x}{n_1}-\frac{y}{n_2}\right|\ge |a-c|)$

We reject $H_0$ if $T\le \alpha$

Exact Unconditional Test

The test statistic is defined as $T=\underset{0\le p\le 1}{\sup }\sum _{x=0}^{n_1}\sum _{y=0}^{n_2}B(x,n_1,p)B(y,n_2,p)I(|z(x,y)|\ge |z(a,c)|)$ where $z(x,y)=\frac{(\frac{x}{n_1}-\frac{y}{n_2})}{\sqrt{\frac{1}{n_1}(\frac{x}{n_1})(1-\frac{x}{n_1})+\frac{1}{n_2}(\frac{y}{n_2})(1-\frac{y}{n_2})}}$

We reject $H_0$ if $T\le \alpha$

Approximate Unconditional Test

The test statistic is defined as $T=\sum _{x=0}^{n_1}\sum _{y=0}^{n_2}B(x,n_1,\hat{p})B(y,n_2,\hat{p})I(|z^{\ast }(x,y)|\ge |z^{\ast }(a,c)|)$ where $z^{\ast }(x,y)=\frac{(\frac{x}{n_1}-\frac{y}{n_2})}{\sqrt{(\frac{1}{n_1}+\frac{1}{n_2})\hat{p}(1-\hat{p})}}$

We reject $H_0$ if $T\le \alpha$

Link to original

Mantel-Haenszel Test

Definition

Consider a sequence of $2\times 2$ contingency tables

	Dead	Alive
Treatment 1	$a_i$	$b_i$	$n_{i1}$
Treatment 2	$c_i$	$d_i$	$n_{i2}$
	$m_{i1}$	$m_{i2}$	$n_i$
where $i=1,2,\dots ,k$ is the indicator of hospital

Under the null hypothesis $H_0:p_{11}=p_{12},p_{21}=p_{22},\dots ,p_{k1}=p_{k2}$, where $p_{i1}=\frac{a_i}{n_{i1}}$ and $p_{i2}=\frac{c_i}{n_{i2}}$, the test statistic for Mantel-Haenszel test is defined as $MH=\frac{\sum _{i=1}^k(a_i-E(a_i))}{\sqrt{\sum _{i=1}^k\mathrm{Var}(a_i)}}\stackrel{D}{\to }N(0,1)$ where $E(a_i)=\frac{m_i{n}_i}{n_i}$ and $\mathrm{Var}(a_i)=\frac{n_{i1}{n}_{i2}{m}_{i1}{m}_{i2}}{n_i^2(n_i-1)}$.

Link to original

Log-Rank Test

Definition

The log-rank test is a test to compare the survival functions of two samples. The test uses the sequence of Mantel-Haenszel statistics of the $2\times 2$ tables at each uncensored event time.

Examples

$z$	$n$	$m_1$	$n_1$	$a$	$E(a)$	$a-E(a)$	$\mathrm{Var}(a)$
$3$	$10$	$1$	$5$	$1$	$0.50$	$0.50$	$0.250$
$5$	$9$	$1$	$4$	$1$	$0.44$	$0.56$	$0.247$
$7$	$8$	$1$	$3$	$1$	$0.38$	$0.62$	$0.234$
$12$	$6$	$1$	$1$	$0$	$0.17$	$-0.17$	$0.139$
$18$	$5$	$1$	$1$	$1$	$0.20$	$0.80$	$0.160$
$19$	$4$	$1$	$0$	$0$	$0$	$0$	$0$
$20$	$3$	$1$	$0$	$0$	$0$	$0$	$0$
The Mantel-Haenszel statistic of the data is calculated as
$MH=\frac{\sum _{i=1}^k[a_i-E(a_i)]}{\sqrt{\sum _{i=1}^k\mathrm{Var}(a_i)}}\approx \frac{2.31}{\sqrt{1.03}}\approx 2.26$
and the p-value is $0.12$

Link to original

Tarone-Ware Test

Definition

The Tarone-Ware test is the generalization of the Mantel-Haenszel Test. $T=\frac{\sum _{i=1}^k{w}_i(a_i-E(a_i))}{\sqrt{\sum _{i=1}^k\mathrm{Var}(a_i)}}\stackrel{D}{\to }N(0,1)$ where $w_i$ is the weight for each table.

If $w_i=1$, then it is MH Statistic, if $w_i=n_i$ then it is Gehan statistic, and if $w_i=\sqrt{n_i}$ then it is Tarone-Ware statistic.

Link to original

Nonparametric Methods: K Samples

Generalized Gehan Test

Definition

For the $i$-th sample, let $T_{i1},T_{i2},\dots ,T_{i{n}_i}$, where $i=1,2,\dots ,K$, be i.i.d. survival time with a distribution $F_i$, and $C_{i1},C_{i2},\dots ,C_{i{n}_i}$ be i.i.d. censoring time with a distribution $G_i$. We can observe $(X_{ij},{\delta }_{ij}),j=1,2,\dots ,n_i$ where $X_{ij}=\min (T_{ij},C_{ij})$ and ${\delta }_{ij}=I(T_{ij}\le C_{ij})$ is the censoring indicator.

The generalized Gehan’s test is an extension of Gehan Test used for more than two sample case. Under the null hypothesis $H_0:F_1=F_2=\cdots =F_K$, the test statistic of generalized Gehen’s test is defined as $W_i=\sum _{i^{\prime }=1\ne i}^K\sum _{j=1}^{n_i}\sum _{j^{\prime }=1}^{n_{i^{\prime }}}U((X_{ij},{\delta }_{ij}),(X_{i^{\prime }{j}^{\prime }},{\delta }_{i^{\prime }{j}^{\prime }}))$ where $U$ is the statistic of thet Gehan Test.

Then, $W=(W_1,W_2,\dots ,W_K)\to N_K(0,N^3{\Sigma }^{\ast })$ where $N=\sum _{i=1}^K{n}_i$, ${\hat{\Sigma }}^{\ast }=\frac{1}{N(N-1)}\sum _{i=1}^K\sum _{j=1}^{n_i}(W_{ij}^{\ast }{)}^2\left(\begin{array}{ccc}{n}_1(N-n_1)& & -n_i{n}_j\\ & \ddots & \\ & & n_K(N-n_K)\end{array}\right)$, and $w_{ij}^{\ast }=\sum _{i^{\prime }=1\ne i}^K\sum _{j^{\prime }=1\ne j}^{K}U((X_{ij},{\delta }_{ij}),(X_{i^{\prime }{j}^{\prime }},{\delta }_{i^{\prime }{j}^{\prime }}))$

Link to original

Generalized Mantel-Haenszel Test

Definition

For the $i$-th sample, let $T_{i1},T_{i2},\dots ,T_{i{n}_i}$, where $i=1,2,\dots ,K$, be i.i.d. survival time with a distribution $F_i$, and $C_{i1},C_{i2},\dots ,C_{i{n}_i}$ be i.i.d. censoring time with a distribution $G_i$. We can observe $(X_{ij},{\delta }_{ij}),j=1,2,\dots ,n_i$ where $X_{ij}=\min (T_{ij},C_{ij})$ and ${\delta }_{ij}=I(T_{ij}\le C_{ij})$ is the censoring indicator.

The generalized Mantel-Haenszel test is an extension of Mantel-Haenszel Test used for more than two sample case.

Let $(z_i,{\zeta }_i),i=1,2,\dots ,N$ be the combined samples

For each uncensored time point, construct a $2\times K$ table.

	$1$	$2$	$\dots$	$K$
Dead	$a_{i1}$	$a_{i2}$	$\dots$	$a_{iK}$	$m_{i1}$
Alive			$\dots$		$m_{i2}$
	$n_{i1}$	$n_{i2}$	$\dots$	$n_{i}K$	$N_i$

Under the null hypothesis $H_0:F_1=F_2=\cdots =F_K$, the test statistic of generalized Mantel-Haenszel test is defined as $MH=(a_{(-1)}-E(a_{(-1)}){)}^t{\Sigma }_{(1)}(a_{(-1)}-E(a_{(-1)}))\stackrel{D}{\to }{\chi }^2(K-1)$ where $a_{(-1)}$, ${\Sigma }_{(-1)}$ are the $K-1$-length vector and $(K-1)\times (K-1)$ matrix, calculated with the first population is deleted data (corner-point constraint).

The $a$ and $\Sigma$ is defined as $a-E(a)=\sum _{i=1}^K(a_i-E(a_i))$ where $E(a_i)=(\frac{m_{i1}{n}_{i1}}{N_i},\dots ,\frac{m_{iK}{n}_{iK}}{N_i})$ $\Sigma =\sum _{i=1}^K\mathrm{Cov}(a_i)$ where $\mathrm{Cov}(a_i)=\frac{m_{i1}{m}_{i2}}{N_i-1}\left(\begin{array}{ccc}\frac{n_{i1}}{N_i}(1-\frac{n_{i1}}{N_i})& & -\frac{n_{ik}{n}_{il}}{N_i^2}\\ & \ddots & \\ & & \frac{n_{iK}}{N_i}(1-\frac{n_{iK}}{N_i})\end{array}\right)$

Link to original

Nonparametric Methods: Regression

Cox Proportional Hazards Model

Cox Proportional Hazards Model

Definition

Cox proportional hazards model assume that covariates affect the Hazard Function.

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator, and have covariates ${\mathbf{x}}_i=(x_{i1},\dots ,x_{ip})$. Then, the Cox proportional hazards model is defined as $\lambda (t:\mathbf{x})={\lambda }_0(t)\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })$ where ${\lambda }_0$ is called the baseline hazard function, i.e. hazard at $\mathbf{x}=0$

Conditional Likelihood

Let $Y_1<Y_2<\cdots <Y_n$ (no ties case), and ${\mathcal{R}}_i$ be the risk set. For each uncensored time $Y_i$, $P\{\text{a death in}[y_i,y_i+\Delta )|{\mathcal{R}}_i\}\approx \sum _{j\in {\mathcal{R}}_i}{\lambda }_0(y_j)\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })\Delta$ Therefore, $P\{\text{a death of }i\text{ at time}{y}_i|\text{one death in }{\mathcal{R}}_i\text{ at time }{y}_i\}=\frac{{\lambda }_0(y_i)\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })\Delta }{\sum _{j\in {\mathcal{R}}_i}{\lambda }_0(y_i)\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })\Delta }\approx \frac{\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })}$ Taking the product of these conditional probabilities gives a conditional likelihood $L_C(\mathbf{\beta })={\prod }_{i:I_u}\frac{\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })}={\prod }_{i=1}^n{\left[\frac{\exp ({\mathbf{x}}_{\mathbf{i}}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })}\right]}^{{\delta }_i}$ where $I_u$ is the indicator set for uncensored samples.

The $L_C(\mathbf{\beta })$ is not a likelihood. However, Cox suggested treating the conditional likelihood as an ordinary likelihood to find the Maximum Likelihood Estimation.

Since there’s no analytic solution for the MLE, iterative methods such as Newton–Raphson method is used to estimate the coefficient $\mathbf{\beta }$.

The hazard ratio $\exp ({\hat{\beta }}_j)$ represents the relative change in Hazard Rate for a one-unit increase in the covariate $x_j$.

Goodness-of-Fit Test

For testing the null hypothesis $H_0:\mathbf{\beta }=0$, Cox suggested the Rao Test. ${\left(\frac{\partial \ln L_C(0)}{\partial \mathbf{\beta }}\right)}^{⊺}{\left(-\frac{{\partial }^2\ln L_C(0)}{\partial {\mathbf{\beta }}^2}\right)}^{-1}\left(\frac{\partial \ln L_C(0)}{\partial \mathbf{\beta }}\right)\stackrel{a}{\sim }{\chi }^2(p)$

Asymptotic Normality of MLE

$\hat{\mathbf{\beta }}\stackrel{D}{\to }{N}_p(\beta ,I^{-1}(\mathbf{\beta }))$ where $I(\mathbf{\beta })$ is the observed Fisher Information

Estimation of Survival Function

Under the Cox proportional hazards model, $S(t;\mathbf{x})=S_0(t{)}^{\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })}$ To estimate $S(t;\mathbf{x})$, we can use $\hat{\mathbf{\beta }}$ for $\mathbf{\beta }$ but we still need to estimate $S_0(t)$, ${\lambda }_0(t)$, or ${\Lambda }_0(t)$.

Breslow suggested the estimators of ${\lambda }_0(t)$ and $S_{0,B}(t)$ as ${\hat{\lambda }}_{0,B}=\frac{1}{(y_{u_i}-y_{u_{i-1}})\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\hat{\mathbf{\beta }})}$ If $Y_{u_{i-1}}<t<Y_{u_i}$

${\hat{S}}_{0,B}(t)={\prod }_{i:y_i\le t}\left(1-\frac{{\delta }_i}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\hat{\mathbf{\beta }})}\right)$ If $\beta =0$, then ${\hat{S}}_{0,B}$ is the Kaplan-Meier Estimator

It has a few drawbacks

• ${\hat{S}}_{0,B}(t)\ne \exp (-{\hat{\Lambda }}_0(t))$
• ${\hat{S}}_{0,B}(t)$ can take negative values.

Tsiatis suggested a non-negative version of ${\hat{S}}_{0,B}(t)$ ${\hat{S}}_{0,T}(t)=\exp (-{\hat{\Lambda }}_{0,T}(t))$ where ${\hat{\Lambda }}_{0,T}(t)=\sum _{y_i\le t}\frac{{\delta }_i}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j^{⊺}\hat{\mathbf{\beta }})}$

Link suggested using the linear smooth of ${\hat{\Lambda }}_{0,T}(t)$.

Discrete on Grouped Data

When data is discrete or grouped, there are ties at each failure. Denote the ordered discrete failure time by $Y_1<Y_2<\cdots <Y_r$ and let ${\mathcal{R}}_i$ be the risk set at $Y_i^{-}$, ${\mathcal{D}}_i$ be the death set at $Y_i$, and $d_i=\#(D_i)$.

Cox suggested combining the all possible permutations. However, it is computationally infeasible. $L_C={\prod }_{i=1}^r\left[\frac{\prod _{j\in {\mathcal{D}}_i}{\psi }_j}{\sum _{{\mathcal{D}}_i^{\ast }}\prod _{j\in {\mathcal{D}}_i^{\ast }}{\psi }_j}\right]$ where ${\psi }_j=\exp ({\mathbf{x}}_j^{⊺}\mathbf{\beta })$, and ${\mathcal{D}}_i^{\ast }$ is the size $d_i$ subset of ${\mathcal{R}}_i$

Peto suggested an alternative likelihood that instead of all possible permutations, use the same contribution. $L_C={\prod }_{i=1}^r\left[\frac{\prod _{j\in {\mathcal{D}}_i}{\psi }_j}{(\sum _{j\in {\mathcal{R}}_i}{\psi }_j{)}^{d_i}}\right]$

Time Dependent Covariates

In the case, the covariate depends on time. We observe ${\mathbf{x}}_i(t)$ and the conditional likelihood defined as $L_C(\mathbf{\beta })={\prod }_{i:I_u}\frac{\exp ({\mathbf{x}}_i(y_i{)}^{⊺}\mathbf{\beta })}{\sum _{j\in {\mathcal{R}}_i}\exp ({\mathbf{x}}_j(y_i{)}^{⊺}\mathbf{\beta })}$ where $I_u$ is the indicator set for uncensored samples, and ${\mathcal{R}}_i$ is the risk set.

Facts

Any two individuals have hazard functions that are constant multiples of the one another.

The Survival Function of the Cox proportional hazard model is a family of Lehmann alternatives. $\forall S\in \mathcal{S},\exists S_0\in \mathcal{S}s.t.S=S_0^{\gamma }$ where $\gamma \in {\mathbb{R}}^{+}$.

If $p=1$, $x_i=I(i\in I_1)$, where $I_1$ is the indicator set for sample 1, and there are no ties, then Cox test is exactly equal to the Mantel-Haenszel Test.

Link to original

Linear Models

Accelerated Life Model

Definition

Consider the random variable $T_0$ represents survival time with $\mathbf{x}=0$ Hazard Function with $f_0(t)$, ${\lambda }_0(t)$, and $S_0(t)$ And assume that the survival time of individual with covariate $\mathbf{x}$ is defined as $T_{\mathbf{x}}=T_0\exp ({\mathbf{x}}^{⊺}\mathbf{\beta })$ If ${\mathbf{x}}^{⊺}\mathbf{\beta }<0$ then the covariate $\mathbf{x}$ accelerates the time to failure. The model based on this assumption is called accelerated failure time (AFT) model.

Under AFT model $S(t:\mathbf{x})=S_0(t\exp (-{\mathbf{x}}^{⊺}\mathbf{\beta }))$ $\lambda (t:\mathbf{x})={\lambda }_0(t\exp (-{\mathbf{x}}^{⊺}\mathbf{\beta }))\exp (-{\mathbf{x}}^{⊺}\mathbf{\beta })$

Let $Y=\ln T_{\mathbf{x}}$, then $E(Y)=E(\ln T_0)+{\mathbf{x}}^{⊺}\mathbf{\beta }\equiv \alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }$ $Y=\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }+ϵ\equiv E(\ln T_0)+{\mathbf{x}}^{⊺}\mathbf{\beta }+\underset{ϵ}{\underbrace{\ln T_0-E(\ln T_0)}}$ where $ϵ=\ln T_0-E(\ln T_0)$.

Assume that $ϵ=\sigma W$, where $W$ is a Random Variable represents error term. Then, the AFT model becomes $Y\equiv \ln T_{\mathbf{x}}=\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }+\sigma W$

Relationship between $Y=\ln T_{\mathbf{x}}$ and $W$. $F_T(t)=P_W(\frac{\ln t-\mu }{\sigma })$ $f_T(t)=\frac{1}{\sigma t}{f}_W(\frac{\ln t-\mu }{\sigma })$ where $\mu =\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta }$

If $W\sim N(0,1)$ then $\ln T\sim N(\mu ,{\sigma }^2)$, and if $W\sim \mathrm{Gumbel}(0,1)$, then $T\sim \mathrm{Exp}(\lambda )$ where $\lambda =\exp [-(\alpha +{\mathbf{x}}^{⊺}\mathbf{\beta })]$.

Link to original

Miller Estimator

Definition

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator

Suppose a Simple Linear Regression model $Y_i=\alpha +\beta X_i+{ϵ}_i,i=1,2,\dots ,n$

With no censoring present, the least squares estimators of the parameters $\alpha ,\beta$ are obtained by minimizing $\frac{1}{n}\sum _{i=1}^n{ϵ}_i^2=\frac{1}{n}\sum _{i=1}^n(T_i-\alpha -\beta x_i{)}^2={\int }_{-\infty }^{\infty }{z}^{2}d{F}_n(z)$ where $F_n(z):=\frac{1}{n}\sum _{i=1}^{n}I({ϵ}_i\le z)$ is the Empirical Distribution Function of $z_1,z_2,\dots ,z_n$ where $z_i=y_i-\alpha -\beta x_i$.

With censoring present, Miller proposed to minimize ${\int }_{-\infty }^{\infty }{z}^{2}d{\hat{F}}_n(z)=\sum _{i=1}^n{\hat{w}}_i(\beta )(Y_i-\alpha -\beta x_i{)}^2$ where $\hat{F}$ is the Kaplan-Meier Estimator based on $(z_i,{\delta }_i)$ and the weights ${\hat{w}}_i(\beta )$ is its jump size.

If the last observation is censored, then $\sum _{i=1}^n{\hat{w}}_i(\beta )<1$. Hence, change the last observation to be uncensored, so that $\sum _{i=1}^n{\hat{w}}_i(\beta )=1$.

Link to original

Buckley-James Estimator

Definition

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator

If we can observe the true survival time $T_i$, we can make a model $E(T_i)=\alpha +\beta x_i$ However, we can’t observe $T_i$, but only censored $Y_i$, and $E(Y_i)\ne \alpha +\beta x_i$ Buckley and James proposed an Unbiased Estimator for $\alpha +\beta x_i$ $Y_i^{\ast }=Y_i{\delta }_i+E(T_i|T_i>Y_i)(1-{\delta }_i)$ Since we also can not observe $y_i^{\ast }$, we estimate it again. ${\hat{y}}_i^{\ast }=y_i{\delta }_i+\left[\hat{\beta }{x}_i+\frac{\sum _{k:{\hat{z}}_k>{\hat{z}}_i}{\hat{w}}_k(\hat{\beta }){\hat{z}}_k}{1-\hat{F}({\hat{z}}_i)}\right](1-{\delta }_i)$ where ${\hat{z}}_i=y_i-\hat{\beta }{x}_i$, $\hat{F}$ is the Kaplan-Meier Estimator based on $({\hat{z}}_i,{\delta }_i)$ and the weights ${\hat{w}}_i(\hat{\beta })$ is its jump size.

The variance of the estimator is estimated by $\widehat{\mathrm{Var}}(\hat{\beta })=\frac{{\hat{\sigma }}_u^2}{\sum _{i:{\delta }_i=1}(x_i-{\bar{x}}_u{)}^2}$ where ${\hat{\sigma }}_u^2=\frac{1}{n_u-2}\sum _{i:{\delta }_i=1}[y_i-{\bar{y}}_u-\hat{\beta }(x_i-{\bar{x}}_u){]}^2$, $n_u=\sum _{i=1}^n{\delta }_i$, ${\bar{x}}_u=\sum _{i:{\delta }_i=1}{x}_i$, and ${\bar{y}}_u=\sum _{i:{\delta }_i=1}{y}_i$

Link to original

Koul-Susarla-Van Ryzin Estimator

Definition

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time. We can observe $(Y_i,{\delta }_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$ and ${\delta }_i=I(T_i\le C_i)$ is censoring indicator

Suppose a Simple Linear Regression model $E(T_i)=\alpha +\beta x_i$

Koul-Susarla-Van Ryzin proposed an Unbiased Estimator for $\alpha +\beta x_i$ $Y_i^{\ast }=\frac{{\delta }_i{Y}_i}{1-G(Y_i)}$ Since $G(Y_i)$ is unknown, it should be estimated. Authors suggested to use Kaplan-Meier Estimator of $G$ with data $(Y_i,1-{\delta }_i),i=1,2,\dots ,n$ as an estimator.

Link to original

Goodness-of-Fit Tests

Graphical Validity Tests for Survival Model

Definition

If the selected model holds, a plot of the data resembles a straight line, and if model fails, a plot resembles a curved line.

There are two types of plots, survival plots $(\hat{S}(t),t)$ and hazard plots $(\hat{\Lambda }(t),t)$.

One Sample

Exponential Distribution

$\ln S(t)=-\lambda t$ $\Lambda (t)=\lambda t$

Weibull Distribution

$\ln \Lambda (t)=\ln \lambda +\alpha \ln t$

Log-Normal Distribution

${\Phi }^{-1}(1-S(t))=-\frac{\mu }{\sigma }+\frac{1}{\sigma }\ln t$ where ${\Phi }^{-1}$ is the Probit.

Two to K Samples

For parametric models, repeat one sample methods on each sample.

The validity of the Cox Proportional Hazards Model can be checked by the Lehmann alternatives property of the Survival Function. $\ln S_i(t)={\gamma }_{ij}\ln S_j(t)$

Regression

Linear Model

Ordinary residual $e_i=y_i-\hat{\beta }{x}_i$ may be used in model checking.

Cox Proportional Hazard Model

$\ln \hat{S}(t)\approx -\lambda t$ where $\hat{S}$ is the Kaplan-Meier Estimator based on $({\hat{\Lambda }}_i,{\delta }_i),i=1,2,\dots ,n$ and ${\hat{\Lambda }}_i=\exp ({\mathbf{x}}_i^{⊺}\hat{\mathbf{\beta }}){\int }_0^{Y_i}\widehat{{\lambda }_0}(u)du$.

Also, the estimated cumulative hazard function and the covariates shouldn’t have any systematic pattern.

Link to original

Goodness-of-Fit Tests for Survival Model

Definition

No Censoring Case

Kolmogorov–Smirnov Test

Definition

The Kolmogorov–Smirnov test (KS test) is a non-parametric test for the equality of continuous, distribution functions.

The Kolmogorov–Smirnov test statistic for a given CDF $F$ is defined as $KS={\sup }_t|F_n(t)-F(t)|={\max }_{1\le i\le n}\left[\frac{1}{n}-F(t_i),F(t_i)-\frac{i-1}{n}\right]$ where $F_n$ is the Empirical Distribution Function based on the i.i.d. random variables $X_i,i=1,2,\dots ,n$.

Link to original

Cramer-von Mises Test

Definition

The Cramer-von Mises test is a non-parametric test for the equality of continuous, distribution functions.

The Cramer-von Mises test statistic for a given CDF $F$ is defined as $CV=n{\int }_{-\infty }^{\infty }[F_n(t)-F(t){]}^{2}dF(t)=\frac{1}{12n}+\sum _{i=1}^n{\left[\frac{2i-1}{2n}-F(t_i)\right]}^2$ where $F_n$ is the Empirical Distribution Function based on the i.i.d. random variables $X_i,i=1,2,\dots ,n$.

Link to original

Censoring Case

Generalized Kolmogorov–Smirnov Test

The generalized Kolmogorov–Smirnov test uses Kaplan-Meier Estimator instead of Empirical Distribution Function used for Kolmogorov–Smirnov Test $GKS=\sqrt{n}|\hat{F}(t)-F(t)|$ where $\hat{F}$ is the Empirical Distribution Function based on the i.i.d. random variables $X_i,i=1,2,\dots ,n$.

Generalized Cramer-von Mises Test

The generalized Cramer-von Mises test uses Kaplan-Meier Estimator instead of Empirical Distribution Function used for Cramer-von Mises Test $GCV=n{\int }_0^1[\hat{F}(t)-F(t){]}^{2}dt$ where $\hat{F}$ is the Empirical Distribution Function based on the i.i.d. random variables $X_i,i=1,2,\dots ,n$.

Link to original

Miscellaneous Topics

Multivariate Survival Model

Kinds

Copula Model

Definition

Let $T_1,T_2,\dots ,T_n$ be i.i.d. survival time with CDF $F(t_1,t_2,\dots ,t_k)$, PDF $f(t_1,t_2,\dots ,t_k)$, and Survival Function $S(t_1,t_2,\dots ,t_k)$, $S_j(t_j)=P(T_j>t_j)$ be marginal survival function, $F_j(t_j)=P(T_j\le t_j)$ be marginal CDF, and $C_1,C_2,\dots ,C_n$ be i.i.d. censoring time.

Copula model define the joint survival function as the copula function whose arguments are marginal survival functions.

When $k=2$, $S(t_1,t_2)=K(S_1(t_1),S_2(t_2))$ where $K$ is the copula function with uniform marginals.

Copula Functions

Clayton’s copula function $S(t_1,t_2)=(S_1(t_1{)}^{-\varphi }+S_2(t_2)-1{)}^{-1/\varphi }$ where $\varphi \in {\mathbb{R}}^{+}$

Crowder’s copula function $S(t_1,t_2)=(1+{\lambda }_1{t}_1^{{\gamma }_1}+{\lambda }_2{t}_2^{{\gamma }_2}{)}^{-\nu }$

Hougaard’s copula function $S(t_1,t_2)=(1+\varphi {\Lambda }_1(t_1)+\varphi {\Lambda }_2(t_2){)}^{-1/\varphi }$

Link to original

Competing Risks Model

Definition

Competing risks (multiple modes of failures) model is designed to accommodate the multiple causes to the same event.

Let $(T_i,K_i),i=1,2,\dots ,n$ be the survival time and failure mode (competing risk type) for each individual, where $K_i\in \{1,2,\dots ,k\}$

The mode-specific hazard function is defined as ${\lambda }_j(t)={\lim }_{\Delta t\to 0}\frac{P(T<t+\Delta t,K=j|T>t)}{\Delta t}$ where $j=1,2,\dots ,k$

The marginal hazard function and marginal cumulative hazard function are defined as $\lambda (t)=\sum _{j=1}^k{\lambda }_j(t)$ $\Lambda (t)=\sum _{j=1}^k{\Lambda }_j(t)$

Likelihood

We can observe $(Y_i,{\delta }_i,K_i),i=1,2,\dots ,n$, where $Y_i=\min (T_i,C_i)$, ${\delta }_i=I(T_i\le C_i)$ is censoring indicator, and $K_i$ is failure mode.

The likelihood function is defined as $L={\prod }_{i=1}^n{\prod }_{j=1}^k{f}_j(y_i{)}^{{\delta }_{ij}}S(y_i{)}^{1-{\delta }_i}$ where ${\delta }_{ij}=I(T_i\le C_i,K_i=j)$

Estimators

The Nelson-Aalen Estimator for competing risk model is defined as ${\hat{\Lambda }}_j(t)=\sum _{i:y_i\le t}\frac{{\delta }_{ij}}{n_i}$ where $j=1,2,\dots ,k$

The survival function is estimated by $\hat{S}(t)=\exp \left[-\sum _{j=1}^k{\hat{\Lambda }}_j(t)\right]$

The sub-distribution function is estimated by ${\hat{F}}_j(t)=\sum _{i:y_i\le t}\hat{S}(t_i)\frac{{\delta }_{ij}}{n_i}$ where $j=1,2,\dots ,k$

Link to original

Examples

When $k=2$, ${\lambda }_j(t)={\lim }_{\Delta t\to 0}\frac{P(T_j<t+\Delta t|T_1>t,T_2>t)}{\Delta t}$ where $j=1,2$

${\lambda }_{1|2}(t_1|t_2)={\lim }_{\Delta t\to 0}\frac{P(T_1<t_1+\Delta t|T_1>t_1,T_2=t_2)}{\Delta t}$ where $t_1>t_2$

${\lambda }_{2|1}(t_2|t_1)={\lim }_{\Delta t\to 0}\frac{P(T_2<t_2+\Delta t|T_1=t_1,T_2>t_2)}{\Delta t}$ where $t_1<t_2$

Link to original

Basic Issues in Clinical Trials

Observational Study

Definition

An observational study draws inferences from a sample to population where the independent variable is not under the control of the researcher.

Types

Case-Control Study

Definition

A case-control study is a retrospective study that compares two groups of people: those with a specific outcome (cases) and similar people without the outcome (control)

Examples

Suppose that researchers want to study the relationship between smoking and lung cancer. They identify 100 lung cancer patients (cases) and 100 matched individuals without lung cancer (controls). They then collect data on past smoking habits for both groups and compare the prevalence of smoking between cases and controls.

Link to original

Cohort Study

Definition

A cohort study is a prospective study that follows a group of individuals (cohort) over time to determine the incidence of a specific outcome.

Examples

Suppose that researchers want to study the relationship between smoking and lung cancer. They follow a group of 10,000 people for 20 years, comparing smokers to non-smokers to determine the incidence of lung cancer.

Link to original

Cross-Sectional Study

Definition

A cross-sectional study analyzes data from a population at a specific point in time. It provides a snapshot of the prevalence in a population

Examples

Suppose that researchers want to study the relationship between smoking and lung cancer. They survey 1,000 people in a city, collecting data on their current smoking habits and presence of lung cancer.

Link to original

Link to original

Relative Risk

Definition

Consider a $2\times 2$ Contingency Table

	Event	Non-event
Group 1	$n_{11}$	$n_{12}$	$n_{1.}$
Group 2	$n_{21}$	$n_{22}$	$n_{2.}$
	$n_{.1}$	$n_{.2}$	$n_{..}$

Point Estimation

The relative risk (RR) is estimated by $RR=\frac{\text{experimental event rate}}{\text{control event rate}}=\frac{n_{11}/n_{1.}}{n_{21}/n_{2.}}$

Confidence Interval

The confidence interval for relative risk is estimated by Delta Method. $\widehat{\mathrm{Var}}(\ln RR)=\frac{1}{n_{12}}(1-\frac{n_{11}}{n_{1.}})+\frac{1}{n_{21}}(1-\frac{n_{21}}{n_{2.}})$ The $100(1-\alpha )\%$ confidence interval for $RR$ is defined as $\exp [\ln RR\pm z_{\alpha /2}\sqrt{\widehat{\mathrm{Var}}(\ln RR)}]$

Facts

Relative risk must be used in Cohort Study or experimental study, can not be used for Case-Control Study.

Link to original

Odds Ratio

Definition

Consider a $2\times 2$ Contingency Table

	Event	Non-event
Group 1	$n_{11}$	$n_{12}$	$n_{1.}$
Group 2	$n_{21}$	$n_{22}$	$n_{2.}$
	$n_{.1}$	$n_{.2}$	$n_{..}$

Point Estimation

The odds ratio (OR) is estimated by $OR=\frac{\text{odds in case group}}{\text{odds in control group}}=\frac{(n_{11}/n_{.1})/(n_{21}/n_{.1})}{(n_{12}/n_{.2})/(n_{22}/n_{.2})}=\frac{n_{11}{n}_{22}}{n_{21}{n}_{12}}$

Confidence Interval

The confidence interval for odds ratio is estimated by Delta Method. $\widehat{\mathrm{Var}}(\ln OR)=\frac{1}{n_{11}}+\frac{1}{n_{12}}+\frac{1}{n_{21}}+\frac{1}{n_{22}}$ The $100(1-\alpha )\%$ confidence interval for $OR$ is defined as $\exp [\ln OR\pm z_{\alpha /2}\sqrt{\widehat{\mathrm{Var}}(\ln OR)}]$

Facts

Odds ratio can be used for Case-Control Study.

Link to original

Tests of Association

Kinds

Independence Test for Two Discrete Variables

Definition

Let category variables $A_1,\dots ,A_r$ and $B_1,\dots ,B_c$ and consider a null hypothesis $H_0:A,B$ are independent. Then, the test Statistic, which follows Chi-squared Distribution, is defined as $Q=\sum _{j=1}^b\sum _{i=1}^a\frac{(X_{ij}-n{\hat{p}}_{ij}{)}^2}{n{\hat{p}}_{ij}}\stackrel{D}{\to }{\chi }^2((r-1)(c-1))$ where $n=\sum _j\sum _i{X}_{ij}$, ${\hat{p}}_{ij}={\hat{p}}_{i.}{\hat{p}}_{.j}=\frac{X_{i.}}{n}\frac{X_{.j}}{n}$, $X_i=\sum _j{X}_{ij}$, $X_j=\sum _i{X}_{ij}$

Link to original

Fisher’s Exact Test

The test statistic is defined as $T=\sum _{j=x_1}^{x_2}\frac{\left(\genfrac{}{}{0ex}{}{n_1}{j}\right)\left(\genfrac{}{}{0ex}{}{n_2}{m_1-j}\right)}{\left(\genfrac{}{}{0ex}{}{n}{m_1}\right)}$ where $x_1=\max (0,m_1-n_2)$, $x_2=\min (n_1,m_1)$

We reject $H_0$ if $T\le \alpha$

Link to original

McNemar's Test

Definition

McNemar’s test is used to analyze paired nominal data (same samples are used for both conditions), particularly in before-and-after studies or matched-pair designs.

The test statistic is defined as $Q=\frac{(n_{12}-n_{21}{)}^2}{n_{12}+n_{21}}\sim {\chi }_2(1)$

Examples

	After treatment / No insomnia	After treatment / Insomnia
Before treatment / No insomnia	45	15
Before treatment / Insomnia	25	15

Consider a null hypothesis $H_0:$ The treatment is not effective. $Q=\frac{(15-25{)}^2}{15+25}=2.5<{\chi }_{0.05}^2(1)$ We can not reject the null hypothesis.

Link to original

Link to original

Confusion Matrix

Definition

	Predicted Positive (PP)	Predicted Negative (PN)
Actual Positive	True Positive (TP)	False Negative (FN)
Actual Negative	False Positive (FP)	True Negative (TN)

Metrics

Accuracy

Definition

$\frac{\text{TP}+\text{TN}}{\text{TP}+\text{TN}+\text{FP}+\text{FN}}$

Link to original

Recall

Definition

$\frac{\text{TP}}{\text{TP}+\text{FN}}$

Recall, Sensitivity, or True positive rate means that the rate of correctly predicted cases out of all the actual positive cases..

Link to original

Precision

Definition

$\frac{\text{TP}}{\text{TP}+\text{FP}}$

Precision means that the rate of actually positive cases out of cases predicted as positive.

Link to original

Specificity

Definition

$\frac{\text{TN}}{\text{FP}+\text{TN}}$

Specificity means that the rate of correctly predicted cases out of all the actual negative cases.

Link to original

Type 1 Error

Definition

$\frac{\text{FP}}{\text{FP}+\text{TN}}$

Link to original

Type 2 Error

Definition

$\frac{\text{FN}}{\text{TP}+\text{FN}}$

Link to original

F-Score

Definition

F1 Score

$F_1=\frac{2}{{\text{recall}}^{-1}+{\text{precision}}^{-1}}=2\frac{\text{precision}\cdot \text{recall}}{\text{precision}+\text{recall}}$ The harmonic mean of Precision and Recall.

F-beta score

$F_{\beta }=(1+{\beta }^2)\frac{\text{precision}\cdot \text{recall}}{({\beta }^2\cdot \text{precision})+\text{recall}}$ where $\beta >0$

Recall is considered $\beta$ times as important as Precision.

Link to original

Positive Predictive Value

Definition

$\frac{\text{Sensitivity}\times \text{Prevalence}}{(\text{Sensitivity}\times \text{Prevalence})+((1-\text{Specificity})\times (1-\text{Prevalence}))}$ where $\text{Prevalence}=\frac{\text{positive cases}}{\text{total population}}$

Positive predictive value (in medical statistics and epidemiology) means that the rate of actually positive cases out of cases predicted as positive.

Link to original

Link to original

Receiver Operating Characteristic Curve

Definition

A receiver operating characteristic curve (ROC curve) is a plot of the True Positive Rate and False Negative Rate at each threshold setting.

AUC

The area under the ROC curve is called AUC (Area under curve)

Link to original