# R-commands for Nelson-Aalen estimates
# =====================================

# We may use the coxph-command (with no covariates) to obtain a 
# Nelson-Aalen estimate (alternatively one may use the survfit-command)

# We will illustrate the commands by computing the Nelson-Aalen estimate 
# with standard errors for the melanoma data

# We first read the data:
melanoma=read.table("http://www.uio.no/studier/emner/matnat/math/STK4080/h12/melanoma.txt", header=T)

# We will use the survival library:
library(survival)

# We will make plots for females and males in the same panel:
par(mfrow=c(1,2))

# We compute and plot the Nelson-Aalen estimate for females:
fit.f=coxph(Surv(lifetime,status==1)~1,data=melanoma,subset=(sex==1))
surv.f=survfit(fit.f)
plot(surv.f,fun="cumhaz", mark.time=F,xlim=c(0,10),ylim=c(0,0.80),
    main="Females",xlab="Years since operation",ylab="Cumulative hazard")

# And then for males:
fit.m=coxph(Surv(lifetime,status==1)~1,data=melanoma,subset=(sex==2))
surv.m=survfit(fit.m)	
plot(surv.m,fun="cumhaz", mark.time=F, xlim=c(0,10),ylim=c(0,0.80),
     main="Males", xlab="Years since operation",ylab="Cumulative hazard")



# Finally we compute and plot Nelson-Aalen estimates for both genders:
par(mfrow=c(1,1))
fit.f=coxph(Surv(lifetime,status==1)~strata(sex),data=melanoma)
surv.f=survfit(fit.f)
plot(surv.f,fun="cumhaz", mark.time=F,xlim=c(0,10),ylim=c(0,0.70),
     xlab="Years since operation",ylab="Cumulative hazard",lty=1:2)