# PRACTICAL EXERCISE 7 # ==================== # Read data and disregard the individuals who smoke pipe or cigar: norw.death=read.table("http://folk.uio.no/borgan/abg-2008/data/causes_death.txt", header=T) norw.death=norw.death[norw.death\$smkgr!=6,] # We will use the survival package, so this has to be loaded library(survival) # Question a # ---------- # As in example 4.1 in the ABG-book, we consider total mortality in this problem. # The following comands reproduce the results of example 4.1 fit0=coxph(Surv(agestart,agestop,dead)~factor(sex)+factor(smkgr), data=norw.death) summary(fit0) # Note that we in the coxph-command need to include a start age, since # the individuals are not all followed from the same age. # Question b # ---------- # We then do a Cox regression with sex, smoking habit, systolic blood (sbp) pressure, # and body mass index (bmi) as covariates: fit.b=coxph(Surv(agestart,agestop,dead)~factor(sex)+factor(smkgr)+sbp+bmi, data=norw.death) summary(fit.b) # We find a clear effect of sbp, but not of bmi. Further the effect of sex is slightly reduced, # while the effects of smoking are slightly increased when correcting for sbp and bmi. # Question c # ---------- # To (e.g) check for interaction between sex and smoking habits, # we fit a model with interaction between the covariates and use the # likelihood ratio test to compare with the model of question b: fit.ss=coxph(Surv(agestart,agestop,dead)~factor(sex)+factor(smkgr)+sbp+bmi+factor(sex):factor(smkgr), data=norw.death) anova(fit.b,fit.ss) # Similar commands can be used to check for interaction between: # - sex and sbp # - sex and bmi # - smoking and sbp # - smoking and bmi # - sbp and bmi # None of the interactiona are significant at the 5% level # (but the one for sbp and bmi is close with a p-value of 5.1%) # Question d # ---------- # The following comands reproduce the results of example 4.4 new.cov=data.frame(sex=c(1,1,2,2),smkgr=c(1,5,1,5)) surv0=survfit(fit0,newdata=new.cov) plot(surv0,mark.time=F,xlim=c(40,70),xlab="Age", ylab="Survival", lty=c(1,1,2,2)) # Question e # ---------- # We will make two plots of estimated survival curves, one for men and one for women. # For each gender we consider the following smoking habits and sbp # (seting bmi=25, which is close to its mean value) # 1) smkgr=1 (never smoked), sbp=120 # 2) smkgr=1 (never smoked), sbp=150 # 3) smkgr=3 (1-9 cigarettes per day), sbp=120 # 4) smkgr=3 (1-9 cigarettes per day), sbp=150 # 5) smkgr=5 (20+ cigarettes per day), sbp=120 # 6) smkgr=5 (20+ cigarettes per day), sbp=150 # par(mfrow=c(1,2)) new.cov.males=data.frame(sex=rep(1,6),smkgr=rep(c(1,3,5),each=2), sbp=rep(c(120,150),3),bmi=rep(25,6)) surv.males=survfit(fit.b,newdata=new.cov.males) plot(surv.males,mark.time=F,xlim=c(40,70),xlab="Age", ylab="Survival", lty=rep(1:2,3),col=rep(c("green","blue","red"),each=2), main="Males") legend("bottomleft",as.character(1:6),lty=rep(1:2,3), col=rep(c("green","blue","red"),each=2)) new.cov.females=data.frame(sex=rep(2,6),smkgr=rep(c(1,3,5),each=2), sbp=rep(c(120,150),3),bmi=rep(25,6)) surv.females=survfit(fit.b,newdata=new.cov.females) plot(surv.females,mark.time=F,xlim=c(40,70),xlab="Age", ylab="Survival", lty=rep(1:2,3),col=rep(c("green","blue","red"),each=2), main="Females") legend("bottomleft",as.character(1:6),lty=rep(1:2,3), col=rep(c("green","blue","red"),each=2))