[R-sig-Epi] Arjas plot

[Tobias Janstad]

#goal: Do Arjas plot. In particular, reproduce figure 11.15 in Klein 
Moeschberger survival analysis (1997).
library(survival);
library(KMsurv);
data(bmt)
#143 bone marrow patients
#goal: reproduce figure 11.15 in Klein Moeschberger survival analysis 
(1997).

#(i) divide the data material into two subgroups
MTX<-subset(bmt,z8==1)#MTX used as an GVHD prophylactit.
noMTX<-subset(bmt,z8==0)#

#(ii) sort with respect to survival time within each group
#allready done

#(iii)
#N_g(t_i), following notation from page 370 in K.M. 1997
#d3 is the event variable.
MTXcumdelta<-cumsum(MTX$d3)
noMTXcumdelta<-cumsum(noMTX$d3)
coxMTX=coxph(Surv(MTX$t1,MTX$d3)~z1+z2+z3+z4+z5+z6+z7,data=MTX,method="breslow")
print("Cox-model fitted for MTX-grupp.")
coxnoMTX=coxph(Surv(noMTX$t1,noMTX$d3)~z1+z2+z3+z4+z5+z6+z7,data=noMTX,method="breslow")
print("Cox-model fitted for noMTX-grupp.")

MTXbase<-basehaz(coxMTX,centered=TRUE)#estimate the cumulative baseline 
hazard (I guess)
noMTXbase<-basehaz(coxnoMTX,centered=TRUE)

MTXcumbase<-cumsum(MTXbase$hazard)
noMTXcumbase<-cumsum(noMTXbase$hazard)

plot(stepfun(MTXcumdelta[1:length(MTXcumbase)],c(0,MTXcumbase)))
lines(stepfun(noMTXcumdelta[1:length(noMTXcumbase)],c(0,noMTXcumbase)))
lines(c(0,20),c(0,20))
#Why do length of cumulative hazard and events differ?
#The plot does not look like figure 11.15. Can someone help me correct it?

Best regards, Tobias Janstad