##############################################
#
# HOW TO USE THIS FILE
# This file (ex6det.R) can just be sourced
# into the R program. But it is really intended to be used interactively
# 
#
######### ANY TEXT AFTER THE # CHARACTER ARE TREATED AS COMMENTS
#
# data stored in c:/GR/aprojects/peas/web/exemp6/data/ex6.Rdata
# this file is c:/GR/aprojects/peas/web/exemp3/program_code/ex3_R.html
#
##############################################
#
# data are in a data frame data
# to get the names of the variables
#
#
#
######################################################################
# mice would not work for this example not quite 
# sure why ? to do with too many categories
#
#  and NORM is going wrong too
#
######################################################################

library(norm)
#
#
rngseed(100)
datam<-as.matrix(data[,-1])  # exclude CASEID
 prelim.datam<-prelim.norm(datam)
 theta.datam<-em.norm(prelim.datam)
 imputed.datam<-imp.norm(prelim.datam,theta.datam,datam)

#
# now set volume variety to zero if either <=0
#

for (i in 1:6)  { volltz<-  (imputed.datam[,18+i]<=0)
                  varltz<-  (imputed.datam[,12+i]<=0)
                  anyltz<- volltz | varltz
                  imputed.datam[anyltz, c(12,18)+i]<-0 
                  imputed.datam<-cbind(imputed.datam,1-anyltz)               
}
dimnames(imputed.datam)[[2]][28:33]<-paste('DPREV',1:6,sep="")
pairs(data[,14:19+6]) 
#
# now get means prevalence by gender and  standard error for this imputation
#
imp.boys<-imputed.datam[data$GENDER==2,]
imp.girls<-imputed.datam[data$GENDER==1,]

mean.g<-list(apply(imp.girls,2,mean)[c(13:24,28:33)])
sd.g<-list(sqrt(apply(imp.girls,2,var)/apply(imp.girls,2,length))[c(13:24,28:33)])
mean.b<-list(apply(imp.boys,2,mean)[c(13:24,28:33)])
sd.b<-list(sqrt(apply(imp.boys,2,var)/apply(imp.boys,2,length))[c(13:24,28:33)])
#
#  now next imputations
#
for ( ii in 2:10){
 theta.datam<-em.norm(prelim.datam)
 imputed.datam<-imp.norm(prelim.datam,theta.datam,datam)
#
# now set volume variety to zero if either <=0
#
for (i in 1:6)  { volltz<-  (imputed.datam[,18+i]<=0)
                  varltz<-  (imputed.datam[,12+i]<=0)
                  anyltz<- volltz | varltz
                  imputed.datam[anyltz, c(12,18)+i]<-0 
                  imputed.datam<-cbind(imputed.datam,1-anyltz)               
}
dimnames(imputed.datam)[[2]][28:33]<-paste('DPREV',1:6,sep="")
#
# now get means prevalence by gender and  standard error for this imputation
#
imp.boys<-imputed.datam[data$GENDER==2,]
imp.girls<-imputed.datam[data$GENDER==1,]

mean.g[[ii]]<- apply(imp.girls,2,mean)[c(13:24,28:33)]
sd.g[[ii]]<- sqrt(apply(imp.girls,2,var)/apply(imp.girls,2,length))[c(13:24,28:33)]
mean.b[[ii]]<-apply(imp.boys,2,mean)[c(13:24,28:33)]
sd.b[[ii]]<- sqrt(apply(imp.boys,2,var)/apply(imp.boys,2,length))[c(13:24,28:33)]
}

#print(mean.g)
#print(sd.g)

result<-mi.inference(mean.b,sd.b)
print("boys")
ans<-cbind(result[[1]],result[[2]],result[[8]]*100)
dimnames(ans)[[2]]<-c('mean','se','% info lost')
print(ans)


result<-mi.inference(mean.g,sd.g)
print("girls")
ans<-cbind(result[[1]],result[[2]],result[[8]]*100)
dimnames(ans)[[2]]<-c('mean','se','% info lost')
print(ans)