# PRACTICAL EXERCISE 4
# ====================

# Read leukemia data:
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 compute the estimates of the survival function for the placebo group with the 
# options "kaplan-meier", "fleming-harrington", and "fh2" (the first is default)

surv.km=survfit(Surv(time,status)~1,data=leukemia,subset=(treat==1),type="kaplan-meier")
surv.fl=survfit(Surv(time,status)~1,data=leukemia,subset=(treat==1),type="fleming-harrington")
surv.fh2=survfit(Surv(time,status)~1,data=leukemia,subset=(treat==1),type="fh2")

# We then list the estimates:
cbind(surv.km$time, surv.km$surv, surv.fl$surv, surv.fh2 $surv)

# We see that estimates are close at early times when there is a fair number of individuals at risk. 
# The difference is larger at later times when there are only a few individuals at risk.
# The Kaplan-Meier estimator is smallest, and the Flemming-Harrington is largest.


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

# We then compute the two standard error estimates:
# * error="greenwood" in combination with type="kaplan-meier" (default)
# * error="tsiatis" in combination with type="fh2" 
surv.gw=survfit(Surv(time,status)~1,data=leukemia,subset=(treat==1),
              type="kaplan-meier",error="greenwood",conf.type="none")
surv.ts=survfit(Surv(time,status)~1,data=leukemia,subset=(treat==1),
                type="fh2",error="tsiatis", conf.type="none")

# We may use the summary command to see the estimated standard errors

summary(surv.gw)
summary(surv.ts)

# The see that the standard errors are quite close