﻿# 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".