Description Details Author(s) References See Also Examples
This package supplements the Survival Analysis in R: A Tutorial paper. The tutorial describes how to apply several basic survival analysis techniques in R using the survival
package. Data sets from the KMsurv
package are used in most examples; this package is a supplement to Klein and Moeschberger's textbook (see References).
All code used in the tutorial are included in the examples below.
Package: | OIsurv |
Type: | Package |
Version: | 0.2 |
Date: | 2013-10-16 |
License: | GPL (>= 2) |
URL: | http://www.openintro.org/stat/surv.php |
LazyLoad: | yes |
David M Diez
Maintainer: David M Diez <david.m.diez@gmail.com>
Fox J (2002). "Cox Proportional-Hazards Regression for Survival Data. Appendix to An R and S-PLUS Companion to Applied Regression." Comprehensive R Archive Network. http://cran.r-project.org/doc/contrib/Fox-Companion/appendix-cox-regression.pdf
Klein JP, Moeschberger ML (2003). Survival Analysis: Techniques for Censored and Truncated Data. Springer Verlag, New York.
ReliaSoft Corporation website (2006). "Data Classification." http://www.weibull.com/LifeDataWeb/data_classification.htm
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 | #=====> 2. Three packages: survival, OIsurv, and KMsurv <=====#
# install.packages("OIsurv")
# library(OIsurv)
data(aids)
aids
attach(aids)
infect
detach(aids)
#=====> 3. Survival objects <=====#
data(tongue)
attach(tongue)
mySurvObject <- Surv(time, delta)
mySurvObject
detach(tongue)
# Surv(time, event, type="left")
# Surv(t1, t2, event)
#=====> 4. Kaplan-Meier estimate and pointwise bounds <=====#
data(tongue)
attach(tongue)
mySurv <- Surv(time[type==1], delta[type==1])
(myFit <- survfit(mySurv ~ 1))
summary(myFit)
myFit$surv # outputs the Kaplan-Meier estimate at each t_i
myFit$time # t_i
myFit$n.risk # Y_i
myFit$n.event # d_i
myFit$std.err # standard error of the K-M estimate at t_i
myFit$lower # lower pointwise estimates (alternatively, $upper)
#pdf("kmPlot.pdf", 7, 4.5)
#par(mar=c(3.9, 3.9, 2.5, 1), mgp=c(2.6, 0.7, 0))
plot(myFit, main="Kaplan-Meier estimate with 95% confidence bounds",
xlab="time", ylab="survival function")
#dev.off()
myFit1 <- survfit(Surv(time, delta) ~ type) # 'type' specifies the grouping
detach(tongue)
#=====> 5. Kaplan-Meier confidence bands <=====#
data(tongue)
attach(tongue)
mySurv <- Surv(time[type==1], delta[type==1])
#pdf("confBand.pdf", 7, 4.5)
#par(mar=c(3.9, 3.9, 2.5, 1), mgp=c(2.6, 0.7, 0))
plot(survfit(mySurv ~ 1), xlab='time',
ylab='Estimated Survival Function',
main='Confidence intervals versus confidence bands')
myCB <- confBands(mySurv)
lines(myCB, lty=3)
legend('topright', legend=c('K-M survival estimate',
'pointwise intervals','EP confidence bands'), lty=1:3)
#dev.off()
detach(tongue)
#=====> 6. Cumulative hazard <=====#
data(baboon)
attach(baboon)
mySurv <- Surv(time, observed)
myFit <- summary(survfit(mySurv ~ 1))
Hhat <- -log(c(myFit$surv, tail(myFit$surv, 1)))
hSortOf <- myFit$n.event / myFit$n.risk
Htilde <- c(cumsum(hSortOf), sum(hSortOf))
#pdf("cumHazard.pdf", 7, 4)
#par(mar=c(3.9, 3.9, 2.5, 1), mgp=c(2.6, 0.7, 0))
plot(c(myFit$time, 1200), Hhat, xlab='time', ylab='cumulative hazard',
main='Comparing cumulative hazards', type='s')
lines(c(myFit$time, 1200), Htilde, lty=2, type='s')
legend('topleft', legend=c('MLE', 'Nelson-Aalen'), lty=1:2)
#dev.off()
detach(baboon)
#=====> 7. Mean and median estimates with bounds <=====#
data(drug6mp)
attach(drug6mp)
mySurv <- Surv(t1, rep(1, 21)) # all placebo patients observed
survfit(mySurv ~ 1)
print(survfit(mySurv ~ 1), print.rmean=TRUE)
detach(drug6mp)
#=====> 8. Tests for two or more samples <=====#
data(btrial)
attach(btrial)
survdiff(Surv(time, death) ~ im) # output omitted
survdiff(Surv(time, death) ~ im, rho=1) # some output omitted
detach(btrial)
#=====> 9. Cox proportional hazards model, constant covariates <=====#
data(burn)
attach(burn)
mySurv <- Surv(T1, D1)
(coxphFit <- coxph(mySurv ~ Z1 + as.factor(Z11)))
coxphFit$coefficients # may use coxphFit$coeff instead
coxphFit$var # I^(-1), estimated cov matrix of the beta-hats
coxphFit$loglik # log-likelihood for alt and null MLEs, resp.
mySurvfit <- survfit(coxphFit)
betaHat <- coxphFit$coef
betaCov <- coxphFit$var
anova(coxphFit)
detach(burn)
#=====> 10. Cox PH model, time-dependent covariates <=====#
data(relapse)
relapse
attach(relapse)
N <- dim(relapse)[1]
t1 <- rep(0, N+sum(!is.na(inter))) # Initialize start times at 0
t2 <- rep(NA, length(t1)) # The end times for each record
e <- rep(NA, length(t1)) # Was the event censored?
g <- rep(NA, length(t1)) # Gender
PI <- rep(FALSE, length(t1)) # Initialize intervention at FALSE
R <- 1 # Row of new record
for(ii in 1:dim(relapse)[1]){
if(is.na(inter[ii])){ # no intervention, copy survival record
t2[R] <- event[ii]
e[R] <- delta[ii]
g[R] <- gender[ii]
R <- R+1
} else { # intervention, split records
g[R+0:1] <- gender[ii] # gender is same for each time
e[R] <- 0 # no relapse observed pre-intervention
e[R+1] <- delta[ii] # relapse occur post-intervention?
PI[R+1] <- TRUE # Intervention covariate, post-intervention
t2[R] <- inter[ii]-1 # End of pre-intervention
t1[R+1] <- inter[ii]-1 # Start of post-intervention
t2[R+1] <- event[ii] # End of post-intervention
R <- R+2 # Two records added
}
}
mySurv <- Surv(t1, t2, e)
coxphFit <- coxph(mySurv ~ g + PI)
detach(relapse)
#=====> 11. Accelerated failure-time models <=====#
data(larynx)
attach(larynx)
srFit <- survreg(Surv(time, delta) ~ as.factor(stage) + age, dist='weibull')
summary(srFit)
srFitExp <- survreg(Surv(time, delta) ~ as.factor(stage) + age, dist='exponential')
summary(srFitExp)
#===> Output is omitted from the commands below
srFitExp$coeff # covariate coefficients
srFitExp$icoef # intercept and scale coefficients
srFitExp$var # variance-covariance matrix
srFitExp$loglik # log-likelihood
srFit$scale # not using srFitExp (defaulted to 1)
detach(larynx)
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.