This file gives some R code useful for some homework problems from "Multiple Linear and 1D Regression" by DJ Olive. #HW16.43 Miller p. 49 time until relapse under maintenance chemotherapy #AML leukemia data library(survival) time <- c(9,13,13,18,23,28,31,34,45,48,161) status <- c(1,1,0,1,1,0,1,1,0,1,0) z <- survfit(Surv(time,status),conf.type="plain") summary(z) #plot KM estimator hat(S)_K(t) and the CIs for S(t) plot(z) #HW16.44a type library(survival) if necessary zy <- rexp(20) zz <- rexp(20,rate=.2) zt <- pmin(zy,zz) zdelta <- as.numeric(zz>zy) out <- survfit(Surv(zt, zdelta), conf.type = "plain") st <- exp( - out$time) plot(out) points(out$time,st) #HW16.44b type library(survival) if necessary #repeat with n = 200 zy <- rexp(200) zz <- rexp(200,rate=.2) zt <- pmin(zy,zz) zdelta <- as.numeric(zz>zy) out <- survfit(Surv(zt, zdelta), conf.type = "plain") st <- exp( - out$time) plot(out) points(out$time,st) #HW16.46a library(survival) z <- coxph(Surv(futime,fustat)~age + ecog.ps, ovarian) z plot(survfit(z)) #16.46b plot(survfit(z,newdata=c(66,4))) #HW16.48 library(survival) zfull <- coxph(Surv(time,status)~strata(sex) + age +ph.ecog + ph.karno+pat.karno+wt.loss, data = lung,na.action=na.omit) zred1 <- coxph(Surv(time,status)~strata(sex) + ph.ecog + ph.karno+pat.karno+wt.loss, data = lung,na.action=na.omit) zred2 <- coxph(formula = Surv(time, status) ~ strata(sex) + ph.ecog + +pat.karno + wt.loss, data = lung, na.action = na.omit) #HW 16.49 need to download bphgfit function, Smith (2002, p. 174) Lawless (1982) #leukemia data lweeks <- c(1, 3, 3, 6, 7, 7, 10, 12, 14, 15, 18, 19, 22, 26, 28, 29, 34, 40, 48, 49, 1, 1, 2, 2, 3, 4, 5, 8, 8, 9, 11, 12, 14, 16, 18, 21, 27, 31, 38, 44 ) lstatus <- 1 + 0*1:40 lstatus[c(15,19,20,37,39)] <- 0 lg <- 0*1:20 lg2 <- lg + 1 lgroup <- c(lg,lg2) leuk <- cbind(lweeks,lstatus,lgroup) rm(lweeks,lstatus,lg,lg2,lgroup) zout <- coxph(Surv(leuk[,1], leuk[,2]) ~ leuk[,3]) zout bphgfit(leuk[,1],leuk[,2],leuk[,3]) #16.50 a) status = 1 + usual status variable flung<-lung[lung$sex==1,] nflung<-na.omit(flung[,c("time","status","ph.ecog","ph.karno", "pat.karno","wt.loss")]) out <- coxph(Surv(time, status) ~ ph.ecog + ph.karno + pat.karno + wt.loss,data = nflung) x <- as.matrix(nflung[,c(3,4,5,6)]) ESP <- x %*% out$coef Time <- as.vector(nflung[,1]) status <- as.vector(nflung[,2]) plot(ESP,Time,pch = ' ') points(ESP[status == 2], Time[status == 2], pch = 3) points(ESP[status == 1], Time[status == 1], pch = 1) title("ESSP for males") #16.50b) flung<-lung[lung$sex==2,] nflung<-na.omit(flung[,c("time","status","ph.ecog","ph.karno", "pat.karno","wt.loss")]) out <- coxph(Surv(time, status) ~ ph.ecog + ph.karno + pat.karno + wt.loss,data = nflung) x <- as.matrix(nflung[,c(3,4,5,6)]) ESP <- x %*% out$coef time <- as.vector(nflung[,1]) status <- as.vector(nflung[,2]) plot(ESP,time,pch = ' ') points(ESP[status == 2], time[status == 2], pch = 3) points(ESP[status == 1], time[status == 1], pch = 1) title("ESSP for Females") #16.51a) Splus p. 269-275 library(survival) nlung<-na.omit(lung[,c("time","status","sex","ph.ecog","ph.karno", "pat.karno","wt.loss")]) par(mfrow=c(2,2)) attach(nlung) fit1 <- coxph(Surv(time, status) ~ strata(sex)+ ph.karno + pat.karno + wt.loss,data = nlung) scatter.smooth(ph.ecog,resid(fit1)) fit2 <- coxph(Surv(time, status) ~ strata(sex)+ ph.ecog + pat.karno + wt.loss,data = nlung) scatter.smooth(ph.karno,resid(fit2)) fit3 <- coxph(Surv(time, status) ~ strata(sex)+ ph.ecog + ph.karno + wt.loss,data = nlung) scatter.smooth(pat.karno,resid(fit3)) fit4 <- coxph(Surv(time, status) ~ strata(sex)+ ph.ecog + ph.karno + pat.karno,data = nlung) scatter.smooth(wt.loss,resid(fit4)) #16.51c) lung.fit2 <- coxph(Surv(time, status) ~ strata(sex)+ ph.ecog + ph.karno + pat.karno + wt.loss,data = nlung) plot(cox.zph(lung.fit2)) cox.zph(lung.fit2) #HW16.53 ovarian data Collett (2003, p. 335-6) library(survival) time <- c(156, 1040, 59, 421, 329, 769, 365, 770, 1227, 268, 475, 1129, 464, 1206, 638, 563, 1106, 431, 855, 803, 115, 744, 477, 448, 353, 377) status <- c(1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0) treat <- c(1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0) age <- c(66, 38, 72, 53, 43, 59, 64, 57, 59, 74, 59, 53, 56, 44, 56, 55, 44, 50, 43, 39, 74, 50, 64, 56, 63, 58) ovar <- cbind(time,status,treat,age) rm(time,status,treat,age) z <- survreg(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4],dist="weibull") summary(z) zfit <- 11.548 -.561*ovar[,3] -.079*ovar[,4] plot(zfit,log(ovar[,1]),pch = ' ') points(zfit[ovar[,2]==0],log(ovar[ovar[,2]==0,1]),pch=1) points(zfit[ovar[,2]==1],log(ovar[ovar[,2]==1,1]),pch=3) abline(0,1) abline( lsfit( zfit[ovar[,2]==1],log(ovar[ovar[,2]==1,1]) )$coef) # HW16.55a ovarian data Collett (2003, p. 335-6) library(survival) time <- c(156, 1040, 59, 421, 329, 769, 365, 770, 1227, 268, 475, 1129, 464, 1206, 638, 563, 1106, 431, 855, 803, 115, 744, 477, 448, 353, 377) status <- c(1, 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 0, 0, 1, 0) treat <- c(1, 1, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0) age <- c(66, 38, 72, 53, 43, 59, 64, 57, 59, 74, 59, 53, 56, 44, 56, 55, 44, 50, 43, 39, 74, 50, 64, 56, 63, 58) ovar <- cbind(time,status,treat,age) rm(time,status,treat,age) z <- survreg(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4],dist="weibull") zx <- cbind(ovar[,3],ovar[,4]) zc <- coxph(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4]) ESPC <- zx%*%zc$coef ESPW <- zx%*%z$coef[-1] sighat<-z$scale plot( - ESPW/sighat, ESPC) abline(0,1) title("a) Weibull") cor( - ESPW/sighat, ESPC) #HW 16.55b z <- survreg(Surv(ovar[,1],ovar[,2])~ovar[,3]+ovar[,4],dist="exponential") ESPE <- zx%*%z$coef[-1] shat<-z$scale plot( - ESPE, ESPC) abline(0,1) title("b) Exponential") #HW16.58 lung data, Splus p. 269-275 library(survival) sfit <- coxph(Surv(time, status)~ strata(sex)+ ph.ecog + pat.karno + wt.loss,data = lung,na.action=na.omit) summary(sfit) plot(survfit(sfit), lty = 2:3) title("F Survival is Top Curve, Ave. Covariates")