#####--------------------------------------------------------------------------- #####--------------------------------------------------------------------------- ## Daniel Wollschläger ## ergänzende Befehle zum Buch "Grundlagen der Datenanalyse mit R" (Springer) ## http://www.uni-kiel.de/psychologie/dwoll/r/ ## Stand: 19.12.12 ## ## Vgl. Abschnitt 4.1, um die folgenden Befehle auszuführen. ## ## Um alle Teile des Skripts ausführen zu können, muss das Paket ## survival ## samt Abhängigkeiten installiert sein, vgl. Abschnitt 1.2.7: ## ## install.packages("survival", repos="http://ftp5.gwdg.de/pub/misc/cran/") #####--------------------------------------------------------------------------- #####--------------------------------------------------------------------------- ###*--------------------------------------------------------------------------- ##### 9 Analyse von Überlebenszeiten #####--------------------------------------------------------------------------- #####--------------------------------------------------------------------------- ##### 9.1 Verteilung von Überlebenszeiten #####--------------------------------------------------------------------------- ##### 9.2 Zensierte und gestutzte Überlebenszeiten ##### 9.2.1 Zeitlich konstante Prädiktoren N <- 180 P <- 3 sex <- factor(sample(c("f", "m"), N, replace=TRUE)) X <- rnorm(N, 0, 1) IV <- factor(rep(LETTERS[1:P], each=N/P)) IVeff <- c(0, -1, 1.5) Xbeta <- 0.7*X + IVeff[unclass(IV)] + rnorm(N, 0, 2) weibA <- 1.5 weibB <- 100 U <- runif(N, 0, 1) eventT <- ceiling((-log(U)*(weibB^weibA)*exp(-Xbeta))^(1/weibA)) obsLen <- 120 plot(ecdf(eventT), xlim=c(0, 200), main="Kumulative Überlebenszeit-Verteilung", xlab="t", ylab="F(t)", cex.lab=1.4) abline(v=obsLen, col="blue", lwd=2) text(obsLen-5, 0.2, adj=1, labels="Ende Beobachtungszeit", cex=1.4) censT <- rep(obsLen, N) obsT <- pmin(eventT, censT) status <- eventT <= censT dfSurv <- data.frame(obsT, status, sex, X, IV) ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) Surv(obsT, status) ##### 9.2.2 Daten in Zählprozess-Darstellung ##### 9.2.2.1 Wiederkehrende Ereignisse ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) dfSurvCP <- survSplit(Surv(obsT, status) ~ ., cut=seq(30, 90, by=30), start="start", id="ID", zero=0, data=dfSurv) idxOrd <- order(dfSurvCP$ID, dfSurvCP$start) head(dfSurvCP[idxOrd, ], n=7) ##### 9.2.2.2 Zeitabhängige Prädiktoren ##### 9.3 Kaplan-Meier-Analyse ##### 9.3.1 Survival-Funktion und kumulatives hazard schätzen ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) KM0 <- survfit(Surv(obsT, status) ~ 1, type="kaplan-meier", conf.type="log", data=dfSurv) (KM <- survfit(Surv(obsT, status) ~ IV, type="kaplan-meier", conf.type="log", data=dfSurv)) quantile(KM0, probs=c(0.25, 0.5, 0.75), conf.int=TRUE) summary(KM0, times=c(20, 50, 100) plot(KM0, main=expression(paste("KM-Schätzer ", hat(S)(t), " mit CI")), xlab="t", ylab="Survival", lwd=2) plot(KM0, main=expression(paste("KM-Schätzer ", hat(S)(t), " mit CI")), xlab="t", ylab="Survival", fun=function(x) { 1-x }, lwd=2) plot(KM, main=expression(paste("KM-Schätzer ", hat(S)[g](t), " für Gruppen")), xlab="t", ylab="Survival", lwd=2, col=1:3) legend(x="topright", col=1:3, lwd=2, legend=LETTERS[1:3]) plot(KM0, main=expression(paste("KM-Schätzer ", hat(Lambda)(t))), xlab="t", ylab="kumulatives hazard", fun="cumhaz", lwd=2) ##### 9.3.2 Test auf gleiche Survival-Funktionen ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") survdiff(Surv(obsT, status) ~ IV, data=dfSurv) survdiff(Surv(obsT, status) ~ IV + strata(sex), data=dfSurv) ##### 9.4 Cox Proportional Hazards Modell ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) (fitCPH <- coxph(Surv(obsT, status) ~ X + IV, data=dfSurv)) coxph(Surv(start, obsT, status) ~ X + IV, data=dfSurvCP) summary(fitCPH) ##### 9.4.1 Anpassungsgüte und Modelltests ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) extractAIC(fitCPH) LLf <- fitCPH$loglik[2] LL0 <- fitCPH$loglik[1] as.vector( 1 - (LLf / LL0)) # R^2 McFadden as.vector( 1 - exp((2/N) * (LL0 - LLf))) # R^2 Cox & Snell as.vector((1 - exp((2/N) * (LL0 - LLf))) / (1 - exp(LL0)^(2/N))) fitCPH1 <- coxph(Surv(obsT, status) ~ X, data=dfSurv) anova(fitCPH1, fitCPH) ##### 9.4.2 Schätzung der Survival-Funktion und des baseline hazards ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) (CPH <- survfit(fitCPH, conf.type="log")) quantile(CPH, probs=c(0.25, 0.5, 0.75), conf.int=FALSE) dfNew <- data.frame(sex=factor(c("f", "f"), levels=levels(dfSurv$sex)), X=c(-2, -2), IV=factor(c("A", "C"), levels=levels(dfSurv$IV))) CPHnew <- survfit(fitCPH, newdata=dfNew) dev.new(width=14, height=7) par(mfrow=c(1, 2), mar=c(5, 4.5, 4, 2)+0.1) plot(CPH, main=expression(paste("Cox PH-Schätzung ", hat(S)(t), " mit CI")), xlab="t", ylab="Survival", lwd=2) lines(CPHnew$time, CPHnew$surv[ , 1], lwd=2, col="blue") lines(CPHnew$time, CPHnew$surv[ , 2], lwd=2, col="red") legend(x="topright", lwd=2, col=c("black", "blue", "red"), legend=c("pseudo-Beobachtung", "sex=f, X=-2, IV=A", "sex=f, X=-2, IV=C")) expCoef <- exp(coef(fitCPH)) Lambda0A <- basehaz(fitCPH, centered=FALSE) Lambda0B <- expCoef[2]*Lambda0A$hazard Lambda0C <- expCoef[3]*Lambda0A$hazard plot(hazard ~ time, main=expression(paste("Cox PH-Schätzung ", hat(Lambda)[g](t), " pro Gruppe")), type="s", ylim=c(0, 5), xlab="t", ylab="kumulatives hazard", lwd=2, data=Lambda0A) lines(Lambda0A$time, Lambda0B, lwd=2, col="red") lines(Lambda0A$time, Lambda0C, lwd=2, col="green") legend(x="bottomright", lwd=2, col=1:3, legend=LETTERS[1:3]) ##### 9.4.3 Modelldiagnostik ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) (czph <- cox.zph(fitCPH)) plot(czph) dfbetas <- residuals(fitCPH, type="dfbetas") plot(dfbetas[ , 1], type="h", main="DfBETAS für X", ylab="DfBETAS", lwd=2) plot(dfbetas[ , 2], type="h", main="DfBETAS für IV-B", ylab="DfBETAS", lwd=2) plot(dfbetas[ , 3], type="h", main="DfBETAS für IV-C", ylab="DfBETAS", lwd=2) resMart <- residuals(fitCPH, type="martingale") plot(dfSurv$X, resMart, main="Martingal-Residuen bzgl. X", xlab="X", ylab="Residuen", pch=20) lines(loess.smooth(dfSurv$X, resMart), lwd=2, col="blue") legend(x="bottomleft", col="blue", lwd=2, legend="LOESS fit") ##### 9.4.4 Vorhersage und Anwendung auf neue Daten ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) predict(fitCPH, type="risk") ##### 9.5 Parametrische Proportional Hazards Modelle ##### 9.5.1 Anpassungsgüte und Modelltests ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") library(survival) fitWeib <- survreg(Surv(obsT, status) ~ X + IV, dist="weibull", data=dfSurv) summary(fitWeib) fitExp <- survreg(Surv(obsT, status) ~ X + IV, dist="exponential", data=dfSurv) anova(fitExp, fitWeib) fitR <- survreg(Surv(obsT, status) ~ X, dist="weibull", data=dfSurv) anova(fitR, fitWeib) ##### 9.5.2 Survival-Funktion schätzen ## für Daten des Buches diese Datei laden: # load("data_survival.Rdata") dfNew <- data.frame(sex=factor(c("m", "m"), levels=levels(dfSurv$sex)), X=c(0, 0), IV=factor(c("A", "C"), levels=levels(dfSurv$IV))) percs <- (1:99)/100 FWeib <- predict(fitWeib, newdata=dfNew, type="quantile", p=percs, se=TRUE) matplot(cbind(FWeib$fit[1, ], FWeib$fit[1, ] - 2*FWeib$se.fit[1, ], FWeib$fit[1, ] + 2*FWeib$se.fit[1, ]), 1-percs, type="l", main=expression(paste("Weibull-Fit ", hat(S)(t), " mit SE")), xlab="t", ylab="Survival", lty=c(1, 2, 2), lwd=2, col="blue") matlines(cbind(FWeib$fit[2, ], FWeib$fit[2, ] - 2*FWeib$se.fit[2, ], FWeib$fit[2, ] + 2*FWeib$se.fit[2, ]), 1-percs, col="red", lwd=2) legend(x="topright", lwd=2, lty=c(1, 2, 1, 2), col=c("blue", "blue", "red", "red"), legend=c("sex=m, X=0, IV=A", "+- 2*SE", "sex=m, X=0, IV=C", "+- 2*SE")) ## entferne die (ggf. automatisch) geladenen Zusatzpakete, sofern möglich try(detach(package:survival)) try(detach(package:splines))