######### Motor Data: Accelerated Life Testing ###### # The dataset "motors" contains the results of an accelerated # life test experiment with 10 replicates at each of four # temperatures reported by Kalbfleisch and Prentice (1980). # The variables are as follows: # temp: the temperature (degrees C) of the test # time: the time in hours to failure or censoring at 8064 hours # (=336 days). # cens: an indicator variable for failure, 1=failed library(survival) library(MASS) data(motors) attach(motors) # Here is the dataset > motors temp time cens 1 150 8064 0 2 150 8064 0 3 150 8064 0 4 150 8064 0 5 150 8064 0 6 150 8064 0 7 150 8064 0 8 150 8064 0 9 150 8064 0 10 150 8064 0 11 170 1764 1 12 170 2772 1 13 170 3444 1 14 170 3542 1 15 170 3780 1 16 170 4860 1 17 170 5196 1 18 170 5448 0 19 170 5448 0 20 170 5448 0 21 190 408 1 22 190 408 1 23 190 1344 1 24 190 1344 1 25 190 1440 1 26 190 1680 0 27 190 1680 0 28 190 1680 0 29 190 1680 0 30 190 1680 0 31 220 408 1 32 220 408 1 33 220 504 1 34 220 504 1 35 220 504 1 36 220 528 0 37 220 528 0 38 220 528 0 39 220 528 0 40 220 528 0 # Here are the survival times, 8064+ means that the survival time is # censored at time 8064 (so the actual but unobserved survival time is # greater than 8064), etc. Surv(time, cens) [1] 8064+ 8064+ 8064+ 8064+ 8064+ 8064+ 8064+ 8064+ 8064+ 8064+ 1764 2772 [13] 3444 3542 3780 4860 5196 5448+ 5448+ 5448+ 408 408 1344 1344 [25] 1440 1680+ 1680+ 1680+ 1680+ 1680+ 408 408 504 504 504 528+ [37] 528+ 528+ 528+ 528+ # Computes an estimate of the survival curve for censored data # using the Kaplan-Meier method fit1 <- survfit(Surv(time,cens) ~ factor(temp)) # The following summeary shows the observed survival times, # number at risk at each time, number of events at each time, # estimated survival probability and its standard error, # the lower and upper confidence limit of the estimated survival # probabilities, for each level of "temp". # Note: for level "temp=150", all survival times are censored, # so we omit the corresponding estimates. > summary(fit1) Call: survfit(formula = Surv(time, cens) ~ factor(temp)) factor(temp)=170 time n.risk n.event survival std.err lower 95% CI upper 95% CI 1764 10 1 0.9 0.0949 0.732 1.000 2772 9 1 0.8 0.1265 0.587 1.000 3444 8 1 0.7 0.1449 0.467 1.000 3542 7 1 0.6 0.1549 0.362 0.995 3780 6 1 0.5 0.1581 0.269 0.929 4860 5 1 0.4 0.1549 0.187 0.855 5196 4 1 0.3 0.1449 0.116 0.773 factor(temp)=190 time n.risk n.event survival std.err lower 95% CI upper 95% CI 408 10 2 0.8 0.126 0.587 1.000 1344 8 2 0.6 0.155 0.362 0.995 1440 6 1 0.5 0.158 0.269 0.929 factor(temp)=220 time n.risk n.event survival std.err lower 95% CI upper 95% CI 408 10 2 0.8 0.126 0.587 1.00 504 8 3 0.5 0.158 0.269 0.93 # Kaplan-Meier estimators par(mfrow=c(2,1)) plot(fit1, xlab="time", ylab="survival", main="Kaplan-Meier estimators (Motor data)", lty=1:4) legend(6000,0.8, c("temp=150","temp=170","temp=190","temp=220"),lty=1:4) plot(fit1, xlab="time", ylab="survival", main="Kaplan-Meier estimators (Motor data, with 95% confidence bands)", conf=T, lty=1:4) legend(6000,0.8, c("temp=150","temp=170","temp=190","temp=220"),lty=1:4) # There may be significant difference between group "temp=170" # and group "temp=190". #### VA data: Veterans Administration lung cancer trial #### # stime: survival time of cancer patients # status: censoring indicator # treat: treatment (standard or test) # age: in years # Karn: Karnofsky score of performance on scale 0-100, # with high values for relatively well patients # diag.time: time since diagnosis in months at entry to trial # cell: one of four cell types > data(VA) > VA stime status treat age Karn diag.time cell prior 1 72 1 1 69 60 7 1 0 2 411 1 1 64 70 5 1 10 3 228 1 1 38 60 3 1 0 4 126 1 1 63 60 9 1 10 5 118 1 1 65 70 11 1 10 6 10 1 1 49 20 5 1 0 7 82 1 1 69 40 10 1 10 8 110 1 1 68 80 29 1 0 9 314 1 1 43 50 18 1 0 10 100 0 1 70 70 6 1 0 11 42 1 1 81 60 4 1 0 12 8 1 1 63 40 58 1 10 13 144 1 1 63 30 4 1 0 14 25 0 1 52 80 9 1 10 15 11 1 1 48 70 11 1 10 ........... attach(VA) Surv(stime, status) [1] 72 411 228 126 118 10 82 110 314 100+ 42 8 144 25+ 11 ...... fit2 <- survfit(Surv(stime,status) ~ treat) > summary(fit2) Call: survfit(formula = Surv(stime, status) ~ treat) treat=1 time n.risk n.event survival std.err lower 95% CI upper 95% CI 3 69 1 0.9855 0.0144 0.95771 1.000 4 68 1 0.9710 0.0202 0.93223 1.000 7 67 1 0.9565 0.0246 0.90959 1.000 8 66 2 0.9275 0.0312 0.86834 0.991 10 64 2 0.8986 0.0363 0.83006 0.973 ....... # Kaplan-Meier estimators par(mfrow=c(2,1)) plot(fit2, xlab="time", ylab="survival", main="Kaplan-Meier estimators (VA data)",lty=1:2) legend(600,0.8, c("treatment: standard", "treatment: test"), lty=1:2) # The two survival curves seem crossing plot(fit2, xlab="time", ylab="survival", conf=T, main="Kaplan-Meier estimators (VA data, with 95% confidence band)",lty=1:2) legend(600,0.8, c("treatment: standard", "treatment: test"), lty=1:2) # There may be no significance difference between the two groups