Übung zur Survival-Analysis

Survival bei konstantem Hazard

Wir nehmen eine über die Zeit konstante Hazard-Funktion an: \[h(t)=h=4\]

Berechnen Sie die Dezile der Survival-Zeit. Dazu brauchen Sie zuerst die Survival-Funktion bei konstanter Hazard-Rate. Diese lösen Sie dann nach \(t\) auf und berechnen die gewünschten Quantile (mit den entsprechenden Werten auf \(S(t)\)).
Welcher Verteilung \(F(t)=\Pr(T\leq t)\) entspricht konstanter Hazard? Lösen Sie das Problem zusätzlich mit der entsprechenden q Funktion in R (analog zu qnorm(), qt(), qchisq() usw.).
Berechne die Rate bei bekannter Halbwertszeit t=0.43 (bei Annahme von konstanter Rate.)

Die Survival-Funktion ist \[ S(t)=\exp(-H(t))=\exp\left(-\int_0^t h(u)du\right). \]

Wenn \(h(t)=h\) konstant, dann vereinfacht sich dies zu \[ S(t)=\exp(-ht) \]

und bei Auflösen nach \(t\): \[ t=-\log(S(t))/h \]

p <- seq(1,0,by=-0.1)##probs of survival
h <- 4 #constant hazard
t<- -log(p)/h ##quantiles
data.frame(Survival=p,Quantil=t)
##    Survival    Quantil
## 1       1.0 0.00000000
## 2       0.9 0.02634013
## 3       0.8 0.05578589
## 4       0.7 0.08916874
## 5       0.6 0.12770641
## 6       0.5 0.17328680
## 7       0.4 0.22907268
## 8       0.3 0.30099320
## 9       0.2 0.40235948
## 10      0.1 0.57564627
## 11      0.0        Inf

h<-4
time<-seq(0,2,by=0.1)
S<-exp(-h*time)
plot(time, S, ylab = "survival", type = "l",ylim=c(0,1))
abline(h=p,v=-log(p)/h,lty=2)
abline(h=0.5,v=-log(0.5)/h,col="red",lty=1)

Verteilung: \(F(t)=1-\exp(-ht)\) ist die Exponentialverteilung und \(f(t)=h\exp(-ht)\) ist die Exponentialdichte. Wir könnten also auch qexp(p=1-p,rate=4) brauchen. Das erste Argument ist 1-p, weil wir Quantile der Survival und nicht der Failure wollen.

qexp(p=1-p,rate=4)
##  [1] 0.00000000 0.02634013 0.05578589 0.08916874 0.12770641 0.17328680
##  [7] 0.22907268 0.30099320 0.40235948 0.57564627        Inf

Schauen wir uns die Verteilung noch graphisch an, mit dexp() und pexp():

curve(dexp(x,rate=4),from=0,to=2,xlab="Time",ylab="Density")
curve(pexp(x,rate=4),from=0,to=2,xlab="Time",ylab="Cumulative Failure")
curve(1-pexp(x,rate=4),from=0,to=2,xlab="Time",ylab="Survival")

Für die Interessierten: Konstante Hazard entsprechen einer gedächtnislosen Verteilung (wie beim radioaktiven Zerfall). Formal bedeutet Gedächtnislosigkeit, dass

\[ \Pr(T>s+t|T>s)=\Pr(T>t), \quad s,t\geq 0 \]

Die Exponentialfunktion ist gedächtnislos: \[ \Pr\left(T > s + t \mid T > s\right) = \frac{\Pr\left(T > s + t \cap T > s\right)}{\Pr\left(T > s\right)} = \frac{\Pr\left(T > s + t \right)}{\Pr\left(T > s\right)} = \frac{e^{-\lambda(s + t)}}{e^{-\lambda s}}= e^{-\lambda t} \\[4pt] = \Pr(T > t). \]

Nach \(h\) auflösen gibt: \(h=-\log(0.5)/t\)

Die Rate ist dann

-log(0.5)/0.43
## [1] 1.61197

Kaplan-Meier

Problem: Überlebensrate bei Patienten mit akuter myeloischer Leukämie (AML). Dabei stellte sich die Frage, ob der Standardverlauf der Chemotherapie um zusätzliche Zyklen erweitert werden sollte (sogenannte “Erhaltungstherapie”), siehe ?aml

library(survival)
aml
##    time status             x
## 1     9      1    Maintained
## 2    13      1    Maintained
## 3    13      0    Maintained
## 4    18      1    Maintained
## 5    23      1    Maintained
## 6    28      0    Maintained
## 7    31      1    Maintained
## 8    34      1    Maintained
## 9    45      0    Maintained
## 10   48      1    Maintained
## 11  161      0    Maintained
## 12    5      1 Nonmaintained
## 13    5      1 Nonmaintained
## 14    8      1 Nonmaintained
## 15    8      1 Nonmaintained
## 16   12      1 Nonmaintained
## 17   16      0 Nonmaintained
## 18   23      1 Nonmaintained
## 19   27      1 Nonmaintained
## 20   30      1 Nonmaintained
## 21   33      1 Nonmaintained
## 22   43      1 Nonmaintained
## 23   45      1 Nonmaintained

Berechnen Sie einen Kaplan-Meier Schätzer für die Survival-Kurve:

overall
für beide Gruppen
Plotten sie den Kaplan-Meier Schätzer

Wie gross sind die Median-Survival für beide Gruppen?
Kann die Nullhypothese, dass sich die beiden Kurven nicht voneinander unterscheiden, verworfen werden? Machen Sie einen Log-Rank Test
Reproduzieren Sie den Zwischengruppentest mit einem Cox-Modell, mit summary(coxph(Surv(time, status)~x,aml)). Um welchen Faktor ist der Hazard bei x=Nonmaintained grösser als bei x=Maintained? Kann der Nulleffekt verworfen werden?

Survival estimation mit Kaplan-Meier

overall

m.KM1 <- survfit(Surv(time, status) ~ 1, data = aml)
plot(m.KM1,,xlab="time",ylab="survival")

by group

m.KM2 <- survfit(Surv(time, status) ~ x, data = aml)
plot(m.KM2,xlab="time",ylab="survival",lty=2:3)
legend(50, 0.9, c("Maintained", "Nonmaintained"), lty = 2:3)

Median-Survival by group

m.KM1
## Call: survfit(formula = Surv(time, status) ~ 1, data = aml)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 23     18     27      18      45
m.KM2
## Call: survfit(formula = Surv(time, status) ~ x, data = aml)
## 
##                  n events median 0.95LCL 0.95UCL
## x=Maintained    11      7     31      18      NA
## x=Nonmaintained 12     11     23       8      NA

Log-Rank Test

survdiff(Surv(time,status)~x,data=aml)
## Call:
## survdiff(formula = Surv(time, status) ~ x, data = aml)
## 
##                  N Observed Expected (O-E)^2/E (O-E)^2/V
## x=Maintained    11        7    10.69      1.27       3.4
## x=Nonmaintained 12       11     7.31      1.86       3.4
## 
##  Chisq= 3.4  on 1 degrees of freedom, p= 0.07

Die Nullhypothese kann nicht verworfen werden.

survminer::ggsurvplot(m.KM2, data = aml, conf.int = TRUE, surv.median.line = "hv", legend.labs = c("Maintained", "Nonmaintained"), pval = TRUE)

Cox-Modell

summary(coxph(Surv(time, status)~x,aml))
## Call:
## coxph(formula = Surv(time, status) ~ x, data = aml)
## 
##   n= 23, number of events= 18 
## 
##                  coef exp(coef) se(coef)     z Pr(>|z|)  
## xNonmaintained 0.9155    2.4981   0.5119 1.788   0.0737 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##                exp(coef) exp(-coef) lower .95 upper .95
## xNonmaintained     2.498     0.4003    0.9159     6.813
## 
## Concordance= 0.619  (se = 0.063 )
## Likelihood ratio test= 3.38  on 1 df,   p=0.07
## Wald test            = 3.2  on 1 df,   p=0.07
## Score (logrank) test = 3.42  on 1 df,   p=0.06

Die Hazard Ratio ist 2.5. Ohne Erhaltungstherapie ist die geschätzte Hazard 2.5 mal grösser als mit Erhaltungstherapie. Der Effekt könnte aber – mit Risiko 0.05 – auch knapp in die andere Richtung gehen (untere Grenze des 95% CI ist 0.916). Der Nulleffekt kann auf \(0.05\)-Niveau nicht verworfen werden.

Sterbetafel

Betrachten Sie folgende Sterbetafel:

Table 1: Life Table
Month	AliveBegin	Deaths	Censored	AtRisk	RiskDeath	RiskSurv	Survival
1	240	12	0	240.0	0.0500000	0.9500000	0.9500000
2	228	9	0	228.0	0.0394737	0.9605263	0.9125000
3	219	17	1	218.5	0.0778032	0.9221968	0.8415046
4	201	36	4	199.0	0.1809045	0.8190955	0.6892726
5	161	6	2	160.0	0.0375000	0.9625000	0.6634249
6	153	18	7	149.5	0.1204013	0.8795987	0.5835476
7	128	13	5	125.5	0.1035857	0.8964143	0.5231005
8	110	11	3	108.5	0.1013825	0.8986175	0.4700672
9	96	14	3	94.5	0.1481481	0.8518519	0.4004276
10	79	13	0	79.0	0.1645570	0.8354430	0.3345345
11	66	15	4	64.0	0.2343750	0.7656250	0.2561280
12	47	6	1	46.5	0.1290323	0.8709677	0.2230792
13	40	6	0	40.0	0.1500000	0.8500000	0.1896173
14	34	4	2	33.0	0.1212121	0.8787879	0.1666334
15	28	5	0	28.0	0.1785714	0.8214286	0.1368774
16	23	7	1	22.5	0.3111111	0.6888889	0.0942933
17	15	12	0	15.0	0.8000000	0.2000000	0.0188587
18	3	3	0	3.0	1.0000000	0.0000000	0.0000000

Schreiben Sie einen (halben) Einzeiler in R, wie man zum geschätzten Survival Wert im Monat fünf kommt, also \(\hat{S}(5)=0.6634249\).

(1-12/240)*(1-9/228)*(1-17/218.5)*(1-36/199)*(1-6/160)
## [1] 0.6634249

oder (bisschen weniger schön…)

0.95*0.9605*0.9222*0.8191*0.9625
## [1] 0.6634127

Allgemein:

\[ \boxed{\hat{S}(t) = \prod_{i=1}^t (1 - \hat{r}_i) = (1 - \hat{r}_1) \times \dots \times (1 - \hat{r}_t)} \]

mit

\[ \boxed{\hat{r}_i = \frac{\text{Anzahl der Patienten, die im Intervall } i \text{ gestorben sind}}{\text{Anzahl der Patienten, die zu Beginn des Intervalls } i \text{ noch lebten}} = \frac{d_i}{n_i}.} \]