4.

Comparison of Survival Curves:

Hypotheses Testing

§4.1. Major procedures

Graphical comparison.
One-sample and Multipe-sample tests.
Tests for trend, Stratified tests and Test for
matched pairs.
More powerful tests for detecting crossing haz-
ard rate functions.
Tests at fixed time point.
§4.2. An example: Bone Marrow Trans-
plantation for Leukemia Patients.
In a cilinical trial, a total of 137 leukemia pa-
tients (99 with acute myeloctic leukemia, AML,
38 with acute lymphoblastic leukemia, ALL) went
1
bone marrow transplantation and then were fol-
lowed up until relapse, death or end of study.
The AML patients are further grouped as AML
low-risk (54 patients) and AML high-risk (45 pa-
tients). Except relapse, death, other events are
also observed. Information other cavariates of
the patients are also collected. For detail of the
trial, see Section 1.3 of Klein and Moeschberger
(2003). For the data and its description, see the
file [Link] on [Link]/stachenz~.
The following R-session produces a graphical
comparison of the survival functions of the three
groups of patients:
bmt=[Link]("[Link]")
[Link]=survfit(Surv(T2,d3)~g,
data=bmt)
plot([Link])

2
The estimated survival functions of
the three groups

3
§4.3. One-sample tests
• Data
Distinct event times: t1 < · · · < tD .
Number of events: d1, . . . , dD .
Number at risk: Y (t1), . . . , Y (tD).

• Hypothesis
H0 : h(t) = h0(t) for all t ≤ τ , where h0(t) is
known.

• Test statistic
Let
Z(τ ) = O(τ ) − E(τ )


D  τ
di
= W (ti) − W (s)h0(s)ds.
Y (ti) 0
i=1

4
The test statistic is defined as
Z(τ )
T (τ ) =  ,
V [Z(τ )]
where  τ
2 h0(s)
V [Z(τ )] = W (s) ds,
0 Y (s)
and W (·) is a weight function which reflects
different emphases over different parts of the
survival time.
Estimate of V [Z(τ )]:


D
h0(ti)
V

[Z(τ )] = 2
W (ti)
Y (ti)
i=1


D
di
or W 2(ti) 2
.
Y (ti)
i=1
Under H0,
T (τ ) → N (0, 1).
5
• The tests
One-sided: H1 : h(t) > h0(t).
Reject H0 at level α, if T (τ ) > z1−α .

Two-sided: H1 : h(t) = h0(t).

Reject H0 at level α, if |T (τ )| > z1−α/2.

• Weight function
Constant weight:
W (t) ≡ 1.
Log-rank weight:
W (t) = Y (t).
Harrington-Fleming weight:
p
W (t) = Y (t)S0 (t)[1 − S0(t)]q .
When W (t) is chosen as Y (t), the test is re-
ferred to as Log-rank test.
6
§4.4. Multi-sample tests
• Data:
K(≥ 2) independent censoring samples.
Pooled distinct event times: t1 < · · · < tD .
Number of events in jth sample: d1j , . . . , dDj .
Number at risk in jth sample: Y1j , . . . , YDj .
K K
di = j=1 dij , Yi = i=1 Yij .

• Hypotheses:
H0 : h1(t) = h2(t) = · · · = hK (t), for all
t ≤ τ,
H1 : at least one of the hj (t)’s is different from
others at some t ≤ τ .

7
• Construction of test statistic
Let W (t) be the weight function.
Measure of the difference of jth population
from the others:

D   
di
Zj (τ ) = W (ti) dij − Yij .
Yi
i=1
K
Note j=1 Zj (τ ) = 0.
Variances and covariances:
σ̂jj = V̂ [Zj (τ )]


D  
2 Yij Yij Yi − di
= W (ti) 1− di ,
Yi Yi Yi − 1
i=1
σ̂jg = Cov[Zj (τ ), Zg (τ )]


D  
2 Yij Yig Yi − di
= − W (ti) di .
Yi Yi Yi − 1
i=1

8
The test statistic:
ΣK−1 = (σ̂jg )j,g=1,...,K−1,

Z = (Z1(τ ), . . . , ZK−1(τ )) ,

2  −1
χ = Z ΣK−1Z.

Under H0, as sample sizes go to infinity,

χ2 → χ2K−1.

A heuristic argument on the statistic.

Consider the contingency table:
Sample
Failure j others Total
Yes dij di· di
No yij − dij Yi· − di· Yi − di
Total Yij Yi· Yi

9
Conditionaling on all the marginals, dij fol-
lows the hypergeometric distribution
k
N −k

dij n−dij
k(dij , n, N, k) = N
,
n
where n = Yij , N = Yi, k = di, with mean µd
and variance σd2 given by
nk di
µd = = Yij ,
N Yi
2 nk(N − k)(N − n)
σd = 2(N − 1)
N
 
Yij Yij Yi − di
= 1− di .
Yi Yi Yi − 1
The above conditional distribution motivates
the form of the test statistic.

10
• One-sided test for two samples
Z1(τ ) Z2(τ )
Z = √ =√
σ̂11 σ̂22
D di
i=1 W (t )[d
i i1 − Y i1 Yi )]
(
=  Y Y Y −d
.
D 2(t ) i1 [1 − i1 ] i i d
i=1 W i Y Y Y −1 i
i i i
Under H0, asymptotically,
Z ∼ N (0, 1).
For alternative H1 : h1(t) > h2(t), H0 is
rejected if Z > z1−α.
For alternative H1 : h1(t) < h2(t), H0 is
rejected if Z < zα.
For alternative H1 : h1(t) = h2(t), H0 is
rejected if |Z| > z1−α/2.

11
• The choice of W (t)
Let Ŝ(t) be the PL estimate obtained by pool-
ing the K samples together. Fleming and Har-
rington’s weight specifies:
Wp,q (ti) = Ŝ p(ti−1)[1− Ŝ(ti−1]q , p, q ≥ 0.
Other weights:
Wα (ti) = Yiα,

Log-rank : α=0
Tarone-Ware : α = 1/2,
Gehan : α = 1.
Define
 di
S̃(t) = (1 − ).
Yi + 1
ti ≤t
Peto-Peto: W (ti) = S̃(ti),
Modified Peto-Peto: W (ti) = S̃(ti)Yi/(Yi+ 1).
12
13
• R function for the tests
survdiff(formula, data, subset,
[Link], rho=0)
– formula: a formula expression as for other
survival models, of the form
Surv(time,status)~predictors.
For a one-sample test, the predictors must
consist of a single ’offset(sp)’ term, where
’sp’ is a vector giving the survival proba-
bility of each subject.
For a k-sample test, each unique combina-
tion of predictors defines a subgroup.
– rho: the exponent in the weight function:
W (t) = Ŝ ρ(t).

14
Example (cont.): Test the difference among
the three groups of Bone Marrow Transplant
patients:
bmt=[Link]("[Link]")
[Link]=survdiff(Surv(T2,d3)~g,
data=bmt)
[Link]
N Obs Exp (O-E)^2/E (O-E)^2/V
g=1 38 24 21.9 0.211 0.289
g=2 54 25 40.0 5.604 11.012
g=3 45 34 21.2 7.756 10.529

Chisq= 13.8 on 2 degrees of freedom,

p= 0.00101
[Link]$var
[,1] [,2] [,3]
[1,] 15.955175 -10.345092 -5.610084
[2,] -10.345092 20.339789 -9.994697
[3,] -5.610084 -9.994697 15.604781
15
§4.5. Tests for trend
• Hypotheses:
H0 : h1(t) = h2(t) = · · · = hK (t), for t ≤ τ ,
H1 : h1(t) ≤ h2(t) ≤ · · · ≤ hK (t), for t ≤ τ ,
with at least one strict inequality.
⇔ S1(t) ≥ · · · ≥ SK (t).

• Test statistic
Let
a1 < a2 < · · · < aK
be positive scores. In particular, aj = j.
The test statistic is given by
K
j=1 aj Zj (τ )
Z =   .
K K
j=1 g=1 aj ag σ̂jg

16
 Under H0, asymptotically,
Z ∼ N (0, 1).
H0 is rejected if Z ≥ z1−α .
Remark: The test is more powerful for detect-
ing the trend. However, it should only be used
when prior information indicates that the al-
ternatives are ordered.

• Computation issue
Use survdiff to obtain

Σ̂K = (σ̂jg ), Z = (Z1(τ ), . . . , ZK (τ )) .
Then compute

aZ
Z= .

a Σ̂a

a = (a1, . . . , aK ) .

17
Example: Bone Marrow Transplant (cont.)
Z=[Link]$[Link]$exp
AA=[Link]$var
a=c(1,2,3)
TT=t(a)%*%Z/sqrt(t(a)%*%AA%*%a)

Z: 2.148285 -14.966116 12.817830

AA:
15.955175 -10.345092 -5.610084
-10.345092 20.339789 -9.994697
-5.610084 -9.994697 15.604781

a: 1 2 3

TT: 1.631266

18
§4.6. Stratified tests
Assume the populations are stratified into M
strata by covariates.
s: index for stratum.
j, g: index for population.
i: index for individual.
• Construction of test statistic
For each stratum, construct the statistics


D   
dis
Zjs(τ ) = W (ti ) dijs − Yijs ,
i=1
Yis


D  
Y Y ijs Yis − dis
W 2(ti )
ijs
σ̂jjs = 1− dis ,
i=1
Y is Y is Y is − 1

D  
Y Y Y − d
W 2(ti )
ijs igs is is
σ̂jgs = − dis.
i=1
Y is Y is Y is − 1

Let


M

M
Zj· (τ ) = Zjs(τ ), σ̂jg· = σ̂jgs.
s=1 s=1

19
The test statistic is computed as

χ2 = (Z1· (τ ), · · · , ZK−1·(τ ))Σ−1
· (Z1· (τ ), · · · , ZK−1·(τ )) ,

where Σ· = (σ̂jg·)j,g=1,...,K−1.
Under H0, χ2 has an asymptotic chi square
distribution with d.f. K − 1.
Two-sample case: Test statistic is
M
Z1s(τ )
Z =  s=1 .
M
s=1 σ̂11s
Under H0, Z ∼ N (0, 1).
Both one- and two-sided alternatives can be
tested with Z.
The stratified tests can be computed by survdiff
with the specification of formula as follows:
Surv(time, status) ~ .+strata(s.v)
20
• Matched pair tests
By treating each matched pair as a stratum,
the stratified tests can be used to test matched
pairs. The two-sample statistic Z reduces to
D1 − D 2
Z= ,
D1 + D2
D1: the number of matched pairs in which
the individual from sample 1 expereience
the event first.
D2: the number of matched pairs in which
the individual from sample 2 expereience
the event first.
The matched pairs, where the smaller of the
two times is a censored observation, make no
contribution to the test statistic.

21
§4.7. More powerful tests for crossing
hazard rate fuctions
When hazard functions cross each other, the
tests in previous sections have little power. The
tests in this section are designed to test crossing
hazard functions with more power.

• Renyi Type tests

H0 : h1(t) = h2(t), t < τ ,
Notation:
t1 < · · · < tD : pooled distinct event times.
dij , Yij : # of events and at risk for group j.
di = di1 + di2, Yi = Yi1 + Yi2.
Weight function: W (ti) = 1 or Yi.

22
Construction of test statistic
For each ti, i = 1, . . . , D, compute

  
dk
Z(ti) = W (tk ) dk1 − Yk1 .
Yk
tk ≤ti
Let

Yk1Yk2(Yk − dk )
σ 2(τ ) = 2
W (tk )dk .
tk ≤τ
Yk (Yk − 1)
2

To test the two-sided alternative: h1(t) = h2(t),

the test statistic is given by
Q = max {|Z(ti)|}/σ(τ ).
1≤i≤D
Under H0,
Q ∼ sup{|B(x)|, 0 ≤ x ≤ 1},
where B(x) is the standard Brownian motion
process. Critical values are given in Table C.5
of Klein and Moeschberger (2003).
23
To test one-sided alternative: h1(t) > h2(t),
the test statistic is given by
Q∗ = max {Z(ti)}/σ(τ ).
1≤i≤D
Under H0,
Q∗ ∼ sup{B(x), 0 ≤ x ≤ 1}.
The p-value is determined as
P r(sup B(x) > Q∗) = 2[1 − Φ(Q∗)],
where Φ is the CDF of the standard normal
distribution.
A note on computation
The R function survdiff can be adapted
to obtain Z(ti) and σ(τ ).

24
§4.8. Tests at fixed point in time
It is needed sometimes to compare the survival
functions or cumulative incidence functions at a
fixed time point. The general method of multiple
comparison can be applied.

• A review of general multiple compar-

ison
A vector c = (c1, . . . , cK ) is said to be a con-
trast vector if

K
cj = 0.
j=1
Let C be a q × K matrix consists of q contrast
row vectors,⎛i.e., ⎞

c1
⎜ . ⎟
C = ⎝ . ⎠.

cK
25
 
To test H0 : Cθ = 0, where θ = (θ1, . . . , θK ) ,
the test statistic is formed as

  −1
2
χ = θ̂ C [C Σ̂θ C ] C θ̂,
where θ̂, which is asymptotically normally dis-
tributed, is a consistent estimator of θ, and Σ̂θ
is the estimated variance matrix of θ̂.
Under H0, χ2 has an asymptotic χ2-distribution
with d.f. q.
The hypothesis: θ1 = · · · = θK is equivalent
to Cθ = 0 with
⎛ ⎞
1 0 ··· 0 −1
⎜ 0 1 ··· 0 −1 ⎟
C=⎜
⎝··· ···
⎟.
··· ··· ··· ⎠
0 0 ··· 1 −1

26
• Tests for survival functions
Suppose there are K groups with survival func-
tions S1(t), . . . , SK (t). It is to test whether
the survival functions differ at a particular point
t0 .

θ = (S1(t0), . . . , SK (t0))

θ̂ = (Ŝ1(t0), . . . , ŜK (t0))

Σ̂θ = Diag(V̂ [Ŝ1(t0)], . . . , V̂ [ŜK (t0)]).

Ŝj (t0) : PL estimate of Sj (t0).

27
• Tests for cumulative incidence func-
tions
The purpose of the tests is to see whether the
incidence functions of a certain risk among
K different groups differ at a particular time
point t0. Suppose the risk factor of concern is
indexed by r.

θ = (F1r (t0), . . . , FKr (t0))

θ̂ = (F̂1r (t0), . . . , F̂Kr (t0))

Σ̂θ = Diag(V̂ [F̂1r (t0)], . . . , V̂ [F̂Kr (t0)]).

where Fjr (t0) is the rth cumulative incedence

function in the jth group.

28
• Pairwise comparison with Bonforonni
adjustment
Test Hij0 : θi = θj for M pairs.
Test statistics:
θ̂i − θ̂j
Zij =  .
V̂ [θ̂i − θ̂j ]
Simultaneous critical value:
cα = z1−α/M .
Reject Hij0, if |Zij | > cα.

• Computation
Use survfit and [Link] to obtain θ̂, Σ̂θ
respectively for testing the survival functions
and cumulative incidence functions.

29
• Example: Bone Marrow Transplant (cont.)
The estimated relapse cumulative incidence
functions at t0 = 365 days for ALL, AML
low risk and AML high risk groups and their
corrsponding variances are obtained using
cumu_inci as
F̂ALL = 0.238, V̂ (F̂ALL) = 0.0048;
F̂AMLL = 0.093, V̂ (F̂ALML) = 0.0015;
F̂AMLH = 0.378, V̂ (F̂AMLH) = 0.0056.
Contrast matrix
 
1 0 −1
C= .
0 1 −1
 
 0.01037 0.00559
C Σ̂C = .
0.00559 0.00713

χ2 = 12.59, p-value = 0.0018

30
R-session for the computation:
bmt=[Link]("[Link]")
[Link] = [Link](time=bmt$T2,
status=bmt$d2+bmt$d3,g=bmt$g)
g=[Link]$g

inci1=[Link](K=2,data=[Link],
subset=(g==1),plot=T)
inci2=[Link](K=2,data=[Link],
subset=(g==2),plot=T)
inci3=[Link](K=2,data=[Link],
subset=(g==3),plot=T)

e.time1 = inci1$[Link]
e.time2 = inci2$[Link]
e.time3 = inci3$[Link]
n1 = length(e.time1)
n2 = length(e.time2)
n3 = length(e.time3)
31
t1=max(c(1:n1)[e.time1<=365])
t2=max(c(1:n2)[e.time2<=365])
t3=max(c(1:n3)[e.time3<=365])

snv=rbind(inci1[t1,c(3,5)],
inci2[t2,c(3,5)],
inci3[t3,c(3,5)])

C=matrix(c(1,0,-1,0,1,-1),ncol=3,
byrow=T)
chi=t(C%*%snv[,1])%*%
solve(C%*%diag(snv[,2])%*%t(C))%*%
C%*%snv[,1]
1-pchisq(chi,2)

32