#########################################
# CIF variance based on the delta method
# Author: Melanie Pintilie
#########################################

cifDM=function(time,cens,group=rep(1,length(time)))
{
# this function calculates the CIF estimates and #
# their variance based on the delta method #
lg=labels(table(group))$group
ng=length(lg)
dd=table(time,cens,group)
tt=sort(unique(time))
for (i in 1:ng)
{
dd1=dd[,,i] ## table of censor, ev, cr
dd2=apply(dd1,1,sum) ## sum for each time point
nrisk=sum(dd2)-cumsum(dd2)
nrisk=c(sum(dd2),nrisk[1:(length(nrisk)-1)])
dev=dd1[,2]
dcr=dd1[,3]
dall=dev+dcr
si=(nrisk-dall)/nrisk
s=cumprod(si)
sminus=c(1,s[1:(length(s)-1)])
fi=dev/nrisk*sminus
f=cumsum(fi)
t1i=dall/(nrisk*(nrisk-dall))
t2i=f*t1i
t3i=f^2*t1i
t4i=sminus^2*(nrisk-dev)*dev/nrisk^3
t5i=sminus*dev/nrisk^2
t6i=f*t5i
t1=cumsum(t1i)
t2=cumsum(t2i)
t3=cumsum(t3i)
t4=cumsum(t4i)
t5=cumsum(t5i)
t6=cumsum(t6i)
v=f^2*t1-2*f*t2+t3+t4-2*f*t5+2*t6
res=data.frame(time=tt,cif=f,var=v,group=rep(lg[i],
length(v)))
res1=res[dev>0,]
if (i==1) cifdelta=res1
else cifdelta=list(cifdelta,res1)
}
return(cifdelta)
}