R/eefAnalytics_03_2017.r

Defines functions caceSRTBoot.formula caceSRTBoot.default caceSRTBoot rbd.cace caceMSTBoot.formula caceMSTBoot.default caceMSTBoot crt.cace caceCRTBoot.formula caceCRTBoot.default caceCRTBoot srt.srt srt g.total g.within rbd.rbd rbdP rbd crt.crt mst.perm crt.perm crtP crt bootCompile errantSummary esProb schSummary esSummary betaSummary covSummary erantBAYES mlmBayes.formula mlmBayes.default mlmBayes mstFREQ.formula mstFREQ.default mstFREQ crtFREQ.formula crtFREQ.default crtFREQ srtFREQ.formula srtFREQ.default srtFREQ

Documented in caceCRTBoot caceMSTBoot caceSRTBoot crtFREQ mlmBayes mstFREQ srtFREQ

#' Multisite trial data.
#'
#' A multisite trial dataset containing 54 schools.
#'
#' \itemize{
#'   \item Posttest. posttest scores
#'   \item Prettest. prettest scores
#'   \item Intervention. indicator for intervention groups in two arm trial 
#'   \item Intervention2. a simulated indicator for intervention groups three arm trial 
#'   \item Compliance. percentage of sessions attended by pupils
#'coded as 1 for intervention group and 0 for control group.
#'   \item Compliance. percentage of sessions attended by pupils
#'   \item School. numeric school identifier
#' }
#'
#' @format A data frame with 210 rows and 5 variables
#' @name mstData
NULL

###iwq

#' Cluster randomised trial data.
#'
#' A cluster randomised trial dataset containing 22 schools.
#'
#' \itemize{
#'   \item Posttest. posttest scores
#'   \item Prettest. prettest scores
#'   \item Intervention. indicator for intervention groups 
#'   \item Intervention2. a simulated indicator for intervention groups three arm trial 
#'coded as 1 for intervention group and 0 for control group.
#'   \item Compliance. percentage of sessions attended by pupils
#'   \item School. numeric school identifier
#' }
#'
#' @format A data frame with 265 rows and 5 variables
#' @name crtData
NULL


## IMPORTS ##
#' @importFrom lme4 lmer ranef
#' @importFrom geoR rinvchisq
#' @importFrom mvtnorm rmvnorm
#' @importFrom metafor forest
#' @importFrom graphics abline barplot hist legend lines par plot points text
#' @importFrom stats confint lm model.frame model.matrix model.response quantile rnorm
NULL


#############################################################################
############# SRT main functions ################################################


#' Analysis of Simple Randomised Education Trial using Linear Regression Model.
#' 
#' \code{srtFREQ} perfoms analysis of educational trials under the assumption of independent errors among pupils. 
#' This can also be used with schools as fixed effects.
#' 
#' @export
#' @param formula the model to be analysed. It is of the form y~x1+x2+.... Where y is the outcome variable and Xs are the predictors.
#' @param intervention a string variable specifying the "intervention variable" as appeared in the formula. See example below
#' @param nBoot number of bootstraps required to generate bootstrap confidence interval. Default is NULL.
#' @param nPerm number of permutations required to generate permutated p-value.  Default is NULL.
#' @param data data frame containing the data to be analysed. 
#' @return S3 object; a list consisting of
#' \itemize{
#' \item \code{Beta}. Estimates and confidence intervals for the predictors specified in the model.
#' \item \code{ES}. Hedges' g effect size for the intervention(s). If nBoot is not specified, the confidence intervals are 95% CIs based on standard errors. If nBoot is specified, they are non-parametric bootstrapped confidence intervals.
#' \item \code{sigma2}. Residual variance.
#' \item \code{Perm}. A vector containing permutated effect sizes under null hypothesis. It is produced only if \code{nPerm} is specified.
#' \item \code{Bootstrap}. A vector containing bootstrapped effect sizes. It is prduced only if \code{nBoot} is specified.
#' }
#' @example inst/examples/srtExample.R
srtFREQ <- function(formula,intervention,nBoot=NULL,nPerm=NULL,data)UseMethod("srtFREQ")

#' @export
srtFREQ.default <- function(formula,intervention,nBoot=NULL,nPerm=NULL,data){stop("No correct formula input given.")}

#' @export
srtFREQ.formula <- function(formula,intervention,nBoot=NULL,nPerm=NULL,data=data){
	
	if(!is.null(nPerm) & !is.null(nBoot)){stop("Either nPerm or nBoot must be specified")}
	if(is.null(nPerm)){nPerm <-0}
	if(is.null(nBoot)){nBoot <-0}
	tmp3 <- which(colnames(data)==intervention)
	data[,tmp3] <- as.factor(data[,tmp3])
	mf <- model.frame(formula=formula, data=data)
	fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
	tmp <- colnames(fixedDesignMatrix )
	tmp[1]  <- "Intercept"
	colnames(fixedDesignMatrix)<- tmp
	posttest <- model.response(mf)
      intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]

	if(length(tmp3)!= 1){stop("Intervention variable misspecified")}

	output <- srt(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt)	


	if(nPerm>0){

		if(nPerm<1000){stop("nPerm must be greater than 1000")}

		permES<- matrix(NA,nPerm,(length(unique(trt))-1))
		
		set.seed(1020252)
		for (i in 1:nPerm){
				
			data[,which(colnames(data)==intervention)]<- sample(trt)
			fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))

			p2CRTFREQ <-srt(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt)	

			permES[i,]  <-  p2CRTFREQ$ES[,1]
		}
		output$Perm<- data.frame(permES)
	}

	if(nBoot > 0){

		if(nBoot<1000){stop("nBoot must be greater than 1000")}


		tid <- c(1:nrow(fixedDesignMatrix))

		set.seed(1020252)
		bootSamples <- sapply(c(1:nBoot),function(x)sample(tid,replace=TRUE))

		bootResults <- data.frame(apply(bootSamples ,2,function(bt)srt.srt(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,bt=bt)))
		bootResults2 <- t(matrix(unlist(bootResults ),(length(unique(trt))-1),nBoot))
		bootES <- apply(bootResults2,2,function(x)quantile(x,prob=c(0.025,0.975),na.rm=TRUE))
		bootES2 <- t(rbind(output$ES[,1],bootES ))
		row.names(bootES2)<- row.names(output$ES)
		colnames(bootES2)<- colnames(output$ES)
		colnames(bootResults2)<- row.names(output$ES)
		output$ES <-bootES2
		output$Bootstrap <- bootResults2 

	}
	output$Method <- "LM"
	class(output) <- "eefAnalytics"
	return(output) 
}



#############################################################################
############# CRT main functions ################################################
#' Analysis of Cluster Randomised Education Trials using Multilevel Model.
#' 
#' \code{crtFREQ} performs Analysis of cluster randomised education trial using multilevel model under the frequentist framework.
#' 
#' @export
#' @param formula the model to be analysed. It is of the form y ~ x1+x2+.... Where y is the outcome variable and Xs are the predictors.
#' @param random a string variable specifying the "clustering variable" as contained in the data. See example below
#' @param intervention a string variable specifying the "intervention variable" as appeared in the formula. See example below
#' @param nBoot number of bootstraps required to generate bootstrap confidence interval. Default is NULL.
#' @param nPerm number of permutations required to generate permutated p-value.  Default is NULL.
#' @param data data frame containing the data to be analysed.
#' @return S3 object; a list consisting of
#' \itemize{
#' \item \code{Beta}. Estimates and confidence intervals for preditors specified in the model.
#' \item \code{ES}. Hedges' g effect size for the intervention(s). If nBoot is not specified, the confidence intervals are 95% CIs based on standard errors. If nBoot is specified, they are non-parametric bootstrapped confidence intervals.
#' \item \code{covParm}. Vector of variance decomposition into between cluster variance (Schools) and within cluster variance (Pupils). It also contains the intral-cluster correlation (ICC).
#' \item \code{SchEffects}. Random intercepts for clusters, e.g schools.
#' \item \code{Perm}. A "nPerm x w" matrix containing permutated effect sizes using residual variance and total variance. "w" denotes number of intervention. "w=1" for two arm trial and "w=2" for three arm trial excluding the control group. It is produced only when \code{nPerm} is specified.
#' \item \code{Bootstrap}. A "w x nBoot" matrix containing the bootstrapped effect sizes using residual variance (Within) and total variance (Total). "w=1" for two arm trial and "w=2" for three arm trial excluding the control group. It is only prduced when \code{nBoot} is specified.
#' }
#' @example inst/examples/crtExample.R
crtFREQ<- function(formula,random,intervention,nPerm=NULL,nBoot=NULL,data)UseMethod("crtFREQ")

#' @export
crtFREQ.default<- function(formula,random,intervention,nPerm=NULL,nBoot=NULL,data){stop("No correct formula input given.")}


#' @export
crtFREQ.formula <- function(formula,random,intervention,nPerm=NULL,nBoot=NULL,data){
  
	data <- data[order(data[,which(colnames(data)==random)]),]
    
	intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]
	tmp2 <- which(colnames(data)==random)
	cluster2 <-  data[,tmp2]

	chk <- sum(rowSums(table(cluster2,trt)!=0)>1)
	if(chk >0){stop("This is not a CRT design")}
 	stp <- as.character(row.names(table(cluster2,trt)))
	stp2 <- as.numeric(apply(table(cluster2,trt),1,function(x)colnames(table(cluster2,trt))[x!=0]))

	if(!is.null(nPerm) & !is.null(nBoot)){stop("Either nPerm or nBoot must be specified")}
	if(is.null(nPerm)){nPerm <-0}
	if(is.null(nBoot)){nBoot <-0}
	tmp3 <- which(colnames(data)==intervention)
	data[,tmp3] <- as.factor(data[,tmp3])

	mf <- model.frame(formula=formula, data=data)
	mf <- mf[order(cluster2),]
	cluster <- cluster2[order(cluster2)]
	trt <- trt[order(cluster2)]
	fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
	tmp <- colnames(fixedDesignMatrix )
	tmp[1]  <- "Intercept"
	colnames(fixedDesignMatrix)<- tmp
	posttest <- model.response(mf)


	if(length(tmp2)!= 1){stop("Cluster variable misspecified")}
	if(length(tmp3)!= 1){stop("Intervention variable misspecified")}


	output <- crt(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,cluster=cluster)	


	if(nPerm>0){
		if(nPerm<999){stop("nPerm must be greater than 1000")}
		output$Perm<- crt.perm(formula,data,stp,stp2,intervention,cluster,nPerm,random)
		output$Perm <- round(data.frame(output$Perm),2)
	}


	if(nBoot >0){


		tid <- c(1:nrow(fixedDesignMatrix))

		set.seed(1020252)
		bootSamples <- NULL

		for(ii in 1:length(unique(cluster))){
			selID <- tid[cluster==unique(cluster)[ii]]
			if(length(selID)>0){
				selID2<- sapply(c(1:nBoot),function(x)sample(selID,length(selID),replace=TRUE))
				bootSamples <- rbind(bootSamples ,selID2)
			}
				
		}

		
		bootResults <- apply(bootSamples ,2,function(bt)crt.crt(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,cluster=cluster,bt=bt))
		bootES <- bootCompile(output=output,trt=trt,bootResults=bootResults)
		output$ES <- bootES
		output$Bootstrap <- round(data.frame(bootResults),2)  
	}
	output$Method <- "MLM"
	
	class(output) <- "eefAnalytics"
	return(output) 
}




#############################################################################
############# MST main functions ################################################

#############################################################################
############# MST main functions ################################################
#' Analysis of Multisite Randomised Education Trials using Multilevel Model.
#' 
#' \code{mstFREQ} performs analysis of multisite randomised education trial using multilevel model within the frequentist framework.
#' 
#' @export
#' @param formula the model to be analysed. It is of the form y ~ x1+x2+.... Where y is the outcome variable and Xs are the predictors.
#' @param random a string variable specifying the "clustering variable" as contained in the data. See example below
#' @param intervention a string variable specifying the "intervention variable" as appeared in the formula. See example below
#' @param nBoot number of bootstraps required to generate bootstrap confidence interval. Default is NULL.
#' @param nPerm number of permutations required to generate permutated p-value.  Default is NULL.
#' @param data data frame containing the data to be analysed.
#' @return S3 object; a list consisting of
#' \itemize{
#' \item \code{Beta}. Estimates and confidence intervals for preditors specified in the model.
#' \item \code{ES}. Hedges' g effect size for the intervention(s). If nBoot is not specified, the confidence intervals are 95% CIs based on standard errors. If nBoot is specified, they are non-parametric bootstrapped confidence intervals.
#' \item \code{covParm}. Vector of variance decomposition into between cluster variance (Schools), clustering by intervention interacttion (Intervention:School) and within cluster variance (Pupils). It also contains the intral-cluster correlation (ICC).
#' \item \code{SchEffects}. Random intercepts for clusters, e.g schools.
#' \item \code{Perm}. A "nPerm x w" matrix containing permutated effect sizes using residual variance and total variance. "w" denotes number of intervention. "w=1" for two arm trial and "w=2" for three arm trial excluding the control group. It is produced only when \code{nPerm} is specified.
#' \item \code{Bootstrap}. A "w x nBoot" matrix containing the bootstrapped effect sizes using residual variance (Within) and total variance (Total). "w=1" for two arm trial and "w=2" for three arm trial excluding the control group. It is only prduced when \code{nBoot} is specified.
#' }
#' @example inst/examples/mstExample.R
mstFREQ<- function(formula,random,intervention,nPerm=NULL,data,nBoot=NULL)UseMethod("mstFREQ")

#'@export
mstFREQ.default<- function(formula,random,intervention,nPerm=NULL,data,nBoot=NULL){stop("No correct formula input given.")}


#' @export
mstFREQ.formula <- function(formula,random,intervention,nPerm=NULL,data,nBoot=NULL){

	data <- data[order(data[,which(colnames(data)==random)],data[,which(colnames(data)==intervention)]),]
	trt <- data[,which(colnames(data)==intervention)]
	tmp2 <- which(colnames(data)==random)
	cluster2 = data[,tmp2]

	chk <- sum(rowSums(table(cluster2,trt)!=0)>1)
	if(chk ==0){stop("This is not a MST design")}
 	
	if(!is.null(nPerm) & !is.null(nBoot)){stop("Either nPerm or nBoot must be specified")}
	if(is.null(nPerm)){nPerm <-0}
	if(is.null(nBoot)){nBoot <-0}


	tmp3 <- which(colnames(data)==intervention)
	data[,tmp3] <- as.factor(data[,tmp3])
	mf <- model.frame(formula=formula, data=data)
	mf <- mf[order(cluster2),]
	cluster <- cluster2[order(cluster2)]
	trt <- trt[order(cluster2)]
	fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
	tmp <- colnames(fixedDesignMatrix )
	tmp[1]  <- "Intercept"
	colnames(fixedDesignMatrix)<- tmp
	posttest <- model.response(mf)
      intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]

	if(length(tmp2)!= 1){stop("Cluster variable misspecified")}
	if(length(tmp3)!= 1){stop("Intervention variable misspecified")}


	output <- rbd(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt,cluster=cluster)	

	if(nPerm > 0){

		if(nPerm<999){stop("nPerm must be greater than 1000")}
		output$Perm <- mst.perm(formula,data,trt,intervention,nPerm,random,cluster)
		output$Perm <- round(data.frame(output$Perm),2)
	}

	if(nBoot>0){


		tid <- c(1:nrow(fixedDesignMatrix))

		set.seed(1020252)
		bootSamples <- NULL

		for(ii in 1:length(unique(cluster))){
			selID <- tid[cluster==unique(cluster)[ii]]
			if(length(selID)>0){
				selID2<- sapply(c(1:nBoot),function(x)sample(selID,length(selID),replace=TRUE))
				bootSamples <- rbind(bootSamples ,selID2)
			}
				
		}
	
		bootResults <- apply(bootSamples ,2,function(bt)rbd.rbd(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt,cluster=cluster,bt=bt))

		bootES <- bootCompile(output=output,trt=trt,bootResults=bootResults)
		output$ES <- bootES
		output$Bootstrap <- round(data.frame(bootResults),2)  

	}
	output$Method <- "MLM"
	class(output) <- "eefAnalytics"
	return(output) 
}

###################################################################################################
############# Bayesian Multilevel Analysis of Randomised Education Trials
##################################################################################################

#' Bayesia multilevel analysis of randomised educatuon trials using vague priors.
#' 
#' \code{mlmBayes} performs analysis of randomised eduation trials using multilevel model under Bayesian framework 
#' assuming vague priors.
#' 
#' @export
#' @param formula the model to be analysed. It is of the form y ~ x1+x2+.... Where y is the outcome variable and Xs are the predictors.
#' @param random a string variable specifying the "clustering variable" as contained in the data. See example below
#' @param intervention a string variable specifying the "intervention variable" as appeared in the formula. See example below
#' @param nSim number of MCMC iterations. A minimum of 10,000 is recommended.
#'
#' @param data data frame containing the data to be analysed.
#' @return S3 object; a list consisting of
#' \itemize{
#' \item \code{Beta}. Estimates and confidence intervals for preditors specified in the model.
#' \item \code{ES}. Hedges' g effect size for the intervention(s). If nBoot is not specified, the confidence intervals are 95% CIs based on standard errors. If nBoot is specified, they are non-parametric bootstrapped confidence intervals.
#' \item \code{covParm}. Vector of variance decomposition into between cluster variance (Schools) and within cluster variance (Pupils). It also contains the intral-cluster correlation (ICC).
#' \item \code{SchEffects}. Random intercepts for clusters, e.g schools.
#' \item \code{ProbES}. A maxtrix containing the probability of observing effect size greater than a pre-specified threshold. 
#' }
#' @example inst/examples/bayesExample.R
mlmBayes <- function(formula,random,intervention,nSim,data)UseMethod("mlmBayes")

#' @export
mlmBayes.default <- function(formula,random,intervention,nSim=nSim,data){stop("No correct formula input given.")}

#' @export
mlmBayes.formula <- function(formula,random,intervention,nSim=nSim,data){
	data <- data[order(data[,which(colnames(data)==random)],data[,which(colnames(data)==intervention)]),]
	tmp3 <- which(colnames(data)==intervention)
	data[,tmp3] <- as.factor(data[,tmp3])
	tmp2 <- which(colnames(data)==random)
	cluster2 = data[,tmp2]

	mf <- model.frame(formula=formula, data=data)
	mf <- mf[order(cluster2),]
	fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
	tmp <- colnames(fixedDesignMatrix )
	tmp[1]  <- "Intercept"
	colnames(fixedDesignMatrix)<- tmp
	posttest <- model.response(mf)
      intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]
	tmp2 <- which(colnames(data)==random)
	cluster <- cluster2[order(cluster2)]
	nsim=nSim
	if(length(tmp2)!= 1){stop("Cluster variable misspecified")}
	if(length(tmp3)!= 1){stop("Intervention variable misspecified")}

	if(nsim < 10000){stop("nsim >= 10000 is recommended")}

	BayesOutput <- erantBAYES(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,cluster=cluster, nsim=nsim)
	output  <- errantSummary(bayesObject=BayesOutput,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention)
	output$SchEffects <- data.frame(Schools=unique(cluster),output$SchEffects)
	colnames(output$SchEffects ) <- c("Schools","Estimate","95% LB", "95% UB")
	output$Method <- "MLM"
	class(output) <- "eefAnalytics"
	return(output) 
}




#############################################################################
############# Bayesian functions ################################################

## Gibbs sampling - internal
erantBAYES<- function(posttest,fixedDesignMatrix,intervention,cluster, nsim) {
      
	cluster <- as.factor(cluster)   
	randomDesignMatrix <- as.matrix(as.data.frame(lm(posttest~cluster-1,x=TRUE)$x))
	nschools <- ncol(randomDesignMatrix )
	npar <- ncol(fixedDesignMatrix)
	totalObservations <- nrow(randomDesignMatrix )
	
	lgrp <- nchar(paste(intervention))
	sgrp <- substring(colnames(fixedDesignMatrix),1,lgrp)
	ssgrp <- colnames(fixedDesignMatrix)[sgrp==intervention]
	iindex <- which(sgrp==intervention)	

	ncluster <- tapply(cluster,cluster,length)

	###########
  	# Initial values  #
  	###########
	beta <- summary(lm(posttest~fixedDesignMatrix))$coef[,1]
	sigma2.ui <- 1
	b <- rnorm(nschools,0,1)
	bmat<-matrix(b,ncol=1,byrow=T)

  	#################
  	# Output matrix #
  	#################

  	nused <- floor(nsim/2)
  	nburnin <-  nsim-nused
  	outputbik <- matrix(0,nused,nschools)
  	outputCov <- matrix(0,nused,4)
      outputES <- matrix(0,nused,3*length(ssgrp) )
      outputBeta <- matrix(0,nused,npar)

	for (i in 1:nsim){

    		set.seed(123*i+i)


		# Update sigma2.e
		
		Ve <- totalObservations
		s2e0 <- fixedDesignMatrix %*% t(t(beta))
            s2e1 <-  t(matrix(bmat,nrow=nschools,ncol=totalObservations))
   		s2e2 <- randomDesignMatrix*s2e1 
           	s2e3 <- t(t(posttest))-s2e0 -rowSums(s2e2) 
		s2e <- (t(s2e3)%*%s2e3)/Ve
   		sigma2.e<-rinvchisq(1, Ve, s2e)


 		# Update sigma2.ui
		Vui <- nschools 
		s2ui <- (t(bmat)%*%bmat)/Vui
            sigma2.ui<-rinvchisq(1,Vui,s2ui)


 		# Update b
		usim1 <- t(t(posttest))-s2e0
  		usim2<- tapply(usim1,cluster,mean)
		uV0 <- (ncluster +  as.numeric(sigma2.e/sigma2.ui))^-1
		uV <-  diag(as.numeric(sigma2.e)*uV0 )
            ui <- ncluster*uV0*usim2
   		b <-rmvnorm(1,ui,uV)
   		bmat<-matrix(b,ncol=1,byrow=T)

		# Update Beta 
            beta1 <-solve(t(fixedDesignMatrix) %*% fixedDesignMatrix )
		beta2 <-  t(matrix(bmat,nrow=nschools,ncol=totalObservations))
   		beta3 <- randomDesignMatrix*beta2
		beta4 <- t(t(posttest)) -rowSums(beta3)
            beta5 <- t(fixedDesignMatrix)%*% beta4
		beta6 <- beta1 %*% beta5 
   		vbeta<- beta1*as.numeric(sigma2.e )
   		beta<-c(rmvnorm(1,beta6,vbeta))
            
            treatment=beta[iindex]


   		if(i > nburnin ){ 
  			outputbik[(i-nburnin),] <- as.numeric(b)
  			outputCov[(i-nburnin),] <- c(sigma2.e,sigma2.ui,sigma2.e+sigma2.ui,(sigma2.ui)/(sigma2.e+sigma2.ui))
  			outputBeta[(i-nburnin),] <- as.numeric(beta)
  			outputES[(i-nburnin),] <- c(sapply(treatment,function(x)x/sqrt(c(sigma2.e,(sigma2.ui),sigma2.e+sigma2.ui))))
		}
	}


      esF <- c("Within","Between","Total")

      es.names <- c(sapply(ssgrp,function(x)paste(x,esF,sep="") ))
      colnames(outputBeta)<- colnames(fixedDesignMatrix)
	colnames(outputCov) <- c("Within","Between","Total","ICC")
	colnames(outputES) <- es.names
	colnames(outputbik) <- unique(cluster )

 	Output <- list(randomEffects=outputbik,CovParameters=outputCov,Beta=outputBeta,ES=outputES)
  	return(Output)
}


## summarise covariance parameters - internal

covSummary <- function(bayesObject){

	covParm <- bayesObject$CovParameters
      covParm2 <- colMeans(covParm) # remove all roundings that appear in the middle of functions 
      covParm3 <- t(apply(covParm,2,function(x)quantile(x,prob=c(0.025,0.975))))
	covParm4 <- data.frame(cbind(covParm2,covParm3))
	covParm5 <- sqrt(covParm4)
	colnames(covParm4) <- c("Variance","95% LB","95% UB")
	rownames(covParm4) <- c("Pupils","Schools","Total","ICC")
	covParm5 <- covParm4[,1]
	return(covParm5)
}


## summarise beta parameters - internal
betaSummary <- function(bayesObject){

	betaParm <- bayesObject$Beta
      betaParm2 <- colMeans(betaParm)
      betaParm3 <- t(apply(betaParm,2,function(x)quantile(x,prob=c(0.025,0.975))))
	betaParm4 <- data.frame(cbind(betaParm2,betaParm3))
	colnames(betaParm4) <- c("Estimate","95% LB","95% UB")
	rownames(betaParm4) <- colnames(betaParm)
	return(betaParm4)
}

## summarise Effect sizes - internal
esSummary <- function(bayesObject,fixedDesignMatrix,intervention){

	esParm <- bayesObject$ES
      esParm2 <- colMeans(esParm)
      esParm3 <- t(apply(esParm,2,function(x)quantile(x,prob=c(0.025,0.975))))
	esParm4 <- data.frame(cbind(esParm2,esParm3))
	colnames(esParm4) <- c("Estimate","95% LB","95% UB")
	btp <- nrow(esParm4)
      if(btp <=3){return(round(esParm4,2))}

	if(btp >3){
		btp2 <- seq(3,btp,3)
		btp3 <- seq(1,btp,3)
		ouptut <- list()
		for(i in 1:length(btp2)){
			ouptut[[i]] <- round(esParm4[btp3[i]:btp2[i],],2)
		}
	rname <- colnames(fixedDesignMatrix)[substring(colnames(fixedDesignMatrix),1,nchar(intervention))==intervention]
	names(ouptut)<- rname
	return(ouptut)
	}
	
}

## summarise random effects - internal

schSummary <- function(bayesObject){

	schParm <- bayesObject$randomEffects
      schParm2 <- colMeans(schParm)
      schParm3 <- t(apply(schParm,2,function(x)quantile(x,prob=c(0.025,0.975))))
	schParm4 <- data.frame(cbind(schParm2,schParm3))
	colnames(schParm4) <- c("Estimate","95% LB","95% UB")
	schParm5<- schParm4
	return(schParm5)

}


##summarise minimum expected effect size - internal
esProb <- function(bayesObject,esOutput){
      es <- seq(0,1,0.1)
	esParm <- bayesObject$ES 
      esParm2 <- sapply(es,function(x)colMeans(esParm>=x))
	esParm4 <- data.frame(cbind(ES=es,t(esParm2)))
	rownames(esParm4) <- NULL
	return(esParm4)
	
}

##summarise all bayesian parameters - internal
errantSummary <- function(bayesObject,fixedDesignMatrix,intervention){
	covValues <- covSummary(bayesObject=bayesObject)
	covValues <- data.frame(covValues[c(2,1,3,4)])
	row.names(covValues ) <- c("Schools","Pupils","Total","ICC")
	covValues <- t(covValues )
	row.names(covValues ) <- NULL
	betaValues <- betaSummary(bayesObject=bayesObject)
	esValues <- esSummary(bayesObject,fixedDesignMatrix,intervention)
	schValues <-schSummary(bayesObject=bayesObject)
	es.prob <- esProb(bayesObject=bayesObject,esOutput=esValues)
      output <- list(Beta=round(betaValues,2),covParm= round(covValues,2),ES=esValues,ProbES=round(data.frame(es.prob),2),SchEffects=round(schValues,2))
}




## compile bootstrap results - internal
bootCompile <- function(output,trt,bootResults){

	withinBoot <- matrix(NA,nrow=length(bootResults),ncol=(length(unique(trt ))-1))
	totalBoot <- matrix(NA,nrow=length(bootResults),ncol=(length(unique(trt ))-1))
	for(k in 1:length(bootResults)){
		tmp <- bootResults[[k]]
		tmpR <- NULL
		for(j in 1:length(tmp)){
					
			tmpR  <- c(tmpR,tmp[[j]][,1])
		}
		withinBoot[k,] <- tmpR[seq(1,2*length(tmp),2)]
		totalBoot[k,] <- tmpR[seq(2,2*length(tmp),2)]
	}

	withinCI <- apply(withinBoot,2,function(x)quantile(x,prob=c(0.025,0.975)))
	TotalCI <- apply(totalBoot,2,function(x)quantile(x,prob=c(0.025,0.975)))
	tmpES <- list()
	for(kk in 1:length(output$ES)){
		tmp1 <- rbind(withinCI[,kk], TotalCI[,kk])
		tmp2 <- cbind(output$ES[[kk]][,1],tmp1)
		colnames(tmp2 )<- c("Estimate","95% LB","95% UB")
		row.names(tmp2)	<- c("Within","Total")
		tmpES[[kk]] <- round(tmp2,2)
	}
	names(tmpES) <- names(output$ES)
	return(tmpES)

}

## random intercept model - internal
crt <- function(posttest,fixedDesignMatrix,intervention,cluster){
	
	freqFit <- lmer(posttest~ fixedDesignMatrix-1+ (1|cluster))
		
      np<- row.names(summary(freqFit)$coef)
      cit <- confint(freqFit,np)
	betaB <- data.frame(cbind(summary(freqFit)$coefficients[,1],cit))
	row.names(betaB)<- colnames(fixedDesignMatrix)
	colnames(betaB) <- c("Estimate","95% LB ","95% UB")
	betaB <- betaB
	var.B<- as.numeric(summary(freqFit)$varcor)
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.B,var.W,var.B+var.W,(var.B/(var.B+var.W)))
	sigmaBE <- sigmaBE
	names(sigmaBE)<- c("Schools","Pupils","Total","ICC")
      schRand <- data.frame(unique(cluster),ranef(freqFit)$cluster)
	names(schRand)<- c("Schools","Estimate")
	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- g.within(var.w=var.W, beta=beta, icc=ICC , group=group, schoolID=cluster)
  		esTotal <- g.total(var.tt=var.tt, beta=beta, icc=ICC , group=group, schoolID=cluster)
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		colnames(output1) <- c("Estimate","95% LB","95% UB")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- round(output1,2)
	}
	names(output2) <- row.names(betaB)[btp]
      output <- list(Beta=round(betaB,2),covParm=round(sigmaBE,2),ES=output2,SchEffects=round(schRand,2))

	return(output)

}

## internal
crtP <- function(posttest,fixedDesignMatrix,intervention,cluster){
	
	freqFit <- lmer(posttest~ fixedDesignMatrix-1+ (1|cluster))
		
	betaB <- data.frame(summary(freqFit)$coefficients[,1])
	row.names(betaB)<- colnames(fixedDesignMatrix)
	#colnames(betaB) <- c("Estimate","95% LB ","95% UB")
	betaB <- betaB
	var.B<- as.numeric(summary(freqFit)$varcor)
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.B,var.W)
	sigmaBE <- sigmaBE
	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- g.within(var.w=var.W, beta=beta, icc=ICC , group=group, schoolID=cluster)
  		esTotal <- g.total(var.tt=var.tt, beta=beta, icc=ICC , group=group, schoolID=cluster)
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		colnames(output1) <- c("Estimate","95% LB","95% UB")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- round(output1,2)
	}
	names(output2) <- row.names(betaB)[btp]
      output <- list(ES=output2)

	return(output)

}


## - internal

crt.perm <- function(formula,data,stp,stp2,intervention,cluster,nPerm,random){

		data2 <- data[,-which(colnames(data)==intervention)]
		
		g <- matrix(NA,nPerm,2*(length(unique(stp2))-1))

		for(i in 1:nPerm){
			set.seed(12890*i+1)
			tp3 <- data.frame(stp,sample(stp2))
			names(tp3) <- c(paste(random),paste(intervention))
			data.tp4 <- merge(data2,tp3,by=random)
			data.tp4 <- data.tp4[order(data.tp4[,which(colnames(data.tp4)==random)]),]
			cluster = data.tp4[,which(colnames(data.tp4)==random)]
			tmp34 <- which(colnames(data.tp4)==intervention)
			data.tp4[,tmp34] <- as.factor(data.tp4[,tmp34])
			mf <- model.frame(formula=formula, data=data.tp4)
			fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data.tp4)))
			tmp <- colnames(fixedDesignMatrix )
			tmp[1]  <- "Intercept"
			colnames(fixedDesignMatrix)<- tmp
			posttest <- model.response(mf)
      		intervention <- intervention
				
			p2CRTFREQ <-crtP(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,cluster=cluster)	
			chkppp <- unlist(p2CRTFREQ$ES)
			chkppp2 <- c(seq(1,6*(length(unique(tp3[,2]))-1),6),seq(2,6*(length(unique(tp3[,2]))-1),6))
			chkppp3 <- chkppp2[order(chkppp2)] 
			g[i,]  <-  chkppp[chkppp3]
		}
			ntpp <- rep(names(p2CRTFREQ$ES),2)
			ntpp <- ntpp[order(ntpp )]
			wt <- rep(c("Within","Total"),length(names(p2CRTFREQ$ES)))
			colnames(g) <- paste(ntpp ,wt,sep="")
			return(g)
}


## - internal


mst.perm <- function(formula,data,trt,intervention,nPerm,random,cluster){

		data2 <- data
		g <- matrix(NA,nPerm,2*(length(unique(trt))-1))
		for(i in 1:nPerm){
			set.seed(12890*i+1)
			data2[,which(colnames(data)==intervention)]<-unlist(tapply(trt,cluster,function(x)sample(x)))
			data3 <- data2[order(data2[,which(colnames(data2)==random)],data2[,which(colnames(data2)==intervention)]),]
			cluster = data3[,which(colnames(data3)==random)]
			mf <- model.frame(formula=formula, data=data3)
			fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data3)))
			tmp <- colnames(fixedDesignMatrix )
			tmp[1]  <- "Intercept"
			colnames(fixedDesignMatrix)<- tmp
			posttest <- model.response(mf)
      		intervention <- intervention
			trt2 <- data3[,which(colnames(data3)==intervention)]	
			p2CRTFREQ <-rbdP(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt2,cluster=cluster)	

			chkppp <- unlist(p2CRTFREQ$ES)
			chkppp2 <- c(seq(1,6*(length(unique(trt))-1),6),seq(2,6*(length(unique(trt))-1),6))
			chkppp3 <- chkppp2[order(chkppp2)] 
			g[i,]  <-  chkppp[chkppp3]
		}
			ntpp <- rep(names(p2CRTFREQ$ES),2)
			ntpp <- ntpp[order(ntpp )]
			wt <- rep(c("Within","Total"),length(names(p2CRTFREQ$ES)))
			colnames(g) <- paste(ntpp ,wt,sep="")
			return(g)

}



## - internal
crt.crt<- function(posttest,fixedDesignMatrix,intervention,cluster,bt){
	
	posttest2 <- posttest[bt]
	fixedDesignMatrix2 <- fixedDesignMatrix[bt,]
	cluster2 <- cluster[bt]


	freqFit <- try(lmer(posttest2~ fixedDesignMatrix2-1+(1|cluster2)),silent=TRUE)
	output2 <- NULL
     if(attr(freqFit,"class")!="try-error"){

	betaB <- data.frame(summary(freqFit)$coefficients[,1])
	row.names(betaB)<- colnames(fixedDesignMatrix)
	betaB <- betaB
	var.B2<- as.matrix(summary(freqFit)$varcor)
	var.B3 <- c(matrix(attr(var.B2[[1]],"stddev")))
	var.B <- var.B3^2
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.B,var.W)
	sigmaBE <- sigmaBE
	names(sigmaBE)<- c("Schools","Pupils")

	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- beta/sqrt(var.W)
  		esTotal <- beta/sqrt(var.tt) 
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		names(output1) <- c("Estimate")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- output1
	}
	names(output2) <- row.names(betaB)[btp]
	
	}
	return(output2)
}

## - internal

rbd <- function(posttest,fixedDesignMatrix,intervention,trt,cluster){

	freqFit <- lmer(posttest~ fixedDesignMatrix-1+(1|trt:cluster)+(1|cluster))
      np<- row.names(summary(freqFit)$coef)
      cit <- confint(freqFit,np)
	betaB <- data.frame(cbind(summary(freqFit)$coefficients[,1],cit))
	row.names(betaB)<- colnames(fixedDesignMatrix)
	colnames(betaB) <- c("Estimate","95% LB ","95% UB")
	betaB <- betaB
	var.B2<- as.matrix(summary(freqFit)$varcor)
	var.B3 <- c(c(matrix(attr(var.B2[[1]],"stddev"))),c(matrix(attr(var.B2[[2]],"stddev"))))
	var.B3 <- var.B3^2
	var.sch <- var.B3[1]
	var.schTrt <- var.B3[2]	
	var.B <- sum(var.B3)
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.sch,var.schTrt,var.W,var.tt,ICC )
	sigmaBE <- sigmaBE
	names(sigmaBE)<- c("Schools","Intervention:School","Pupils","Total","ICC")
      schRand <- data.frame(unique(cluster),ranef(freqFit)$cluster)
	names(schRand)<- c("Schools","Estimate")
	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- g.within(var.w=var.W, beta=beta, icc=ICC , group=group, schoolID=cluster)
  		esTotal <- g.total(var.tt=var.tt, beta=beta, icc=ICC , group=group, schoolID=cluster)
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		colnames(output1) <- c("Estimate","95% LB","95% UB")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- round(output1,2)
	}
	names(output2) <- row.names(betaB)[btp]
      output <- list(Beta=round(betaB,2),covParm=round(sigmaBE,2),ES=output2,SchEffects=round(schRand,2))


	return(output)
}


## - internal

rbdP <- function(posttest,fixedDesignMatrix,intervention,trt,cluster){

	freqFit <- lmer(posttest~ fixedDesignMatrix-1+(1|trt:cluster)+(1|cluster))

	betaB <- data.frame(summary(freqFit)$coefficients)
	row.names(betaB)<- colnames(fixedDesignMatrix)
	#colnames(betaB) <- c("Estimate","95% LB ","95% UB")
	betaB <- betaB
	var.B2<- as.matrix(summary(freqFit)$varcor)
	var.B3 <- c(c(matrix(attr(var.B2[[1]],"stddev"))),c(matrix(attr(var.B2[[2]],"stddev"))))
	var.B3 <- var.B3^2
	var.sch <- var.B3[1]
	var.schTrt <- var.B3[2]	
	var.B <- sum(var.B3)
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.sch,var.schTrt,var.W)
	sigmaBE <- sigmaBE
	names(sigmaBE)<- c("Schools","Intervention:School","Pupils")
      schRand <- data.frame(ranef(freqFit)$cluster)
	names(schRand)<- "Estimate"
	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- g.within(var.w=var.W, beta=beta, icc=ICC , group=group, schoolID=cluster)
  		esTotal <- g.total(var.tt=var.tt, beta=beta, icc=ICC , group=group, schoolID=cluster)
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		colnames(output1) <- c("Estimate","95% LB","95% UB")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- round(output1,2)
	}
	names(output2) <- row.names(betaB)[btp]
      output <- list(ES=output2)


	return(output)
}


## - internal

rbd.rbd <- function(posttest,fixedDesignMatrix,intervention,trt,cluster,bt){
	
	posttest2 <- posttest[bt]
	fixedDesignMatrix2 <- fixedDesignMatrix[bt,]
	trt2 <- trt[bt]
	cluster2 <- cluster[bt]


	freqFit <- try(lmer(posttest2~ fixedDesignMatrix2-1+(1|trt2:cluster2)+(1|cluster2)),silent=TRUE)
	output2 <- NULL
     if(attr(freqFit,"class")!="try-error"){

	betaB <- data.frame(summary(freqFit)$coefficients[,1])
	row.names(betaB)<- colnames(fixedDesignMatrix)
	betaB <- betaB
	var.B2<- as.matrix(summary(freqFit)$varcor)
	var.B3 <- c(c(matrix(attr(var.B2[[1]],"stddev"))),c(matrix(attr(var.B2[[2]],"stddev"))))
	var.B3 <- var.B3^2
	var.sch <- var.B3[1]
	var.schTrt <- var.B3[2]	
	var.B <- sum(var.B3)
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.sch,var.schTrt,var.W)
	sigmaBE <- sigmaBE
	names(sigmaBE)<- c("Schools","Intervention:School","Pupils")

	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- beta/sqrt(var.W)
  		esTotal <- beta/sqrt(var.tt) 
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		names(output1) <- c("Estimate")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- output1
	}
	names(output2) <- row.names(betaB)[btp]
	
	}
	return(output2)
}


## - internal

g.within <- function(var.w, beta, icc, group, schoolID){
   t <- group; id <- schoolID
   d.w <- (beta/sqrt(var.w))
   n.it <- table(id[t==1]); n.ic <- table(id[t==0])
   m.t <- length(unique(id[t==1])); m.c <- length(unique(id[t==0]))
   M <- (m.t + m.c)
   N.t <- sum(table(id[t==1])); N.c <- sum(table(id[t==0]))
   N <- (N.t + N.c)
   n.sim.1 <- ((N.c * sum(n.it^2))/(as.numeric(N.t)*as.numeric(N)))
   n.sim.2 <- ((N.t * sum(n.ic^2))/(as.numeric(N.c)*as.numeric(N)))
   n.sim <- (n.sim.1 + n.sim.2)
   vterm1 <- ((N.t+N.c)/(as.numeric(N.t)*as.numeric(N.c)))
   vterm2 <- (((1+(n.sim-1)*icc))/(1-icc))
   vterm3 <- ((d.w^2)/(2*(N-M)))
   se <- sqrt(vterm1*vterm2+vterm3)
   LB <- (d.w-1.96*se); UB <- (d.w+1.96*se)
   output <- data.frame(d.w, LB, UB)
   names(output) <- c("g", "LB", "UB")
   return(output)
}

## - internal

g.total <- function(var.tt, beta, icc, group, schoolID){
   t <- group; id <- schoolID
   n.it <- table(id[t==1]); n.ic <- table(id[t==0])
   m.t <- length(unique(id[t==1])); m.c <- length(unique(id[t==0]))
   M <- (m.t + m.c)
   N.t <- sum(table(id[t==1])); N.c <- sum(table(id[t==0]))
   N <- (N.t + N.c)
   n.ut <- ((N.t^2-sum(n.it^2))/(as.numeric(N.t)*as.numeric(m.t-1)))
   n.uc <- ((N.c^2-sum(n.ic^2))/(as.numeric(N.c)*as.numeric(m.c-1)))
   dt.1 <- (beta/sqrt(var.tt))
   dt.2 <- sqrt(1-icc*(((N-n.ut*m.t-n.uc*m.c)+n.ut+n.uc-2)/(N-2)))
   d.t <- (dt.1*dt.2)
   
   n.sim.1 <- ((as.numeric(N.c) * sum(n.it^2))/(as.numeric(N.t)*as.numeric(N)))
   n.sim.2 <- ((as.numeric(N.t) * sum(n.ic^2))/(as.numeric(N.c)*as.numeric(N)))
   n.sim <- (n.sim.1 + n.sim.2)
   B <- (n.ut*(m.t-1)+n.uc*(m.c-1))
   A.t <- ((as.numeric(N.t)^2*sum(n.it^2)+(sum(n.it^2))^2-2*as.numeric(N.t)*sum(n.it^3))/as.numeric(N.t)^2)
   A.c <- ((as.numeric(N.c)^2*sum(n.ic^2)+(sum(n.ic^2))^2-2*as.numeric(N.c)*sum(n.ic^3))/as.numeric(N.c)^2)
   A <- (A.t + A.c)
   
   vterm1 <- (((N.t+N.c)/(as.numeric(N.t)*as.numeric(N.c)))*(1+(n.sim-1)*icc))
   vterm2 <- (((N-2)*(1-icc)^2+A*icc^2+2*B*icc*(1-icc))*d.t^2)
   vterm3 <- (2*(N-2)*((N-2)-icc*(N-2-B)))
   se <- sqrt(vterm1+vterm2/vterm3)
   LB <- (d.t-1.96*se); UB <- (d.t+1.96*se)
   output <- data.frame(d.t, LB, UB)
   names(output)<- c("g", "LB", "UB")
   return(output)
}


## - internal
srt <- function(posttest,fixedDesignMatrix,intervention,trt){
	
	freqFit <- lm(posttest~ fixedDesignMatrix-1)
      cit <- confint(freqFit)
	citt <- rowSums(is.na(cit))
	betaB <- data.frame(cbind(summary(freqFit)$coefficients[which(citt==0),1],cit[which(citt==0),]))
	row.names(betaB)<- colnames(fixedDesignMatrix)[which(citt==0)]
	colnames(betaB) <- c("Estimate","95% LB ","95% UB")
	betaB <- betaB

	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	btp2 <- which(substring(colnames(fixedDesignMatrix),1,nchar(intervention))==intervention)

	tmpTRT <- table(trt)      
      output2 <- NULL
	for( i in 1:length(btp )){
	
  		beta <- betaB[btp[i],1]
   		sd.pool <- summary(freqFit)$sigma
   		cd <- (beta/sd.pool)
             trt2 <- unique(fixedDesignMatrix[,btp2[i]])
		n.c <- tmpTRT[names(tmpTRT)==trt2[2]] 
		n.t <- tmpTRT[names(tmpTRT)==trt2[1]] 
   		var.cd <- ((n.t+n.c)/(n.t*n.c)+cd^2/(2*(n.t+n.c)))
   		se.cd <- sqrt(var.cd)
   		cd.lb <- (cd - 1.96*se.cd)
   		cd.ub <- (cd + 1.96*se.cd)
   		j.df <- (1 - (3/(4*(n.t+n.c-2)-1)))
   		g <- (j.df*cd)
   		var.g <- (j.df^2 * var.cd)
   		se.g <- sqrt(var.g)
   		g.lb <- (g - 1.96*se.g)
   		g.ub <- (g + 1.96*se.g)
  		gtmp <- round(data.frame(g, g.lb, g.ub),2)
		output2 <- rbind(output2,gtmp)
	}
	names(output2)<- c("Estimate","95% LB","95% UB")
	row.names(output2) <- row.names(betaB)[btp]
      output <- list(Beta=round(betaB,2),ES=round(output2,2),sigma2=round(summary(freqFit)$sigma^2,2))
	return(output)

}


## - internal
srt.srt<- function(posttest,fixedDesignMatrix,intervention,bt){
	
	posttest2 <- posttest[bt]
	fixedDesignMatrix2 <- fixedDesignMatrix[bt,]



	freqFit <- try(lm(posttest2~ fixedDesignMatrix2-1),silent=TRUE)
	output2 <- NULL
      if(attr(freqFit,"class")!="try-error"){

		betaB <- data.frame(summary(freqFit)$coefficients[,1])
            ntpp <- as.character(sapply(as.character(rownames(betaB)),function(x)substring(x,(nchar("fixedDesignMatrix2")+1),nchar(x))))
		row.names(betaB)<- ntpp
		betaB <- betaB

		sd.pool <- summary(freqFit)$sigma

		btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      		output1 <- NULL
		for( i in 1:length(btp )){
	
  			beta <- betaB[btp[i],1]	
     
  	
			output1[i] <-beta/sd.pool

		}
		names(output1) <- row.names(betaB)[btp]
	
	}
	output1<- matrix(output1,1,length(btp))

	return(output1)
}



#################
############# CRT main functions ################################################
#'Complier Average Causal Effect (CACE) Analysis of Cluster Randomised Education Trials using Multilevel Model.
#' 
#' \code{caceCRTBoot} performs exploratoty CACE analysis of cluster randomised education trials.
#' 
#' @export
#' @param formula the model to be analysed. It is of the form y ~ x1+x2+.... Where y is the outcome variable and Xs are the predictors.
#' @param random a string variable specifying the "clustering variable" as contained in the data. See example below
#' @param intervention a string variable specifying the "intervention variable" as appeared in the formula. See example below
#' @param compliance a string variable specifying the "compliance variable" as contained in the data.  The data must be in percentages ranging from 0 - 100.
#' @param nBoot number of bootstraps required to generate bootstrap confidence interval. Default is NULL.
#' @param data data frame containing the data to be analysed. 
#' @return S3 object; a list consisting of
#' \itemize{
#' \item \code{CACE}. Estimates of CACE adjusted effect sizes based on pre-specified thresholds. Only produced for threshold with at least 50% compliance rate.
#' \item \code{Compliers}. Percentage of pupils that achieved a pre-specified threshold of compliance.
#' }
#' @example inst/examples/crtCACEExample.R
caceCRTBoot<- function(formula,random,intervention,compliance,nBoot,data)UseMethod("caceCRTBoot")


#' @export
caceCRTBoot.default <- function(formula,random,intervention,compliance,nBoot,data){stop("No correct formula input given.")}

#' @export
caceCRTBoot.formula <- function(formula,random,intervention,compliance,nBoot,data){
  
	data <- data[order(data[,which(colnames(data)==random)]),]
       
	intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]

	if(length(table(trt))!=2){stop("Not applicable to multi-arm trials. Re-define the treatment group as dummy variables and run the function for each dummy variable separately")}
	comp1 <- data[,which(colnames(data)==compliance)]
	tmp2 <- which(colnames(data)==random)
	cluster2 = data[,tmp2]
	tmp3 <- which(colnames(data)==intervention)
	data[,tmp3] <- as.factor(data[,tmp3])
	mf <- model.frame(formula=formula, data=data)
	mf <- mf[order(cluster2),]
	comp2 <- comp1[order(cluster2)]
	cluster <- cluster2[order(cluster2)]
	fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
	tmp <- colnames(fixedDesignMatrix )
	tmp[1]  <- "Intercept"
	colnames(fixedDesignMatrix)<- tmp
	posttest <- model.response(mf)

	output <- crt.cace(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,cluster=cluster)	
	
	theta <- seq(0,ifelse(max(comp2)==100,90,max(comp2)),10)
	alpha_ITT_E <- output$ES[[1]][2,1]

	d <- sapply(theta,function(x)ifelse(comp2 > x,1,0))
	ibit <- data[,tmp3]
      tc.pp <- apply(d,2,function(x) table(x,ibit))
	tc.pp2 <- apply(d,2,function(x) nrow(table(x,ibit)))
	if(sum(tc.pp2==2) != length(tc.pp2)){
		t.tpp1 <- c.tpp1 <- p1 <- NULL
		for (j in 1:length(tc.pp2)){
			if(tc.pp2[j]==1){
				t.tpp1[j] <- 1
				c.tpp1[j] <- 0
				p1[j] <- 1
			}
			if(tc.pp2[j]==2){	
				tppk <- unlist(tc.pp[j])
      			c.tpp1[j] <- tppk[2]/(tppk[1]+tppk[2])
				t.tpp1[j] <-tppk[4]/(tppk[3]+tppk[4])
				p1[j] <- t.tpp1[j] - c.tpp1[j] 
			}
		}
	}

	if(sum(tc.pp2==2) == length(tc.pp2)){

      	c.tpp1 <- tc.pp[2,]/(tc.pp[1,]+tc.pp[2,])
		t.tpp1 <-tc.pp[4,]/(tc.pp[3,]+tc.pp[4,])
		p1 <- t.tpp1 - c.tpp1

	}

	caceES  <- alpha_ITT_E/(p1)
	p1Table <- rbind(t.tpp1,c.tpp1,p1)
	p1Table <- data.frame(round(p1Table,2))
	row.names(p1Table )<- c("pT","pC","P=PT-pC")
	colnames(p1Table ) <- paste("P > ", theta,sep="")
	

	p2Table <- p1Table[,round(p1,2)>0.5]
	cutTheta <- theta[round(p1,2)>0.5]
	caceES2 <- caceES[round(p1,2) > 0.5] 

	tid <- c(1:nrow(fixedDesignMatrix))
	set.seed(1020252)
	bootSamples <- NULL

	for(ii in 1:length(unique(cluster))){
		selID <- tid[cluster==unique(cluster)[ii]]
		if(length(selID)>0){
			selID2<- sapply(c(1:nBoot),function(x)sample(selID,length(selID),replace=TRUE))
			bootSamples <- rbind(bootSamples ,selID2)
		}
				
	}
	row.names(bootSamples ) <- NULL
	bootSamples <- data.frame(bootSamples )

	bESOutput <- NULL

	for(i in 1:ncol(bootSamples )){

		bData <- data[bootSamples[,i],]

		bcluster <- bData[,tmp2]
		bmf <- model.frame(formula=formula, data=bData)
		bfixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(bmf, "terms"), data=bData)))
		bposttest <- model.response(bmf)
		boutput <- crt.cace(posttest=bposttest,fixedDesignMatrix=bfixedDesignMatrix,intervention=intervention,cluster=bcluster)	
	
		balpha_ITT_E <- boutput$ES[[1]][2,1]

		bcomp <- bData[,which(colnames(bData)==compliance)]
		bd <- sapply(cutTheta,function(x)ifelse(bcomp > x,1,0))
		bibit <- bData[,tmp3]
      	btc.pp <- apply(bd,2,function(x) table(x,bibit))
		btc.pp2 <- apply(bd,2,function(x) nrow(table(x,bibit)))
		if(sum(btc.pp2==2) != length(btc.pp2)){
			t.tpp1 <- c.tpp1 <-bp1 <- NULL
			for (j in 1:length(btc.pp2)){
				if(btc.pp2[j]==1){
					t.tpp1[j] <- 1
					c.tpp1[j] <- 0
					bp1[j] <- 1
				}
				if(btc.pp2[j]==2){	
					tppk <- unlist(btc.pp[j])
      				c.tpp1[j] <- tppk[2]/(tppk[1]+tppk[2])
					t.tpp1[j] <-tppk[4]/(tppk[3]+tppk[4])
					bp1[j] <- t.tpp1[j] - c.tpp1[j] 
				}
			}
		}

		if(sum(btc.pp2==2) == length(btc.pp2)){

      		c.tpp1 <- btc.pp[2,]/(btc.pp[1,]+btc.pp[2,])
			t.tpp1 <-btc.pp[4,]/(btc.pp[3,]+btc.pp[4,])
			bp1 <- t.tpp1 - c.tpp1

		}
		bcaceES  <- balpha_ITT_E/(bp1)
		bESOutput <- rbind(bESOutput ,bcaceES)  

	}

	bootES <- apply(bESOutput ,2,function(x)quantile(x,prob=c(0.025,0.975),na.rm=TRUE))
	output <- data.frame(Compliance=paste("P>",as.character(cutTheta)),ES=round(caceES2,2),LB=round(bootES[1,],2),UB=round(bootES[2,],2))
	output2 <- list(CACE=output,Compliers=p1Table )
	output2$Method <- "MLM"
	class(output2) <- "eefAnalytics"
	return(output2 ) 
}

############


## - internal
crt.cace <- function(posttest,fixedDesignMatrix,intervention,cluster){
	
	freqFit <- lmer(posttest~ fixedDesignMatrix-1+ (1|cluster))
		
	betaB <- data.frame(cbind(summary(freqFit)$coefficients[,1]))
	row.names(betaB)<- colnames(fixedDesignMatrix)
	var.B<- as.numeric(summary(freqFit)$varcor)
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.B,var.W)
	sigmaBE <- sigmaBE
	names(sigmaBE)<- c("Schools","Pupils")
      #schRand <- data.frame(ranef(freqFit)$cluster)
	#names(schRand)<- "Estimate"
	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- g.within(var.w=var.W, beta=beta, icc=ICC , group=group, schoolID=cluster)
  		esTotal <- g.total(var.tt=var.tt, beta=beta, icc=ICC , group=group, schoolID=cluster)
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		colnames(output1) <- c("Estimate","95% LB","95% UB")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- round(output1,2)
	}
	names(output2) <- row.names(betaB)[btp]
      output <- list(ES=output2)

	return(output)
}

#####################
#################

#############################################################################
############# MST main functions ################################################
#' CACE Analysis of Multisite Randomised Education Trials.
#' 
#' \code{caceMSTBoot} performs exploratory CACE analysis of multisite randomised education trials.
#' 
#' @export
#' @param formula the model to be analysed. It is of the form y ~ x1+x2+.... Where y is the outcome variable and Xs are the predictors.
#' @param random a string variable specifying the "clustering variable" as contained in the data. See example below 
#' @param intervention a string variable specifying the "intervention variable" as appeared in the formula. See example below
#' @param compliance a string variable specifying the "compliance variable" as contained in the data. The data must be in percentages ranging from 0 - 100.
#' @param nBoot number of bootstraps required to generate bootstrap confidence interval. Default is NULL.
#' @param data data frame containing the data to be analysed.
#' @return S3 object; a list consisting of
#' \itemize{
#' \item \code{CACE}. Estimates of CACE adjusted effect sizes based on pre-specified thresholds. Only produced for threshold with at least 50% compliance rate.
#' \item \code{Compliers}. Percentage of pupils that achieved a pre-specified threshold of compliance.
#' }
#' @example inst/examples/mstExample.R
caceMSTBoot<- function(formula,random,intervention,compliance,nBoot,data)UseMethod("caceMSTBoot")


#' @export
caceMSTBoot.default <- function(formula,random,intervention,compliance,nBoot,data){stop("No correct formula input given.")}
  

#' @export
caceMSTBoot.formula <- function(formula,random,intervention,compliance,nBoot,data){

	data <- data[order(data[,which(colnames(data)==random)],data[,which(colnames(data)==intervention)]),]
	intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]
	tmp2 <- which(colnames(data)==random)
	cluster2 = data[,tmp2]

	chk <- sum(rowSums(table(cluster2,trt)!=0)>1)
	if(chk ==0){stop("This is not a MST design")}
 	
	tmp3 <- which(colnames(data)==intervention)
	data[,tmp3] <- as.factor(data[,tmp3])
	mf <- model.frame(formula=formula, data=data)
	mf <- mf[order(cluster2),]
	cluster <- cluster2[order(cluster2)]
	trt <- trt[order(cluster2)]
	fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
	tmp <- colnames(fixedDesignMatrix )
	tmp[1]  <- "Intercept"
	colnames(fixedDesignMatrix)<- tmp
	posttest <- model.response(mf)
      intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]
	if(length(table(trt))!=2){stop("Not applicable to multi-arm trials. Re-define the treatment group as dummy variables and run the function for each dummy variable separately")}
	comp1 <- data[,which(colnames(data)==compliance)]
	comp2 <- comp1[order(cluster2)]

	if(length(tmp2)!= 1){stop("Cluster variable misspecified")}
	if(length(tmp3)!= 1){stop("Intervention variable misspecified")}


	output <- rbd.cace(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt,cluster=cluster)	


	theta <- seq(0,ifelse(max(comp2)==100,90,max(comp2)),10)
	alpha_ITT_E <- output$ES[[1]][2,1]

	d <- sapply(theta,function(x)ifelse(comp2 > x,1,0))
	ibit <- data[,tmp3]
      tc.pp <- apply(d,2,function(x) table(x,ibit))
	tc.pp2 <- apply(d,2,function(x) nrow(table(x,ibit)))
	if(sum(tc.pp2==2) != length(tc.pp2)){
		t.tpp1 <- c.tpp1 <- p1 <- NULL
		for (j in 1:length(tc.pp2)){
			if(tc.pp2[j]==1){
				t.tpp1[j] <- 1
				c.tpp1[j] <- 0
				p1[j] <- 1
			}
			if(tc.pp2[j]==2){	
				tppk <- unlist(tc.pp[j])
      			c.tpp1[j] <- tppk[2]/(tppk[1]+tppk[2])
				t.tpp1[j] <-tppk[4]/(tppk[3]+tppk[4])
				p1[j] <- t.tpp1[j] - c.tpp1[j] 
			}
		}
	}

	if(sum(tc.pp2==2) == length(tc.pp2)){

      	c.tpp1 <- tc.pp[2,]/(tc.pp[1,]+tc.pp[2,])
		t.tpp1 <-tc.pp[4,]/(tc.pp[3,]+tc.pp[4,])
		p1 <- t.tpp1 - c.tpp1

	}

	caceES  <- alpha_ITT_E/(p1)
	p1Table <- rbind(t.tpp1,c.tpp1,p1)
	p1Table <- data.frame(round(p1Table,2))
	row.names(p1Table )<- c("pT","pC","P=PT-pC")
	colnames(p1Table ) <- paste("P > ", theta,sep="")

	p2Table <- p1Table[,round(p1,2)>0.5]
	cutTheta <- theta[round(p1,2)>0.5]
	caceES2 <- caceES[round(p1,2) > 0.5] 

	tid <- c(1:nrow(fixedDesignMatrix))
	set.seed(1020252)

	bootSamples <- NULL

	for(ii in 1:length(unique(cluster))){
		selID <- tid[cluster==unique(cluster)[ii]]
		if(length(selID)>0){
			selID2<- sapply(c(1:nBoot),function(x)sample(selID,length(selID),replace=TRUE))
			bootSamples <- rbind(bootSamples ,selID2)
		}
				
	}
	row.names(bootSamples ) <- NULL
	bootSamples <- data.frame(bootSamples )

	bESOutput <- NULL

	for(i in 1:ncol(bootSamples )){

		bData <- data[bootSamples[,i],]

		bcluster <- bData[,tmp2]
		bmf <- model.frame(formula=formula, data=bData)
		bfixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(bmf, "terms"), data=bData)))
		bposttest <- model.response(bmf)
		btrt <- bData[,which(colnames(bData)==intervention)]
		boutput <-rbd.cace(posttest=bposttest,fixedDesignMatrix=bfixedDesignMatrix,intervention=intervention,trt=btrt,cluster=bcluster)
	
		balpha_ITT_E <- boutput$ES[[1]][2,1]

		bcomp <- bData[,which(colnames(bData)==compliance)]
		bd <- sapply(cutTheta,function(x)ifelse(bcomp > x,1,0))
		bibit <- bData[,tmp3]
      	btc.pp <- apply(bd,2,function(x) table(x,bibit))
		btc.pp2 <- apply(bd,2,function(x) nrow(table(x,bibit)))
		if(sum(btc.pp2==2) != length(btc.pp2)){
			t.tpp1 <- c.tpp1 <- bp1 <- NULL
			for (j in 1:length(btc.pp2)){
				if(btc.pp2[j]==1){
					t.tpp1[j] <- 1
					c.tpp1[j] <- 0
					bp1[j] <- 1
				}
				if(btc.pp2[j]==2){	
					tppk <- unlist(btc.pp[j])
      				c.tpp1[j] <- tppk[2]/(tppk[1]+tppk[2])
					t.tpp1[j] <-tppk[4]/(tppk[3]+tppk[4])
					bp1[j] <- t.tpp1[j] - c.tpp1[j] 
				}
			}
		}

		if(sum(btc.pp2==2) == length(btc.pp2)){

      		c.tpp1 <- btc.pp[2,]/(btc.pp[1,]+btc.pp[2,])
			t.tpp1 <-btc.pp[4,]/(btc.pp[3,]+btc.pp[4,])
			bp1 <- t.tpp1 - c.tpp1

		}
		bcaceES  <- balpha_ITT_E/(bp1)
		bESOutput <- rbind(bESOutput ,bcaceES)  


	}

	bootES <- apply(bESOutput ,2,function(x)quantile(x,prob=c(0.025,0.975),na.rm=TRUE))
	output <- data.frame(Compliance=paste("P>",as.character(cutTheta)),ES=round(caceES2,2),LB=round(bootES[1,],2),UB=round(bootES[2,],2))
	output2 <- list(CACE=output,Compliers=p1Table )
	output2$Method <- "MLM"
	
	class(output2) <- "eefAnalytics"
	return(output2 ) 

}


##########


##  - internal

rbd.cace <- function(posttest,fixedDesignMatrix,intervention,trt,cluster){

	freqFit <- lmer(posttest~ fixedDesignMatrix-1+(1|trt:cluster)+(1|cluster))

	betaB <- data.frame(cbind(summary(freqFit)$coefficients[,1]))
	row.names(betaB)<- colnames(fixedDesignMatrix)
	betaB <- betaB
	var.B2<- as.matrix(summary(freqFit)$varcor)
	var.B3 <- c(c(matrix(attr(var.B2[[1]],"stddev"))),c(matrix(attr(var.B2[[2]],"stddev"))))
	var.B3 <- var.B3^2
	var.sch <- var.B3[1]
	var.schTrt <- var.B3[2]	
	var.B <- sum(var.B3)
	var.W<- summary(freqFit)$sigma^2
	var.tt <- var.W+var.B
	ICC <- var.B/var.tt
	sigmaBE <- c(var.sch,var.schTrt,var.W)
	sigmaBE <- sigmaBE
	names(sigmaBE)<- c("Schools","Intervention:School","Pupils")
      #schRand <- data.frame(ranef(freqFit)$cluster)
	#names(schRand)<- "Estimate"
	btp <- which(substring(row.names(betaB),1,nchar(intervention))==intervention)
	      
      
      output2 <- list()
	for( i in 1:length(btp )){
	
		beta <- betaB[btp[i],1]
      	group <- fixedDesignMatrix[,btp[i]]

  		esWithin <- g.within(var.w=var.W, beta=beta, icc=ICC , group=group, schoolID=cluster)
  		esTotal <- g.total(var.tt=var.tt, beta=beta, icc=ICC , group=group, schoolID=cluster)
   	
		output1 <- data.frame(rbind(esWithin,esTotal))
		colnames(output1) <- c("Estimate","95% LB","95% UB")
		rownames(output1) <- c("Within","Total")
		output2[[i]] <- round(output1,2)
	}
	names(output2) <- row.names(betaB)[btp]
      output <- list(ES=output2)


	return(output)
}



####

#' CACE Analysis of Simple Randomised Education Trials.
#' 
#' \code{caceSRTBoot} performs exploraty CACE analysis of simple randomised education trials.
#' 
#' @export
#' @param formula the model to be analysed. It is of the form y ~ x1+x2+.... Where y is the outcome variable and Xs are the predictors.
#' @param intervention a string variable specifying the "intervention variable" as appeared in the formula. See example below
#' @param compliance a string variable specifying the "compliance variable" as contained in the data.  The data must be in percentages ranging from 0 - 100.
#' @param nBoot number of non-parametric bootstraps.  Default is NULL. 
#' @param data data frame containing the data to be analysed.
#' @return S3 \code{mcpi} object; a list consisting of
#' \itemize{
#' \item \code{CACE}. Estimates of CACE adjusted effect sizes based on pre-specified thresholds. Only produced for threshold with atleast 50% compliance rate.
#' \item \code{Compliers} Percentage of pupils that achieved a pre-specified threshold of compliance.
#' }
#' @example inst/examples/srtCACEExample.R
caceSRTBoot <- function(formula,intervention,compliance,nBoot,data){UseMethod("caceSRTBoot")}
  

#' @export
caceSRTBoot.default <- function(formula,intervention,compliance,nBoot,data){stop("No correct formula input given.")}
  

#' @export
caceSRTBoot.formula <- function(formula,intervention,compliance,nBoot,data){

	tmp3 <- which(colnames(data)==intervention)
	data[,tmp3] <- as.factor(data[,tmp3])
	mf <- model.frame(formula=formula, data=data)
	fixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(mf, "terms"), data=data)))
	tmp <- colnames(fixedDesignMatrix )
	tmp[1]  <- "Intercept"
	colnames(fixedDesignMatrix)<- tmp
	posttest <- model.response(mf)
      intervention <- intervention
	trt <- data[,which(colnames(data)==intervention)]

	comp2 <- data[,which(colnames(data)==compliance)]


	output <- srt(posttest=posttest,fixedDesignMatrix=fixedDesignMatrix,intervention=intervention,trt=trt)


	theta <- seq(0,ifelse(max(comp2)==100,90,max(comp2)),10)
	alpha_ITT_E <- output$ES[1,1]

	d <- sapply(theta,function(x)ifelse(comp2 > x,1,0))
	ibit <- data[,tmp3]
      tc.pp <- apply(d,2,function(x) table(x,ibit))
	tc.pp2 <- apply(d,2,function(x) nrow(table(x,ibit)))
	if(sum(tc.pp2==2) != length(tc.pp2)){
		t.tpp1 <- c.tpp1 <- p1 <- NULL
		for (j in 1:length(tc.pp2)){
			if(tc.pp2[j]==1){
				t.tpp1[j] <- 1
				c.tpp1[j] <- 0
				p1[j] <- 1
			}
			if(tc.pp2[j]==2){	
				tppk <- unlist(tc.pp[j])
      			c.tpp1[j] <- tppk[2]/(tppk[1]+tppk[2])
				t.tpp1[j] <-tppk[4]/(tppk[3]+tppk[4])
				p1[j] <- t.tpp1[j] - c.tpp1[j] 
			}
		}
	}

	if(sum(tc.pp2==2) == length(tc.pp2)){

      	c.tpp1 <- tc.pp[2,]/(tc.pp[1,]+tc.pp[2,])
		t.tpp1 <-tc.pp[4,]/(tc.pp[3,]+tc.pp[4,])
		p1 <- t.tpp1 - c.tpp1

	}

	caceES  <- alpha_ITT_E/(p1)
	p1Table <- rbind(t.tpp1,c.tpp1,p1)
	p1Table <- data.frame(round(p1Table,2))
	row.names(p1Table )<- c("pT","pC","P=PT-pC")
	colnames(p1Table ) <- paste("P > ", theta,sep="")

	p2Table <- p1Table[,round(p1,2)>0.5]
	cutTheta <- theta[round(p1,2)>0.5]
	caceES2 <- caceES[round(p1,2) > 0.5] 

	tid <- c(1:nrow(fixedDesignMatrix))
	set.seed(1020252)
	tid <- c(1:nrow(fixedDesignMatrix))
	bootSamples <- sapply(c(1:nBoot),function(x)sample(tid,replace=TRUE))

	row.names(bootSamples ) <- NULL
	bootSamples <- data.frame(bootSamples )

	bESOutput <- NULL

	for(i in 1:ncol(bootSamples )){

		bData <- data[bootSamples[,i],]

		bmf <- model.frame(formula=formula, data=bData)
		bfixedDesignMatrix <- as.matrix(data.frame(model.matrix(attr(bmf, "terms"), data=bData)))
		bposttest <- model.response(bmf)
		btrt <- bData[,which(colnames(bData)==intervention)]
		boutput <-srt(posttest=bposttest,fixedDesignMatrix=bfixedDesignMatrix,intervention=intervention,trt=btrt)
	
		balpha_ITT_E <- boutput$ES[1,1]

		bcomp <- bData[,which(colnames(bData)==compliance)]
		bd <- sapply(cutTheta,function(x)ifelse(bcomp > x,1,0))
		bibit <- bData[,tmp3]
      	btc.pp <- apply(bd,2,function(x) table(x,bibit))
		btc.pp2 <- apply(bd,2,function(x) nrow(table(x,bibit)))
		if(sum(btc.pp2==2) != length(btc.pp2)){
			t.tpp1 <- c.tpp1 <- bp1 <- NULL
			for (j in 1:length(btc.pp2)){
				if(btc.pp2[j]==1){
					t.tpp1[j] <- 1
					c.tpp1[j] <- 0
					bp1[j] <- 1
				}
				if(btc.pp2[j]==2){	
					tppk <- unlist(btc.pp[j])
      				c.tpp1[j] <- tppk[2]/(tppk[1]+tppk[2])
					t.tpp1[j] <-tppk[4]/(tppk[3]+tppk[4])
					bp1[j] <- t.tpp1[j] - c.tpp1[j] 
				}
			}
		}

		if(sum(btc.pp2==2) == length(btc.pp2)){

      		c.tpp1 <- btc.pp[2,]/(btc.pp[1,]+btc.pp[2,])
			t.tpp1 <-btc.pp[4,]/(btc.pp[3,]+btc.pp[4,])
			bp1 <- t.tpp1 - c.tpp1

		}
		bcaceES  <- balpha_ITT_E/(bp1)
		bESOutput <- rbind(bESOutput ,bcaceES)  
		
	}

	bootES <- apply(bESOutput ,2,function(x)quantile(x,prob=c(0.025,0.975),na.rm=TRUE))
	output <- data.frame(Compliance=paste("P>",as.character(cutTheta)),ES=round(caceES2,2),LB=round(bootES[1,],2),UB=round(bootES[2,],2))
	output2 <- list(CACE=output,Compliers=p1Table )
	output2$Method <- "LM"
	class(output2) <- "eefAnalytics"
	return(output2 ) 

}



#####

#############################################################################
############# SRT main functions ################################################

Try the eefAnalytics package in your browser

Any scripts or data that you put into this service are public.

eefAnalytics documentation built on May 31, 2017, 4:17 a.m.