# 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))