Nothing
#' 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 ################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.