#Two-sample test using multiple imputation, which is
#implemented by Poor-man's data augmentation
# Wei Pan, Division of Biostatistics, University of Minnesota
# weip@biostat.umn.edu
# Feb 2000

mi2test2<-function(datafile, M=10, path=""){
#datafile: data file name
#M: imputation numbers
#path: directory path from current directory to that of the datafile

foutname<-paste("mi2test2.res",M, sep="")
G0<-G1<-G12<-g0<-g1<-g12<-v0<-v1<-v12<-rep(0,M)

dat<-read.table(datafile)

#calling C to calculate NPMLE:
#Group 1:
lc0<-c(0,dat$V1[dat$V3==0])
rc0<-c(0,dat$V2[dat$V3==0])
N0<-length(lc0)-1
storage.mode(lc0)<-"double"
storage.mode(rc0)<-"double"
storage.mode(N0)<-"integer"
ret0<-.C("IcmNPMLE", x=lc0, S=rc0, n=N0)

#Group 2:
lc1<-c(0,dat$V1[dat$V3==1])
rc1<-c(0,dat$V2[dat$V3==1])
N1<-length(lc1)-1
storage.mode(lc1)<-"double"
storage.mode(rc1)<-"double"
storage.mode(N1)<-"integer"
ret1<-.C("IcmNPMLE", x=lc1, S=rc1, n=N1)

x0<-ret0$x
S0<-ret0$S
n0<-ret0$n
storage.mode(x0)<-"double"
storage.mode(S0)<-"double"
storage.mode(n0)<-"integer"
x1<-ret1$x
S1<-ret1$S
n1<-ret1$n
storage.mode(x1)<-"double"
storage.mode(S1)<-"double"
storage.mode(n1)<-"integer"

for(m in 1:M){
  #impute interval-censored:
  iret0<-.C("impute",x0, S0, n0, y=lc0, d=rc0, N0)
  iret1<-.C("impute",x1, S1, n1, y=lc1, d=rc1, N1)
 
  #comparing 2 imputed samples:
  group<-c(rep(0,N0), rep(1,N1))
  G<-survdiff(Surv(c(iret0$y[-1], iret1$y[-1]),
                   c(iret0$d[-1], iret1$d[-1]))~group, rho=0)
  G0[m]<-sqrt(G$chisq)
  g0[m]<-G$exp-G$obs
  v0[m]<-g0[m]*g0[m]/G$chisq
  G<-survdiff(Surv(c(iret0$y[-1], iret1$y[-1]),
                   c(iret0$d[-1], iret1$d[-1]))~group, rho=1)
  G1[m]<-sqrt(G$chisq)
  g1[m]<-G$exp-G$obs
  v1[m]<-g1[m]*g1[m]/G$chisq
  G<-survdiff(Surv(c(iret0$y[-1], iret1$y[-1]),
                   c(iret0$d[-1], iret1$d[-1]))~group, rho=0.5)
  G12[m]<-sqrt(G$chisq)
  g12[m]<-G$exp-G$obs
  v12[m]<-g12[m]*g12[m]/G$chisq
  }
G0bar<-mean(G0)/sqrt(1+(1+1/M)*var(G0))
G1bar<-mean(G1)/sqrt(1+(1+1/M)*var(G1))
G12bar<-mean(G12)/sqrt(1+(1+1/M)*var(G12))
g0bar<-mean(g0)^2/(mean(v0)+(1+1/M)*var(g0))
g1bar<-mean(g1)^2/(mean(v1)+(1+1/M)*var(g1))
g12bar<-mean(g12)^2/(mean(v12)+(1+1/M)*var(g12))
sink(foutname, append=T)
cat(G0bar^2, G1bar^2, G12bar^2, g0bar, g1bar, g12bar,
    mean(v0)/var(g0), mean(v1)/var(g1), mean(v12)/var(g12), "\n")
sink()
}