Survival Analysis - Universität Rostock

Survival Analysis
(Modul: Lebensdaueranalyse)
R OLAND R AU
Universität Rostock, Sommersemester 2015
16. Juni 2015
© Roland Rau
Survival Analysis
1 / 18
Teams, Themen, & Termine
Zwei Termine mit jeweils 3 “slots” stehen zur Auswahl:
07. Juli & 14. Juli
#
1
2
3
4
5
6
Thema
Versicherungsstatus & Sterblichkeit
SES & Sterblichkeit
BMI & Sterblichkeit
Rauchersterblichkeit
Kosten für med. Hilfe & Sterblichkeit
Self rated health & Sterblichkeit
© Roland Rau
Team
Grützmacher, Smilde-Becker
Lepzien, Schulz
Köpke, Liewert, Sauerberg
Dronia, Eckardt
Epperlein
Wenau
Survival Analysis
Termin / “Slot”
1/2
2/1
1/1
2/3
1/3
2/2
2 / 18
Proportional Hazards Modelle
Grundsätzliche Form:
hi (t|X) = h0 (t)eβ1 x1,i +β2 x2,i +...
Trennung in
Baseline Hazard (h0 (t)) und
Effekt/e (β1 , β2 , . . .) der Kovariaten (x1 , x2 , . . .).
Falls h0 (t) parametrisch spezifiziert ⇒ parametrische
Survival-Regression
Falls h0 (t) ⇒ Cox-Modell
© Roland Rau
Survival Analysis
3 / 18
Zur Erinnerung I
© Roland Rau
Survival Analysis
4 / 18
Zur Erinnerung II
wolpertinger <- read.table("wolpertinger-geschlecht-2015.txt",
header=TRUE)
geschlecht <- ifelse(wolpertinger$sex==1, 1, 0)
loglike <- function(theta, t, covariate,delta) {
lambda <- theta[1]
beta <- theta[2]
return(-(sum( (log(lambda) + beta*covariate)*delta) +
sum(-lambda * exp(beta * covariate) * t)))
}
ergebnis <- optim(par=c(0.5, 1.2), fn=loglike, t=wolpertinger$zeit,
delta=wolpertinger$cens, covariate=geschlecht)
##
##
##
##
Warning:
Warning:
Warning:
Warning:
NaNs
NaNs
NaNs
NaNs
produced
produced
produced
produced
ergebnis
##
##
##
##
##
##
##
##
##
##
##
##
##
$par
[1] 0.1584 0.2470
$value
[1] 929.5
$counts
function gradient
71
NA
$convergence
[1] 0
© Roland Rau
Survival Analysis
5 / 18
Einfacher: I
library(survival)
s.obj <- Surv(wolpertinger$zeit, wolpertinger$cens)
survreg(s.obj ~ geschlecht, dist="exponential")
##
##
##
##
##
##
##
##
##
##
##
##
Call:
survreg(formula = s.obj ~ geschlecht, dist = "exponential")
Coefficients:
(Intercept) geschlecht
1.8423
-0.2469
Scale fixed at 1
Loglik(model)= -929.5
Loglik(intercept only)= -932.1
Chisq= 5.22 on 1 degrees of freedom, p= 0.022
n= 1000
© Roland Rau
Survival Analysis
6 / 18
Zur Erinnerung: Cox-Modell I
coxmodel <- coxph(s.obj ~ geschlecht)
coxmodel
##
##
##
##
##
##
##
##
Call:
coxph(formula = s.obj ~ geschlecht)
coef exp(coef) se(coef)
z
p
geschlecht 0.249
1.28
0.109 2.3 0.022
Likelihood ratio test=5.29
on 1 df, p=0.0214
© Roland Rau
n= 1000, number of events= 343
Survival Analysis
7 / 18
Numerisches Beispiel: Wolpertinger (reloaded)
Wie interpretieren wir nun diese Ergebnisse?
round(coef(coxmodel),4)
## geschlecht
##
0.2492
Da wir die Variable x so gewählt haben, dass für das weibliche Geschlecht x = 0 gilt, so erkennen wir:
βx
h(t) = h0 e
β0
= h0 e
0
= h0 e = h0
dass der Baseline-Hazard λ =0.247 dem Hazard der Frauen entspricht. Dieser wird bei Männern um den Faktor
βx
h(t) = h0 e
β1
= h0 e
0.2492
= h0 e
= h0 × 1.282999
erhöht.
Das Hazard Ratio beträgt damit:
HR =
β x
h0 (t)e 1 ♂
h0
β x
(t)e 1 ♀
=
β x
h0 (t)e 1 ♂
h0 (t)
β x
= e 1 ♂ = 1.282999
Das bedeutet, dass in unserem Modell Männer in jeder Altersstufe ein in etwa 30% höheres Sterberisiko als Frauen
besitzen (28.2999%).
© Roland Rau
Survival Analysis
8 / 18
Das Cox-Modell mit unseren Daten I
Zuerst müssen wir natürlich wieder unsere Daten einlesen und
leicht überarbeiten:
newdata <- read.table("smalldataset2015.txt", header=TRUE, sep=",")
newdata2 <- subset(newdata, status %in% c(0,1))
newdata3 <- subset(newdata2, age %in% 18:84)
newdata4 <- subset(newdata3, yob %in% 1914:1980)
eintrittsalter <- newdata4$age
delta <- newdata4$status
austrittsalter <- ifelse(newdata4$status == 0,
2007 - newdata4$yob, # also falls rechtszensiert
newdata4$year - newdata4$yob+0.5) # falls Ereignis
surv.obj <- Surv(eintrittsalter, austrittsalter, delta)
© Roland Rau
Survival Analysis
9 / 18
Ein erstes Cox-Modell mit “unseren” Daten
Wir berechnen nun mit unseren Daten ein Modell analog zu den
Wolpertingern:
coxmodel1 <- coxph(surv.obj ~ as.factor(newdata4$sex))
coxmodel1
##
##
##
##
##
##
##
##
Call:
coxph(formula = surv.obj ~ as.factor(newdata4$sex))
coef exp(coef) se(coef)
z p
as.factor(newdata4$sex)2 -0.436
0.646
0.0356 -12.3 0
Likelihood ratio test=149
on 1 df, p=0
n= 33364, number of events= 3234
Wir sehen also, dass das Geschlecht 2 (=Frauen) ein niedrigeres
Sterberisiko hat als die Referenzkategorie Männer (β = −0.4362).
Dies können wir als relatives Risiko interpretieren in der Form von
β = e−0.4362 = 0.6465
Wenn die Annahmen unseres Modells also stimmen, so haben
Frauen in jeder Altersstufe ein 1 − 0.6465=0.3535 =35.35%
niedrigeres Sterberisiko.
© Roland Rau
Survival Analysis
10 / 18
Ein zweites Cox-Modell mit “unseren” Daten
Welchen Einfluss hat das Rauchen auf die Sterblichkeit? Wir operationalisieren diese
Frage mittels der Variable eversmoke.
table(newdata4$eversmoke)
##
##
1
2
## 15898 17325
7
86
8
8
9
47
Wir können uns nun entscheiden, ob wir die Kategorien 7–9 ausschliessen oder aber
als Residualkategorie zusammenfassen.
newdata4$raucher <- "Residual"
newdata4$raucher[newdata4$eversmoke==1] <- "Raucher"
newdata4$raucher[newdata4$eversmoke==2] <- "Nichtraucher"
table(newdata4$raucher)
##
## Nichtraucher
##
17325
Raucher
15898
Residual
141
Es bietet sich zudem an, die Variable explizit als kategoriale Variable zu definieren mit
einer Referenzgruppe.
newdata4$raucher <- as.factor(newdata4$raucher)
newdata4$raucher <- relevel(newdata4$raucher, ref="Nichtraucher")
© Roland Rau
Survival Analysis
11 / 18
Ein zweites Cox-Modell mit “unseren” Daten
coxmodel2 <- coxph(surv.obj ~ as.factor(newdata4$sex) + newdata4$raucher)
coxmodel2
##
##
##
##
##
##
##
##
##
##
Call:
coxph(formula = surv.obj ~ as.factor(newdata4$sex) + newdata4$raucher)
coef exp(coef) se(coef)
z
p
as.factor(newdata4$sex)2 -0.293
0.746
0.0366 -8.012 1.1e-15
newdata4$raucherRaucher
0.571
1.769
0.0382 14.932 0.0e+00
newdata4$raucherResidual -0.023
0.977
0.3030 -0.076 9.4e-01
Likelihood ratio test=380
on 3 df, p=0
© Roland Rau
n= 33364, number of events= 3234
Survival Analysis
12 / 18
Kleine Aufgabe
Wie hoch ist mit den gegebenen Schätzwerten der vorherigen Folie das relative
Sterberisiko
von rauchenden Frauen im Vergleich zu nichtrauchenden Männern?
von rauchenden Männern im Vergleich zu nichtrauchenden Frauen?
HR =
h(t)Frau, raucht
ho (t)eβ1 x1,Frau +β2 x2,raucht
1.3205
e−0.293·1+0.571·1
=
≈
≈ 1.3205
=
β
x
+β
x
1
1,Mann
2
2,rauchtnicht
h(t)Mann, raucht nicht
e0·0+0·0
1
ho (t)e
HR =
h(t)Mann, raucht
ho (t)eβ1 x1,Mann +β2 x2,raucht
e0·0+0.571·1
1.77004
= −0.293·1+0·0 ≈
=
≈ 2.3726
β
x
+β
x
h(t)Frau, raucht nicht
e
0.74602
ho (t)e 1 1,Frau 2 2,rauchtnicht
© Roland Rau
Survival Analysis
13 / 18
Wie können wir graphisch diese Ergebnisse
darstellen?
© Roland Rau
Survival Analysis
14 / 18
Nächste Veranstaltung:
Das Cox-Modell mit unseren Daten
Wie überprüfe ich die Proportionalitätsannahme?
Was tun, wenn die Proportionalitätsannahme verletzt ist?
Wie plotte ich Ergebnisse in R?
© Roland Rau
Survival Analysis
15 / 18
Vielen Dank für Ihre Aufmerksamkeit!
© Roland Rau
Survival Analysis
16 / 18
Lizenz
This open-access work is published under the terms of the Creative
Commons Attribution NonCommercial License 2.0 Germany, which
permits use, reproduction & distribution in any medium for non-commercial
purposes, provided the original author(s) and source are given credit.
Für ausführlichere Informationen:
http://creativecommons.org/licenses/by-nc/2.0/de/ (Deutsch)
http://creativecommons.org/licenses/by-nc/2.0/de/deed.en (English)
© Roland Rau
Survival Analysis
17 / 18
Kontakt
Universität Rostock
Institut für Soziologie und Demographie
Lehrstuhl für Demographie
Ulmenstr. 69
18057 Rostock
Germany
Tel.: +49-381-498 4044
Fax.: +49-381-498 4395
Email: [email protected]
Sprechstunde im Sommersemester 2015: Mittwochs, 09:00–10:00
(und nach Vereinbarung)
© Roland Rau
Survival Analysis
18 / 18