STAT3655 Survival Analysis

Yu Gu, PhD
Assistant Professor

Department of Statistics & Actuarial Science

The University of Hong Kong

1 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

2 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

3 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

4 / 46
Survival Data

Survival (failure) times: times to the occurrence of a given event (failure)

measured from a well-defined starting point.
↭ death times of patients enrolled in a clinical trial
↭ times to the occurrence of a disease after exposing to a hazard material
↭ ages at the onset of a genetically related disease
↭ life times of machine components in industrial reliability
↭ times to complete specified tasks in psychological experimentation

Censoring: some subjects are not observed for the full time to failure.

Sources of (right) censoring:

↭ loss to follow-up
↭ alive at the end of the study
↭ death from other causes

5 / 46
Right-Censored Survival Data

Full data of the ith subject:

↭ failure time: Ti
↭ censoring time: Ci
↭ covariates: Xi = (Xi1 , . . . , Xip )T

Observed data of the ith subject:

↭ observation time: Yi = min(Ti , Ci )
↭ failure indicator: ωi = I(Ti → Ci )
↭ covariates: Xi = (Xi1 , . . . , Xip )T

6 / 46
Survival Analysis

Problems
↭ Estimating the failure time distribution
↭ Testing the equality of two failure time distributions
↭ Estimating the e!ects of covariates on the failure time

Naive methods
↭ Ignore censoring statuses: the failure times associated with censored
observations are underrepresented.
↭ Delete censored cases: estimates are biased towards smaller failure times
because larger failure times are more likely to be censored.

The need to accommodate censoring has been a key motivating factor for
the development of specialized statistical methods for analyzing failure
time data.

7 / 46
Example: Carcinogenesis

Table 1.1: Days to vaginal cancer mortality in rats

Group 1 143 164 188 188 190 192 206

209 213 216 220 227 230 234
246 265 304 216* 244*
Group 2 142 156 163 198 205 232 232
233 233 233 233 239 240 261
280 280 296 296 323 204* 344*
* Right-censored failure times.

In this example, the (right) censoring may have arisen because these four rats
died of causes unrelated to carcinogen, or they may not have died by the time
of data analysis.

8 / 46
Example: Randomized Clinical Trial
Table 1.2: Days to severe AGVHD in bone marrow transplant patients

9 / 46
Failure Time Distributions

The failure time T is a continuous nonnegative random variable. We are

interested in the following quantities:

Distribution function: F (t) = Pr(T → t)

Density function: f (t) = dF (t)/dt

Survival function: S(t) = Pr(T > t) = 1 ↑ F (t)

En F
Hazard function: instantaneous risk or rate of failure
1
ω(t) = lim Pr(t → T < t + !t | T ↓ t) = f (t)/S(t)
!t
!t→0

!t
Cumulative hazard function: ”(t) = 0 ω(u)du

10 / 46
Failure Time Distributions (Cont.)

We can easily verify the following relationships:

!↑
(i) S(t) = t
f (u)du

(ii) f (t) = ↑ dS(t)

dt

(iii) ω(t) = ↑ d logdtS(t)

(iv) S(t) = exp{↑”(t)}

(v) f (t) = ω(t) exp{↑”(t)}

11 / 46
Independent Censoring

Recall that Ti and Ci are the ith subject’s failure time and censoring time,
respectively, and Yi = min(Ti , Ci ) is the observation time.

A right-censoring mechanism is said to be independent if subjects

censored at each time t > 0 are “representative” of the subjects under
observation. In other words, subjects cannot be censored because they
appear to be at an unusually high or low risk of failure.

Usually, this means that Ti and Ci are statistically independent given the
covariates Xi .

Let F , G and H be the distribution functions of Ti , Ci and Yi ,

respectively. Under the independent censoring assumption, they satisfy
1 ↑ H = (1 ↑ F )(1 ↑ G).

12 / 46
Other Types of Censoring

Although the main focus of this course is right-censored failure time data, it’s
good to know other types of censoring.

Left censoring: a subject is observed to fail prior to some time t, but the
actual failure time is unknown.

Interval censoring: failure time is only known to fall within some interval (a, b).

Both left censoring and right censoring are special cases of interval censoring.
For right censoring, a is the observation time while b = ↔.
For left censoring, a = 0 while b is the observation time.

13 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

14 / 46
Notation

0 = t0 < t1 < · · · < tL = ε : partition of the study period [0, ε ]

dl : number of failures within the interval [tl↓1 , tl ), l = 1, . . . , L
rl : number of subjects at risk (i.e., under observation at the previous
instant) at tl↓1
T1↔ < T2↔ < · · · < TK↔ : distinct time points of observed failures
Dk : number of failures at Tk↔ , k = 1, . . . , K
Rk : number of subjects at risk at Tk↔ , k = 1, . . . , K

15 / 46
Nelson-Aalen Estimator
When !t is small enough, we have

”(t + !t) ↑ ”(t) ↗ ω(t)!t

↗ Pr(t → T < t + !t | T ↓ t)
= Pr(t → T < t + !t | T ↓ t, C ↓ t )
" #$ %
Y ↗t

Thus,

”(tl ) ↑ ”(tl↓1 ) ↗ Pr(tl↓1 → T < tl | Y ↓ tl↓1 ) ↗ dl /rl

A natural nonparametric estimator of ”(t) is

' '
”(t)
& = dl /rl ↘ Dk /Rk (Nelson-Aalen Estimator)
l:tl ↘t k:Tk→ ↘t

as L ↘ ↔ and max1↘l↘L |tl ↑ tl↓1 | ↘ 0.

16 / 46
Nelson-Aalen Estimator (Cont.)

Figure 1.1: Example of the Nelson-Aalen estimator

17 / 46
Kaplan-Meier Estimator
To estimate the survival function, note that

Pr(T > tl )
= Pr(T > t1 | T > t0 )Pr(T > t2 | T > t1 ) · · · Pr(T > tl | T > tl↓1 )

Thus,
(
S(t) ↗ Pr(T ↓ tl | T ↓ tl↓1 )
l:tl ↘t
(
= {1 ↑ Pr(T < tl | T ↓ tl↓1 )}
l:tl ↘t

A natural nonparametric estimator of S(t) is

( (
S(t)
& = (1 ↑ dl /rl ) ↘ (1 ↑ Dk /Rk ) (Kaplan-Meier Estimator)
l:tl ↘t k:Tk→ ↘t

as L ↘ ↔ and max1↘l↘L |tl ↑ tl↓1 | ↘ 0.

18 / 46
Kaplan-Meier Estimator (Cont.)

We can easily observe that S(T

& ↔ ) = S(T
k
& ↔ )(1 ↑ Dk /Rk ).
k↓1

For uncensored data, S

& reduces to one minus the empirical distribution
function.

An alternative estimator for S(t) is S(t)

) = e ↓&
”(t)
, which is close to S(t)
&
under certain condition.

19 / 46
Kaplan-Meier Estimator (Cont.)

Figure 1.2: Example of the Kaplan-Meier estimator

20 / 46
Exercise
Calculate the Kaplan-Meier estimator for the following data:

1 2 2 4↔ 5↔ 6 7↔ 8↔ 9↔ 10↔

21 / 46
Exercise
Calculate the Kaplan-Meier estimator for the following data:

1 2 2 4↔ 5↔ 6 7↔ 8↔ 9↔ 10↔

Dk
Tk↔ Dk Rk 1↑ Rk S(T
& ↔)
k
1 9 9
1 1 10 1↑ 10 = 10 10
2 7 9 7 7
2 2 9 1↑ 9 = 9 10 ≃ 9 = 10
1 4 7 4 14
6 1 5 1↑ 5 = 5 10 ≃ 5 = 25

21 / 46
Exercise
Calculate the Kaplan-Meier estimator for the following data:

1 2 2 4↔ 5↔ 6 7↔ 8↔ 9↔ 10↔

Dk
Tk↔ Dk Rk 1↑ Rk S(T
& ↔)
k
1 9 9
1 1 10 1↑ 10 = 10 10
2 7 9 7 7
2 2 9 1↑ 9 = 9 10 ≃ 9 = 10
1 4 7 4 14
6 1 5 1↑ 5 = 5 10 ≃ 5 = 25



1 if 0→t<1

9 if 1→t<2
S(t)
& = 10
7

 if 2→t<6
 10

14
25 if t↓6

21 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

22 / 46
Simple Graphical Comparison

Based solely on the above plots, we would conclude that treatment A yields
superior survival compared to treatment B. The hazard rate of group A is
initially lower than that of group B, but subsequently becomes higher than
that of group B.

23 / 46
Simple Graphical Comparison (Cont.)
What if we extend the follow-up period?

It turns out that at 15 time units, treatment B yields a much higher survival
probability compared to treatment A.

24 / 46
Simple Graphical Comparison (Cont.)
What if we extend the follow-up period?

It turns out that at 15 time units, treatment B yields a much higher survival
probability compared to treatment A.

The failure time distribution seems to di#er between the two treatment
groups, but how can we formally test for this di#erence?

24 / 46
Log-Rank Test
H0 : S1 (t) ⇐ S2 (t)

T1↔ < · · · < TK↔ : distinct time points of observed failures in combined
sample

Contingency table at Tk↔ (k = 1, . . . , K ):

Failures Non-failures At risk
Group 1 D1k R1k ↑ D1k R1k
Group 2 D2k R2k ↑ D2k R2k
Total Dk Rk ↑ Dk Rk

Under H0 , the conditional distribution of D1k given (R1k , R2k , Dk ) is

Hypergeometric(Rk , R1k , Dk ), with conditional mean and variance

Dk R1k Dk (Rk ↑ Dk )R1k R2k

E1k = , V1k =
Rk Rk2 (Rk ↑ 1)

25 / 46
Log-Rank Test (Cont.)
-

Thus, the statistic D1k ↑ E1k has conditional mean 0 and variance V1k .

Summing over the K failure times yields the log-rank statistic

.K
(D1k ↑ E1k )
Q = k=1
/. ⇒ N(0, 1)
K
k=1 V1k

For the carcinogenesis data, Q 2 = 3.12, p-value = 0.08.

Group 1 143 164 188 188 190 192 206
209 213 216 220 227 230 234
246 265 304 216↔ 244↔
Group 2 142 156 163 198 205 232 232
233 233 233 233 239 240 261
280 280 296 296 323 204↔ 344↔

26 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

27 / 46
Exponential Distribution

As before, T ↓ 0 is a random variable representing failure time, and t

represents a point in its range.

The one-parameter exponential distribution is obtained by taking the

um

hazard function to be constant, i.e., ω(t) ⇐ ω > 0, for all t ↓ 0.

i
um

Properties:
(i) Survival function: S(t) = e →ωt
(ii) Density function: f (t) = εe →ωt
(iii) Memoryless property: Pr(T ↑ t1 + t2 | T ↑ t1 ) = Pr(T ↑ t2 )
(iv) Expectation: E (T ) = 1/ε
(v) Moment-generating function: MT (s) = E (e sT ) = ω→s ω
for s < ε
(vi) If T ↓ Exp(ε), then for c > 0, cT ↓ Exp(ε/c).
(vii) If T1 and T2 are independently distributed according to Exp(ε1 ) and
Exp(ε2 ), respectively, then min(T1 , T2 ) ↓ Exp(ε1 + ε2 ).

28 / 46
Weibull Distribution

An important generalization of the exponential distribution allows for a

power dependence of the hazard on time. This yields the two-parameter
-

Weibull distribution with hazard function

ω(t) = ωp(ωt)p↓1 , for ω, p > 0.

Properties:
(i) Survival function: S(t) = exp{↔(εt)p }
(ii) Density function: f (t) = εp(εt)p→1 exp{↔(εt)p }
(iii) ε(t) is decreasing if p < 1, constant if p = 1, and increasing if p > 1.
(iv) If T ↓ Weibull(ε, p) and U = T p , then U ↓ Exp(εp ).

29 / 46
Weibull Distribution (Cont.)

Figure 1.3: Hazard functions for Weibull models with di!erent shape parameters

30 / 46
Gamma Distribution

Another -
two-parameter generalization of the exponential model is the
Gamma distribution with density function
-

ω(ωt)k↓1 e ↓ωt
f (t) = , for ω, k > 0.
$(k)

When k = 1, the Gamma distribution reduces to the exponential

distribution.

The moment-generating function is MT (s) = {ω/(ω ↑ s)}k .

Both the survival and hazard functions involve incomplete gamma

functions and cannot be written down explicitly. This restricts the use of
the Gamma distribution in practice.

31 / 46
Log-Normal Distribution

The failure time T has a log-normal distribution if and only if

log T ⇒ N(µ, ϑ 2 ).

The survival and hazard functions have no closed forms which involve
%(·), the distribution function of a standard normal random variable.

The log-normal model is particularly simple to apply if there is no

censoring, but with censoring the computations become more di&cult.

32 / 46
Log-Logistic Distribution

The failure time T has a log-logistic distribution if and only if

log T = µ + ϑZ , where Z follows the standard logistic distribution, with
density and distribution functions
ez ez
fZ (z) = , FZ (z) = , for ↑ ↔ < z < ↔.
(1 + e z )2 1 + ez

Let ω = e ↓µ and p = 1/ϑ, then the density, survival and hazard functions
of T are
ωp(ωt)p↓1 1 ωp(ωt)p↓1
f (t) = , S(t) = , ω(t) = .
{1 + (ωt)p }2 1 + (ωt)p 1 + (ωt)p

Like the exponential and Weibull models, this model has simple algebraic
expressions for the survival and hazard functions, thus is more convenient
than the log-normal distribution in handling censored data, while
providing a good approximation to it except in the extreme tails.

33 / 46
Gompertz Distribution

The relationship in which the hazard function is an exponential function of

the failure time T has been found to be descriptive in many investigations.

This leads to the Gompertz hazard ω(t) = ωe εt , for ω, ϖ > 0.

The survival and density functions of the Gompertz distribution are

 
ω
S(t) = exp (1 ↑ e εt ) ,
ϖ
 
ω
f (t) = ωe εt exp (1 ↑ e εt ) .
ϖ

34 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

35 / 46
Covariates

In the previous section, we consider modeling the survival distribution of a

homogeneous population.

However, the failure time T usually depends on some covariates X . For

example,
↭ COVID-19 vaccination status may a!ect the time of infection.
↭ Smoking may increase the risk of lung cancer.
↭ Some studies show that patients with Type 2 diabetes are at an increased
risk of dementia.

Regression models allow us to evaluate the associations between T and X .

36 / 46
Exponential Regression Models

The exponential distribution can be generalized to obtain a regression

model by allowing the hazard rate to be a function of the covariates:

ω(t; X ) = ω(X ).

Thus the hazard function for a given X is a constant characterizing an

exponential failure time distribution.

The function ω(·) can be parameterized in many ways, say

ω(t; X ) = ωg(ϖ T X ),

where ϖ = (ϖ1 , . . . , ϖp )T is a vector of regression parameters, ω is a

positive constant, and g(·) is a specified non-negative functional form.

37 / 46
Exponential Regression Models (Cont.)
The choice of g may depend on the particular data being considered.
Some examples are:
(i) g(u) = 1 + u, which can be interpreted as the hazard rate being a linear
function of X .
(ii) g(u) = (1 + u)→1 , which can be interpreted as the mean survival time
being a linear function of X .
(iii) g(u) = e u , which specifies that the log hazard rate is a linear function of
X.

Both (i) and (ii) su#er from the disadvantage that the set of ϖ values
considered must be restricted to guarantee that g(ϖ T X ) > 0.

In many ways, (iii) is the most natural form since it takes only positive
values. In this case, the hazard function of T given X is

ω(t; X ) = ω exp(ϖ T X ) (1)

38 / 46
Weibull Regression Models

The Weibull distribution can be generalized to the regression situation in

essentially the same way. For example, the hazard function can be
specified as
1
ω(t; X ) = ωp(ωt)p↓1 exp(ϖ T X ) (2)

Under both models (1) and (2), the e#ect of the covariates is to act
multiplicatively on the hazard function.

A more general regression model of this feature is the Cox model.

39 / 46
Cox Model
The Cox model specifies

ω(t; X ) = ω0 (t) exp(ϖ T X ), (3)

where ω0 (·) is an arbitrary unspecified baseline hazard function for

continuous T . Thus, this model is semiparametric.

It is easy to see that models (1) and (2) are both special cases of the Cox
model by setting ω0 (t) = ω and ω0 (t) = ωp(ωt)p↓1 , respectively.

The conditional survival and density functions of T given X are

  t 
T
S(t; X ) = exp ↑ exp(ϖ T X ) ω0 (u)du = {S0 (t)}exp(ε X ) ,
0
  t 
T T
f (t; X ) = ω0 (t) exp(ϖ X ) exp ↑ exp(ϖ X ) ω0 (u)du ,
0

where S0 (t) = S(t; X = 0) is the baseline survival function.

40 / 46
Accelerated Failure Time Model
To postulate a direct relationship between T and X , we consider another
semiparametric class of log-linear models for T :

log T = ϖ T X + W ,

where W is an error variable with unspecified density f .

Exponentiation of the above model gives T = exp(ϖ T X )S, where

S = exp(W ) > 0 has hazard function ω0 (s) that is independent of ϖ.

It follows that the survival and hazard functions for T can be written in
terms of this baseline hazard ω0 (·) according to
 T

S(t; X ) = exp ↑”0 (te ↓ε X ) ,
T
X
ω(t; X ) = ω0 (te ↓ε ) exp(↑ϖ T X ),
!t
where ”0 (t) = 0
ω0 (u)du.

41 / 46
Cox versus AFT Models

Figure 1.4: Baseline hazard function (X = 0) and hazard functions for X = 1 under
the Cox and AFT models

42 / 46
Table of Contents

1 Chapter 1: Overview
Introduction
Estimation of the Failure Time Distribution
Comparison of Survival Curves
Parametric Failure Time Models
Regression Models
Maximum Likelihood Estimation

43 / 46
Data and Likelihood
Observed data:
 
Yi = min(Ti , Ci ), ϱi = I(Ti → Ci ), Xi , i = 1, . . . , n

Model: any parametric (regression) model previously discussed

Unknown parameters: ς
Assumption: Ti ⇑ ⇑ Ci | Xi (independent censoring)
Likelihood contribution from the ith subject:
↭ f (Yi ; ϑ) if ωi = 1
↭ S(Yi ; ϑ) if ωi = 0
Likelihood:
n
(
Ln (ς) = f (Yi ; ς)ϑi S(Yi ; ς)1↓ϑi
i=1

Log-likelihood:
n 
' 
φn (ς) = ϱi log f (Yi ; ς) + (1 ↑ ϱi ) log S(Yi ; ς)
i=1

44 / 46
Maximum Likelihood Estimation

Let φ̇n (ς) and φ̈n (ς) be the gradient and Hessian of φn (ς) with respect to
ς, respectively.

The maximum likelihood estimator (MLE) for ς solves the score equation
φ̇n (ς) = 0.

If the solution to the score equation does not have an explicit expression,
we can approximate it using numerical methods such as the
Newton-Raphson algorithm, which iteratively updates the estimator ς&(k)
at the kth iteration by

ς&(k+1) = ς&(k) ↑ {φ̈n (ς&(k) )}↓1 φ̇n (ς&(k) )

until convergence.

45 / 46
Inference and Hypothesis Testing

Under some mild regularity conditions, the MLE ς&n satisfies

p
↭ consistency: ϑ
&n ↗ ϑ
↘ d
↭ asymptotic normality: n(ϑ&n ↔ ϑ) ↗ N(0, I(ϑ)→1 ), where I(ϑ) is the
Fisher information
↘ matrix. d
↭ Delta method: n{g(ϑ&n ) ↔ g(ϑ)} ↗ N(0, ġ(ϑ)T I(ϑ)→1 ġ(ϑ)), for any
continuous function g of ϑ.

Three commonly used tests for the null hypothesis H0 : ς = ς↔ (suppose ς

is an r -dimensional vector):
d
↭ Wald test: W = (ϑ
n &n ↔ ϑ↑ )T {nI(ϑ↑ )}(ϑ&n ↔ ϑ↑ ) ↗ ϖ2r
d
↭ Score test: SCn = ϱ̇n (ϑ↑ )T {nI(ϑ↑ )}→1 ϱ̇n (ϑ↑ ) ↗ ϖ2r
d
↭ Likelihood ratio test: LRCn = 2{ϱn (ϑ&n ) ↔ ϱn (ϑ↑ )} ↗ ϖ2r

Inference for semiparametric models such as the Cox and AFT models will
be studied in future chapters.

46 / 46