compCP=function(time,cens,group=rep(1,length(time)))
{
########################################################
# this function compares the CP-s of two groups        #
# it is based on the test presented in the paper       #
# from Statistics in Medicine, 1993, by Pepe and Mori  #
########################################################
ttau=unlist(by(time,group,max))
tau=min(ttau)
ng=table(group)
lg=labels(ng)$group
dd=table(time,cens,group)
if (dim(dd)[3]!=2) stop("Pepe-Mori test is for two groups")
if (dim(dd)[2]>3) stop("All competing risks should be grouped under  code 2")
if (dim(dd)[2]<=2) stop("There are either no censored obs or \nno competing risks or no events of interest")
tt=sort(unique(time))
tt1=c(tt[2:length(tt)],NA)
nt=table(tt1<=tau)[2]
deltat=tt1-tt
dd=dd[1:nt,,]
deltat=deltat[1:nt]
tt=tt[1:nt]

for (i in 1:2)
{

dd1=dd[,,i] ## table of censor, ev, cr
dd2=apply(dd1,1,sum) ## sum for each time point
nrisk=ng[i]-cumsum(c(0,dd2[1:(nt-1)]))
dev=dd1[,2]
dcr=dd1[,3]
dcens=dd1[,1]+dcr
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)
fcri=dcr/nrisk*sminus
fcr=cumsum(fcri)
cp=f/(1-fcr)
Ci=(nrisk-dcens)/nrisk
C=cumprod(Ci)
Cminus=c(1,C[1:(length(C)-1)])

if (i==1)
{nrisk1=nrisk
s1=s
f1=f
fcr1=fcr
C1=C
CP1=cp
C1minus=Cminus
dall1=dall
dev1=dev}
if (i==2) 
{nrisk2=nrisk
s2=s
f2=f
fcr2=fcr
C2=C
CP2=cp
C2minus=Cminus
dall2=dall
dev2=dev}

}

wi=C1minus*C2minus*sum(ng)/(ng[1]*C1minus+ng[2]*C2minus)
si=wi*(CP1-CP2)*deltat
s=sqrt(ng[1]*ng[2]/sum(ng))*sum(si)

# for group 1
temp=wi*s1*deltat/(1-fcr1)**2
temp[is.na(temp)]=0
tparti=rev(cumsum(rev(temp)))

sigma1i=tparti**2*(dev1*(1-fcr1)**2+(dall1-dev1)*f1**2)/(nrisk1*(nrisk1-1))
sigma1=sum(sigma1i,na.rm=T)
# for group 2
temp=wi*s2*deltat/(1-fcr2)**2
temp[is.na(temp)]=0
tparti=rev(cumsum(rev(temp)))

sigma2i=tparti**2*(dev2*(1-fcr2)**2+(dall2-dev2)*f2**2)/(nrisk2*(nrisk2-1))
sigma2=sum(sigma2i,na.rm=T)
sigma=ng[1]*ng[2]*(sigma1+sigma2)/sum(ng)
z=s^2/sigma
pvalue=1-pchisq(z,1)
r=data.frame(chisquare=z,pvalue=pvalue)
row.names(r)=""
return(r)
}