# PRACTICAL EXERCISE 3
# ====================

# Read data:path="http://www.uio.no/studier/emner/matnat/math/STK4080/h14/melanoma.txt"
leukemia=read.table("http://folk.uio.no/borgan/abg-2008/data/leukemia.txt", header=T)

# We use the survival library, so this has to be loaded.


# Question a
# ----------
# We run the commands given in the exercise text (slightly modified):

fit.b=coxph(Surv(time,status)~1,data=leukemia,subset=(treat==1),method="breslow")
surv.b=survfit(fit.b)
fit.e=coxph(Surv(time,status)~1,data=leukemia,subset=(treat==1),method="efron")
surv.e=survfit(fit.e)
relapse.time= surv.b$time
nelson.aalen.b=-log(surv.b$surv)
nelson.aalen.e=-log(surv.e$surv)
cbind(relapse.time,nelson.aalen.b,nelson.aalen.e)


# We may also compute the Nelson-Aalen estimates from first principles 
# from formulas (3.12) and (3.13).

# For formula (3.12) we obtain for times 1, 2, 3, and 4:
naa.3.12=cumsum(c(1/21+1/20,1/19+1/18,1/17, 1/16+1/15))

# For formula (3.13) we obtain for times 1, 2, 3, and 4:
naa.3.13=cumsum(c(2/21,2/19,1/17,2/16))

# We then display the results:
cbind(1:4,naa.3.12,naa.3.13)

# Comparing this with the output above, we see that "efron" 
# corresponds to (3.12) and "breslow" corresponds to (3.13)


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

# We run the commands that give the estimated variances (slightly modified):
var.b=surv.b$std.err^2
var.e= surv.e$std.err^2
cbind(relapse.time,var.b,var.e)

# We then compute the variance estimates (3.15) and (3.16) from first principles.

# For formula (3.15) we obtain for times 1, 2, 3, and 4:
var.3.15=cumsum(c(1/21^2+1/20^2,1/19^2+1/18^2,1/17^2, 1/16^2+1/15^2))

# For formula (3.16) we obtain for times 1, 2, 3, and 4:
var.3.16=cumsum(c((21-2)*2/21^3,(19-2)*2/19^3,(17-1)*1/17^3,(16-2)*2/16^3))

# We also compute the variance estimate that has increment 
# d_j/Y(T_j)^2 at time T_j for the first four failure times:
var.new=cumsum(c(2/21^2,2/19^2,1/17^2,2/16^2))

# We then display the results:
cbind(1:4,var.3.15,var.3.16,var.new)

# Comparing this with the output above, we see that "efron" corresponds
# to the variance estimator (3.15), while "breslow" doesn't give (3.16)
# but the alternative estimator "var.new".