# PRACTICAL EXERCISE 8
# ====================

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

# We fit a Cox model with sex and smoking habits as categorical covariates 
# and systolic blodd pressure and body mass index as numeric covariates: 
fit.b=coxph(Surv(agestart,agestop,dead)~factor(sex)+factor(smkgr)+sbp+bmi, data=norw.death)
summary(fit.b)

# The model assumes a log-linear effect of the numeric covariates and 
# proportionality of all covariates


# Question b
# ----------

# We fit a model with smoothing spline for sbp:
fit.b=coxph(Surv(agestart,agestop,dead)~factor(sex)+factor(smkgr)+pspline(sbp)+bmi, data=norw.death)
print(fit.b)
par(mfrow=c(2,2))
termplot(fit.b,se=T)

# The effect of systolic blood pressure is fairly log-linear.


# Question c
# ----------

# We fit a model with smoothing spline for bmi:

fit.c=coxph(Surv(agestart,agestop,dead)~factor(sex)+factor(smkgr)+sbp+pspline(bmi), data=norw.death)
print(fit.c)
par(mfrow=c(2,2))
termplot(fit.c,se=T)

# The effect of body mass index in not log-linear. It seems difficult to find 
# a useful trasformation of body mass index, so we will work with grouped bmi.

# We make a categorical variable for groupe bmi (cf. problem text) and fit a model with this 
# covariate using the igroup [20,25) as reference:
norw.death$gbmi=cut(norw.death$bmi,c(0,20,25,30,60),right=F)
norw.death$gbmi=relevel(norw.death$gbmi,ref="[20,25)")
fit.c2=coxph(Surv(agestart,agestop,dead)~factor(sex)+factor(smkgr)+sbp+factor(gbmi), data=norw.death)
summary(fit.c2)


# Question d
# ----------
# We perform a test and make plots to check proportionality for the model with grouped body mass index 
cox.zph(fit.c2,transform='log')
par(mfrow=c(3,3))
plot(cox.zph(fit.c2,transform='log'))

# The assumption of proportionality seems reasonable for all covariates
# (There may be some problem with proportionality for smoking group 5, but this
#  is not serious enough to "force us" to modify the model.)