Survival Analysis (Modul: Lebensdaueranalyse) · Nichtparametrische Schätzung: “Kaplan-Meier”...

Survival Analysis(Modul: Lebensdaueranalyse)

ROLAND RAU

Universität Rostock, Sommersemester 2014

20. Mai 2014

© Roland Rau Survival Analysis 1 / 25

“Teams & Themen”


Kovariaten

Einführung zu:

Welchen Einfluss üben Kovariatenauf die Prozesszeit X aus?

Heutige Veranstaltung und kommende Veranstaltungen:

nicht-parametrische Schätzungparametrische Schätzungsemi-parametrische Schätzung

+ kleinere Arbeiten mit dem Datensatz.


Nichtparametrische Schätzung: “Kaplan-Meier”

Ziel des Kaplan-Meier Schätzers:(nichtparametrische) Schätzung der Survival-Funktion:

auch wenn die Beobachtungen nicht vollständig sind (Derursprüngliche “KM” hat lediglich rechts-zensierteBeobachtungen berücksichtigt. Eine Erweiterung umLinks-Trunkierung ist jedoch leicht möglich.)und man zusätzlich den Einfluss von Kovariaten auf denSurvival-Verlauf berücksichtigen möchte.



Der Kaplan-Meier Schätzer ist folgendermaßen definiert für dieProzesszeit x und die Zeitpunkte i, an denen sich mindestensein “Event” ereignet (siehe für Definition Abschnitt 4.2 in Kleinand Moeschberger (2003)1).

S (x) =

{1, falls x < x1∏

xi≤x

(1 − di

Yi

)falls x1 ≤ x

wobei di der Anzahl der Ereignisse (z.B. Todesfälle) zum i-tenZeitpunkt entspricht und Yi der Anzahl der Personen, die umi-ten Zeitpunkt dem Risiko ausgesetzt waren, das Ereignis zuerfahren.

1Es sei darauf hingewiesen, dass Klein and Moeschberger (2003) nicht xsondern t verwenden.



Beispiel 1: komplett beobachtete Überlebenszeiten

7 3 24 19 8 7



Beispiel 1: komplett beobachtete Überlebenszeiten

7 3 24 19 8 7

xi di Yi S(x)3 1 6 0.833337 2 5 0.500008 1 3 0.33333

19 1 2 0.1666724 1 1 0.00000

0 5 10 15 20

0.0

0.2

0.4

0.6

0.8

1.0

x

S(x)



Beispiel 2: komplett beobachtete und rechts-zensierte (“+”)Überlebenszeiten

7 3 24 19 8 7 11+ 19+



Beispiel 2: komplett beobachtete und rechts-zensierte (“+”)Überlebenszeiten

7 3 24 19 8 7 11+ 19+

xi di Yi S(x)3 1 8 0.875007 2 7 0.625008 1 5 0.50000

19 1 3 0.3333324 1 1 0.00000

0 5 10 15 20

0.0

0.2

0.4

0.6

0.8

1.0

x

S(x)



Zur Berechnung von Konfidenz-Intervallen wird typischerweisedie “Greenwood-Formel” für die geschätzte Varianzherangezogen:

V(

S (x))= S (x)2

∑xi≤x

di

Yi (Yi − di)

. . . aber das lassen wir den Computer machen (allerdingsnimmt der eine andere Formel, da es eine Vielzahl anSchätzern für Konfidenzintervalle für denKaplan-Meier-Schätzer gibt).



R ist unser Freund und ist schon hervorragend für die Survival-Analyseaufgestellt.Hierzu benötigen wir das Paket survival, das zu den Standardpaketen vonR gehört.

library(survival)zeiten <- c(7, 3, 24, 19, 8, 7, 11, 19)delta <- c(1, 1, 1, 1, 1, 1, 0, 0)

Zuerst muss man bei der Survival-Analyse R erst einmal mitteilen, was dieProzesszeit ist und welche Variable anzeigt, ob es sich um einenrechtszensierten (δ = 0) oder komplett beobachteten (δ = 1) Fall handelt.

survivalobjekt <- Surv(zeiten, delta) ## Surv, nicht surv (!)


Nichtparametrische Schätzung: “Kaplan-Meier”Den Kaplan-Meier-Schätzer berechnet man folgendermaßen:

km1 <- survfit(survivalobjekt ~ 1)

Die “1” rechts der Tilde bedeutet, dass es keinerlei Kovariaten gibt . . . dazu kommen wir gleich.Die wichtigsten Information erhalten wir ganz einfach mittels:

km1

## Call: survfit(formula = survivalobjekt ~ 1)#### records n.max n.start events median 0.95LCL 0.95UCL## 8.0 8.0 8.0 6.0 13.5 7.0 NA

Etwas ausführlicher gibt es das mittels:

summary(km1)

## Call: survfit(formula = survivalobjekt ~ 1)#### time n.risk n.event survival std.err lower 95% CI upper 95% CI## 3 8 1 0.875 0.117 0.673 1.000## 7 7 2 0.625 0.171 0.365 1.000## 8 5 1 0.500 0.177 0.250 1.000## 19 3 1 0.333 0.180 0.116 0.961## 24 1 1 0.000 NaN NA NA


Nichtparametrische Schätzung: “Kaplan-Meier”Einer der vielen Vorteile von R liegt auch in den grafischen Fähigkeiten:

plot(x=km1, col=c("black", "red", "red"), xlab="x",

ylab="S(x)", lty=1)

0 5 10 15 20

0.0

0.2

0.4

0.6

0.8

1.0

x

S(x

)


Nichtparametrische Schätzung: “Kaplan-Meier”Das Niveau des Konfidenzintervalls lässt sich durch das Argument conf.int bei derSchätzung mittels survfit steuern. Standardmäßig wird ein 95% Niveau geschätzt.

km1.ci90 <- survfit(survivalobjekt ~ 1, conf.int=0.9)plot(x=km1.ci90, col=c("black", "red", "red"), xlab="x",

ylab="S(x)", lty=1)

0 5 10 15 20

0.0

0.2

0.4

0.6

0.8

1.0

x

S(x

)



Wir wollen nun aber nicht nur die allgemeine Survival-Zeit schätzen sondern auchüberprüfen, ob eventuell festgestellte Unterschiede auch statistisch signifikant sind.Wir verwenden hierfür den Datensatz leukemia, der im Paket survival zurVerfügung gestellt wird.

data(leukemia)

Mittels

?leukemia

können Sie mehr über Inhalt und Quelle der Daten erfahren.


Nichtparametrische Schätzung: “Kaplan-Meier”Der erste Schritt der Survival Analyse in R besteht immer darin, ein Survivalobjekt zu definieren:

survleuk <- Surv(leukemia$time, leukemia$status)survleuk

## [1] 9 13 13+ 18 23 28+ 31 34 45+ 48 161+ 5 5 8## [15] 8 12 16+ 23 27 30 33 43 45

Der Kaplan-Meier-Schätzer wird nun wiederum durch survfit berechnet. Allerdings steht nun rechts der Tilde dieVariable, nach der wie eine Unterscheidung treffen wollen; in diesem Falle nennt sich die Variable x, welche anzeigt,ob Chemotherapie durchgeführt wurde oder nicht.

km.leukemia <- survfit(survleuk ~ leukemia$x)km.leukemia

## Call: survfit(formula = survleuk ~ leukemia$x)#### records n.max n.start events median 0.95LCL## leukemia$x=Maintained 11 11 11 7 31 18## leukemia$x=Nonmaintained 12 12 12 11 23 8## 0.95UCL## leukemia$x=Maintained NA## leukemia$x=Nonmaintained NA

summary(km.leukemia)

## Call: survfit(formula = survleuk ~ leukemia$x)#### leukemia$x=Maintained## time n.risk n.event survival std.err lower 95% CI upper 95% CI## 9 11 1 0.909 0.0867 0.7541 1.000## 13 10 1 0.818 0.1163 0.6192 1.000## 18 8 1 0.716 0.1397 0.4884 1.000## 23 7 1 0.614 0.1526 0.3769 0.999## 31 5 1 0.491 0.1642 0.2549 0.946## 34 4 1 0.368 0.1627 0.1549 0.875## 48 2 1 0.184 0.1535 0.0359 0.944#### leukemia$x=Nonmaintained## time n.risk n.event survival std.err lower 95% CI upper 95% CI## 5 12 2 0.8333 0.1076 0.6470 1.000## 8 10 2 0.6667 0.1361 0.4468 0.995## 12 8 1 0.5833 0.1423 0.3616 0.941## 23 6 1 0.4861 0.1481 0.2675 0.883## 27 5 1 0.3889 0.1470 0.1854 0.816## 30 4 1 0.2917 0.1387 0.1148 0.741## 33 3 1 0.1944 0.1219 0.0569 0.664## 43 2 1 0.0972 0.0919 0.0153 0.620## 45 1 1 0.0000 NaN NA NA

Die beiden Überlebenskurven können wir natürlich graphisch zeigen:



plot(km.leukemia, col = c("red", "blue"))legend("topright", col = c("red", "blue"), lty = 1, legend = c("Maintained",

"Nonmaintained"))

0 50 100 150

0.0

0.2

0.4

0.6

0.8

1.0 MaintainedNonmaintained


gruppe1 <- 1:km.leukemia$strata[1]gruppe2 <- (km.leukemia$strata[1] + 1):(km.leukemia$strata[1] + km.leukemia$strata[2])plot(x = km.leukemia$time[gruppe1], y = km.leukemia$surv[gruppe1], type = "s",

xlab = "Alter x", ylab = "S(x)", ylim = c(0, 1), col = "red", lwd = 3)lines(x = km.leukemia$time[gruppe1], y = km.leukemia$upper[gruppe1], type = "s",

col = "red", lty = 2)lines(x = km.leukemia$time[gruppe1], y = km.leukemia$lower[gruppe1], type = "s",

col = "red", lty = 2)lines(x = km.leukemia$time[gruppe2], y = km.leukemia$upper[gruppe2], type = "s",

col = "blue", lty = 2)lines(x = km.leukemia$time[gruppe2], y = km.leukemia$lower[gruppe2], type = "s",

col = "blue", lty = 2)lines(x = km.leukemia$time[gruppe2], y = km.leukemia$surv[gruppe2], type = "s",

col = "blue", lty = 1, lwd = 5)legend("topright", legend = c("Chemo JA", "Chemo NEIN"), col = c("red", "blue"),

lwd = 3, bg = "lightblue")

50 100 150

0.0

0.2

0.4

0.6

0.8

1.0

Alter x

S(x

)

Chemo JAChemo NEIN


Einlesen “unseres” NHIS-Datensatzes

nhis <- read.table("dataset2014.txt", header = TRUE, sep = ",")

head(str(nhis)) # Struktur des Datensatzes

## 'data.frame': 103477 obs. of 40 variables:## $ pid : num 2e+11 2e+11 2e+11 2e+11 2e+11 ...## $ surveyyear : int 1997 1997 1997 1997 1997 1997 1997 1997 1997 1997 ...## $ hhx : int 308 308 308 308 309 309 309 309 310 310 ...## $ fmy : int 1 1 1 1 1 1 2 2 1 1 ...## $ px : int 1 2 3 4 1 2 3 6 1 2 ...## $ sex : int 2 1 1 2 1 2 2 2 1 2 ...## $ age : int 33 36 12 4 42 52 34 2 41 36 ...## $ mob : int 11 11 6 7 8 6 99 7 11 12 ...## $ yob : int 1963 1960 1984 1992 1954 1944 9999 1994 1955 1960 ...## $ hisp : int 1 1 1 1 1 1 1 1 1 1 ...## $ ethnicity : int 1 1 1 1 1 1 1 1 1 1 ...## $ marstat : int 1 1 0 0 1 1 4 0 1 1 ...## $ region : int 4 4 4 4 4 4 4 4 4 4 ...## $ bigcity : int 1 1 1 1 1 1 1 1 1 1 ...## $ adl : int 2 2 2 2 2 2 2 NA 2 2 ...## $ iadl : int 2 2 2 NA 2 2 2 NA 2 2 ...## $ worknow : int 2 2 NA NA 2 2 2 NA 2 2 ...## $ walkprob : int 2 2 2 2 2 2 2 2 2 2 ...## $ memory : int 2 2 2 2 2 2 2 2 2 2 ...## $ anylimit : int 0 0 0 0 0 0 0 0 0 0 ...## $ selfrthealth : int 1 1 1 1 3 2 4 3 1 1 ...## $ meddelay : int 2 2 2 2 2 2 2 2 2 2 ...## $ mednotafford : int 2 2 2 2 2 2 2 2 2 2 ...## $ hospnights : int NA NA NA NA NA NA NA NA NA 1 ...## $ doctor10plus : int 2 2 2 2 2 2 2 2 2 1 ...## $ medicare : int 3 3 3 3 3 3 3 3 3 3 ...## $ medicare2 : int NA NA NA NA NA NA NA NA NA NA ...## $ medicaid : int 3 3 1 1 3 3 1 1 3 3 ...## $ healthcarecost: int NA NA NA NA 3 3 NA NA 3 3 ...## $ usborn : int 2 2 2 1 2 2 2 1 2 2 ...## $ edu : int 9 6 5 96 8 9 9 96 6 6 ...## $ jobstatus : int 1 1 NA NA 1 4 1 NA 1 4 ...## $ earnings : int 2 1 NA NA 6 1 2 NA 6 NA ...## $ famincome : int 3 3 3 3 6 6 2 2 6 6 ...## $ incpov : int 3 3 3 3 9 9 2 2 6 6 ...## $ homestatus : int 4 4 4 4 2 2 3 3 2 2 ...## $ eligible : int 1 1 2 2 1 1 1 2 1 1 ...## $ status : int 0 0 NA NA 1 0 0 NA 0 0 ...## $ quarter : int NA NA NA NA 3 NA NA NA NA NA ...## $ year : int NA NA NA NA 1999 NA NA NA NA NA ...## NULL



# Exklusion derjenigen Faelle, die NICHT eligible fÃ¼r Follow-Up waren:table(nhis$status)

#### 0 1## 61256 6735

nhis2 <- subset(nhis, !is.na(nhis$status))

table(nhis2$age)

#### 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32## 1212 1191 1172 1137 1120 1099 1135 1201 1313 1416 1379 1354 1376 1414 1536## 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47## 1501 1565 1510 1649 1538 1575 1600 1581 1543 1548 1480 1331 1355 1315 1277## 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62## 1273 1286 1290 1023 977 966 1060 840 839 829 801 797 784 654 663## 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77## 656 668 709 677 699 660 644 661 616 605 593 520 543 479 480## 78 79 80 81 82 83 84 85## 399 419 354 304 297 252 253 998



# Ein Problem sind also immer noch die Personen mit dem Alter 85, was ja 85+# entspricht:

nhis3 <- subset(nhis2, age < 85)

# Ein weiteres Problem ist das Geburtsjahr, wie ich nach langer Zeit# herausfinden konnte:table(nhis3$yob)

#### 1913 1914 1915 1916 1917 1918 1919 1920 1921 1922 1923 1924 1925 1926 1927## 375 265 291 330 391 414 414 502 493 525 517 643 612 624 628## 1928 1929 1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940 1941 1942## 653 692 704 658 697 665 633 653 726 777 821 803 826 820 933## 1943 1944 1945 1946 1947 1948 1949 1950 1951 1952 1953 1954 1955 1956 1957## 1030 959 989 1151 1307 1226 1267 1323 1323 1338 1364 1530 1503 1574 1571## 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 1971 1972## 1612 1496 1614 1518 1583 1541 1498 1489 1336 1380 1372 1373 1367 1265 1142## 1973 1974 1975 1976 1977 1978 1979 9997 9998 9999## 1115 1118 1081 1183 1160 1181 651 49 19 310

nhis4 <- subset(nhis3, yob < 9000)


Hausaufgabe

Erstellen Sie ein Survivalobjekt für die Prozesszeit Alter!

s.objekt.alter <- Surv(alter.rein, alter.raus, derstatus)

Wie definiere ich die drei Argumente der Funktion Surv richtig?


Klein, J. P. and M. L. Moeschberger (2003). Survival Analysis : Techniquesfor Censored and Truncated Data. Statistics for Biology and Health. NewYork, NY: Springer.


Lizenz

This open-access work is published under the terms of the CreativeCommons Attribution NonCommercial License 2.0 Germany, whichpermits use, reproduction & distribution in any medium for non-commercialpurposes, 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)


http://creativecommons.org/licenses/by-nc/2.0/de/

http://creativecommons.org/licenses/by-nc/2.0/de/deed.en

Kontakt

Universität RostockInstitut für Soziologie und DemographieLehrstuhl für DemographieUlmenstr. 6918057 RostockGermany

Tel.: +49-381-498 4044Fax.: +49-381-498 4395Email: [email protected]

Sprechstunde im Sommersemester 2014: Mittwochs, 09:00–10:00(und nach Vereinbarung)


Survival Analysis (Modul: Lebensdaueranalyse) · Nichtparametrische Schätzung: “Kaplan-Meier”...

Documents

Transcript of Survival Analysis (Modul: Lebensdaueranalyse) · Nichtparametrische Schätzung: “Kaplan-Meier”...