Appendix B. R Simulation Script
library(lavaan)
library(semTools)
library(spatstat.utils)
library(tcltk)
library(Matrix)
library(gdata)
library(bain)
library(gtools)
library(MIIVsem)
library(simsem)
library(MplusAutomation)
library(metafor)
library(bayesDP)
library(RBesT)
library(MLmetrics)
# Replace function is needed to fill in parameters in Mplus script with specified values
loopReplace <- function(text, replacements) {
for (v in names(replacements)){
text <- gsub(sprintf(“\\[\\[%s\\]\\]”, v), replacements[[v]], text)
}
return(text)
}
set.seed(31092)
simulation1.coef<-NULL
simulation1.se<-NULL
simulation1.cover<-NULL
simulation1.bias<-NULL
n<-200
coef1<-1
coef2<-0.5
#SIMULATION MODEL#
population.model <- ‘ f1 =~ x1 + 0.6*x2 + 0.6*x3 + 0.6*x4 + 0.6*x5
f2 =~ x6 + 0.6*x7 + 0.6*x8 + 0.6*x9 + 0.6*x10
f3 =~ x11 + 0.6*x12 + 0.6*x13 + 0.6*x14 + 0.6*x15
f3 ~ 1*f1 + 0.5*f2′
population.model.hetero <- ‘ f1 =~ x1 + 0.6*x2 + 0.6*x3 + 0.6*x4 + 0.6*x5
f2 =~ x6 + 0.6*x7 + 0.6*x8 + 0.6*x9 + 0.6*x10
f3 =~ x11 + 0.6*x12 + 0.6*x13 + 0.6*x14 + 0.6*x15
f3 ~ 0.5*f1 + 0.25*f2
‘
for(z in 1:100) {
#################
##GENERATE DATA##
#################
sim.data1 <- simulateData(population.model.hetero, sample.nobs=n)
sim.data2 <- simulateData(population.model.hetero, sample.nobs=n)
sim.data3 <- simulateData(population.model.hetero, sample.nobs=n)
sim.data4 <- simulateData(population.model.hetero, sample.nobs=n)
sim.data5 <- simulateData(population.model.hetero, sample.nobs=n)
sim.data6 <- simulateData(population.model, sample.nobs=n)
####### Create MLE Syntax for First Data Set #######
mle.script <- mplusObject(
TITLE = “SEM for data1;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = ml;”,
MODEL = “
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 f2;
output: cinterval
“,
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data1)
##RUN MLE IN MPLUS##
mle.result = mplusModeler(mle.script, “sim.data1”, modelout = “mle.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR MLE##
mle.coef.f1<-mle.result$results$parameters$unstandardized[mle.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & mle.result$results$parameters$unstandardized$param == ‘F1′,3]
mle.se.f1<-mle.result$results$parameters$unstandardized[mle.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & mle.result$results$parameters$unstandardized$param == ‘F1′,4]
mle.coef.f2<-mle.result$results$parameters$unstandardized[mle.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & mle.result$results$parameters$unstandardized$param == ‘F2′,3]
mle.se.f2<-mle.result$results$parameters$unstandardized[mle.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & mle.result$results$parameters$unstandardized$param == ‘F2′,4]
mle.coef.f1
mle.se.f1
mle.coef.f2
mle.se.f2
mle.low.f1<-mle.result$results$parameters$ci.unstandardized[mle.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & mle.result$results$parameters$ci.unstandardized$param==‘F1′,4]
mle.high.f1<-mle.result$results$parameters$ci.unstandardized[mle.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & mle.result$results$parameters$ci.unstandardized$param==‘F1′,8]
mle.cover.f1<-ifelse(coef1>=mle.low.f1 & coef1<=mle.high.f1,1,0)
mle.low.f2<-mle.result$results$parameters$ci.unstandardized[mle.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & mle.result$results$parameters$ci.unstandardized$param==‘F2′,4]
mle.high.f2<-mle.result$results$parameters$ci.unstandardized[mle.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & mle.result$results$parameters$ci.unstandardized$param==‘F2′,8]
mle.cover.f2<-ifelse(coef2>=mle.low.f2 & coef2<=mle.high.f2,1,0)
####### Create Bayes Syntax for First Data Set #######
bayes1.script <- mplusObject(
TITLE = “SEM for data1;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = “
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 f2;
“,
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data1)
##RUN BAYES1 IN MPLUS##
bayes1.result = mplusModeler(bayes1.script, “sim.data1”, modelout = “bayes1.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
bayes1.coef.f1<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param == ‘F1′,3]
bayes1.se.f1<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param == ‘F1′,4]
bayes1.coef.f2<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param == ‘F2′,3]
bayes1.se.f2<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param == ‘F2′,4]
bayes1.coef.f1
bayes1.se.f1
bayes1.coef.f2
bayes1.se.f2
bayes1.var.f1<-bayes1.se.f1^2
bayes1.var.f2<-bayes1.se.f2^2
df.f<-data.frame(bayes1.coef.f1,bayes1.coef.f2,bayes1.var.f1,bayes1.var.f2)
bayes1.low.f1<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param==‘F1′,6]
bayes1.high.f1<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param==‘F1′,7]
bayes1.cover.f1<-ifelse(coef1>=bayes1.low.f1 & coef1<=bayes1.high.f1,1,0)
bayes1.low.f2<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param==‘F2′,6]
bayes1.high.f2<-bayes1.result$results$parameters$unstandardized[bayes1.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes1.result$results$parameters$unstandardized$param==‘F2′,7]
bayes1.cover.f2<-ifelse(coef2>=bayes1.low.f2 & coef2<=bayes1.high.f2,1,0)
####### Create Bayes Syntax for Second Data Set #######
bayes2.script <- mplusObject(
TITLE = “SEM for data2;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = loopReplace(“
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 (beta1);
f3 on f2 (beta2);
MODEL PRIORS:
beta1~N([[bayes1.coef.f1]], [[bayes1.var.f1]]);
beta2~N([[bayes1.coef.f2]], [[bayes1.var.f2]]);
“, df.f),
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data2)
##RUN BAYES2 IN MPLUS##
bayes2.result = mplusModeler(bayes2.script, “sim.data2”, modelout = “bayes2.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
bayes2.coef.f1<-bayes2.result$results$parameters$unstandardized[bayes2.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes2.result$results$parameters$unstandardized$param == ‘F1′,3]
bayes2.se.f1<-bayes2.result$results$parameters$unstandardized[bayes2.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes2.result$results$parameters$unstandardized$param == ‘F1′,4]
bayes2.coef.f2<-bayes2.result$results$parameters$unstandardized[bayes2.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes2.result$results$parameters$unstandardized$param == ‘F2′,3]
bayes2.se.f2<-bayes2.result$results$parameters$unstandardized[bayes2.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes2.result$results$parameters$unstandardized$param == ‘F2′,4]
bayes2.coef.f1
bayes2.se.f1
bayes2.coef.f2
bayes2.se.f2
bayes2.var.f1<-bayes2.se.f1^2
bayes2.var.f2<-bayes2.se.f2^2
df.f<-data.frame(bayes2.coef.f1,bayes2.coef.f2,bayes2.var.f1,bayes2.var.f2)
####### Create Bayes Syntax for Third Data Set #######
bayes3.script <- mplusObject(
TITLE = “SEM for data3;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = loopReplace(“
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 (beta1);
f3 on f2 (beta2);
MODEL PRIORS:
beta1~N([[bayes2.coef.f1]], [[bayes2.var.f1]]);
beta2~N([[bayes2.coef.f2]], [[bayes2.var.f2]]);
“, df.f),
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data3)
##RUN BAYES2 IN MPLUS##
bayes3.result = mplusModeler(bayes3.script, “sim.data3”, modelout = “bayes3.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
bayes3.coef.f1<-bayes3.result$results$parameters$unstandardized[bayes3.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes3.result$results$parameters$unstandardized$param == ‘F1′,3]
bayes3.se.f1<-bayes3.result$results$parameters$unstandardized[bayes3.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes3.result$results$parameters$unstandardized$param == ‘F1′,4]
bayes3.coef.f2<-bayes3.result$results$parameters$unstandardized[bayes3.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes3.result$results$parameters$unstandardized$param == ‘F2′,3]
bayes3.se.f2<-bayes3.result$results$parameters$unstandardized[bayes3.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes3.result$results$parameters$unstandardized$param == ‘F2′,4]
bayes3.coef.f1
bayes3.se.f1
bayes3.coef.f2
bayes3.se.f2
bayes3.var.f1<-bayes3.se.f1^2
bayes3.var.f2<-bayes3.se.f2^2
df.f<-data.frame(bayes3.coef.f1,bayes3.coef.f2,bayes3.var.f1,bayes3.var.f2)
####### Create Bayes Syntax for Fourth Data Set #######
bayes4.script <- mplusObject(
TITLE = “SEM for data4;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = loopReplace(“
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 (beta1);
f3 on f2 (beta2);
MODEL PRIORS:
beta1~N([[bayes3.coef.f1]], [[bayes3.var.f1]]);
beta2~N([[bayes3.coef.f2]], [[bayes3.var.f2]]);
“, df.f),
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data4)
##RUN BAYES2 IN MPLUS##
bayes4.result = mplusModeler(bayes4.script, “sim.data4”, modelout = “bayes4.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
bayes4.coef.f1<-bayes4.result$results$parameters$unstandardized[bayes4.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes4.result$results$parameters$unstandardized$param == ‘F1′,3]
bayes4.se.f1<-bayes4.result$results$parameters$unstandardized[bayes4.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes4.result$results$parameters$unstandardized$param == ‘F1′,4]
bayes4.coef.f2<-bayes4.result$results$parameters$unstandardized[bayes4.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes4.result$results$parameters$unstandardized$param == ‘F2′,3]
bayes4.se.f2<-bayes4.result$results$parameters$unstandardized[bayes4.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes4.result$results$parameters$unstandardized$param == ‘F2′,4]
bayes4.coef.f1
bayes4.se.f1
bayes4.coef.f2
bayes4.se.f2
bayes4.var.f1<-bayes3.se.f1^2
bayes4.var.f2<-bayes3.se.f2^2
df.f<-data.frame(bayes4.coef.f1,bayes4.coef.f2,bayes4.var.f1,bayes4.var.f2)
####### Create Bayes Syntax for Fifth Data Set #######
bayes5.script <- mplusObject(
TITLE = “SEM for bayes5;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = loopReplace(“
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 (beta1);
f3 on f2 (beta2);
MODEL PRIORS:
beta1~N([[bayes4.coef.f1]], [[bayes4.var.f1]]);
beta2~N([[bayes4.coef.f2]], [[bayes4.var.f2]]);
“, df.f),
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data5)
##RUN BAYES2 IN MPLUS##
bayes5.result = mplusModeler(bayes5.script, “sim.data5”, modelout = “bayes5.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
bayes5.coef.f1<-bayes5.result$results$parameters$unstandardized[bayes5.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes5.result$results$parameters$unstandardized$param == ‘F1′,3]
bayes5.se.f1<-bayes5.result$results$parameters$unstandardized[bayes5.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes5.result$results$parameters$unstandardized$param == ‘F1′,4]
bayes5.coef.f2<-bayes5.result$results$parameters$unstandardized[bayes5.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes5.result$results$parameters$unstandardized$param == ‘F2′,3]
bayes5.se.f2<-bayes5.result$results$parameters$unstandardized[bayes5.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes5.result$results$parameters$unstandardized$param == ‘F2′,4]
bayes5.coef.f1
bayes5.se.f1
bayes5.coef.f2
bayes5.se.f2
bayes5.var.f1<-bayes5.se.f1^2
bayes5.var.f2<-bayes5.se.f2^2
df.f<-data.frame(bayes5.coef.f1,bayes5.coef.f2,bayes5.var.f1,bayes5.var.f2)
####### Create Bayes Syntax for Sixth Data Set #######
bayes6.script <- mplusObject(
TITLE = “SEM for bayes6;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = loopReplace(“
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 (beta1);
f3 on f2 (beta2);
MODEL PRIORS:
beta1~N([[bayes5.coef.f1]], [[bayes5.var.f1]]);
beta2~N([[bayes5.coef.f2]], [[bayes5.var.f2]]);
“, df.f),
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data6)
##RUN BAYES2 IN MPLUS##
bayes6.result = mplusModeler(bayes6.script, “sim.data6”, modelout = “bayes6.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
bayes6.coef.f1<-bayes6.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param == ‘F1′,3]
bayes6.se.f1<-bayes6.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param == ‘F1′,4]
bayes6.coef.f2<-bayes6.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param == ‘F2′,3]
bayes6.se.f2<-bayes6.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param == ‘F2′,4]
bayes6.coef.f1
bayes6.se.f1
bayes6.coef.f2
bayes6.se.f2
bayes6.var.f1<-bayes6.se.f1^2
bayes6.var.f2<-bayes6.se.f2^2
df.f<-data.frame(bayes6.coef.f1,bayes6.coef.f2,bayes6.var.f1,bayes6.var.f2)
bayes6.low.f1<-bayes1.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param==‘F1′,6]
bayes6.high.f1<-bayes1.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param==‘F1′,7]
bayes6.cover.f1<-ifelse(coef1>=bayes6.low.f1 & coef1<=bayes6.high.f1,1,0)
bayes6.low.f2<-bayes6.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param==‘F2′,6]
bayes6.high.f2<-bayes6.result$results$parameters$unstandardized[bayes6.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & bayes6.result$results$parameters$unstandardized$param==‘F2′,7]
bayes6.cover.f2<-ifelse(coef2>=bayes6.low.f2 & coef2<=bayes6.high.f2,1,0)
###DATA POOLING###
pooled.data<-rbind(sim.data1,sim.data2,sim.data3,sim.data4,sim.data5,sim.data6)
####### Create ML Syntax for Pooled Data Set #######
ml.pooled.script <- mplusObject(
TITLE = “SEM for pooled data;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = ml;”,
MODEL = “
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 f2;
output: cinterval”,
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=pooled.data)
##RUN BAYES2 IN MPLUS##
ml.pooled.result = mplusModeler(ml.pooled.script, “pooled.data”, modelout = “mlpooled.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
ml.pooled.coef.f1<-ml.pooled.result$results$parameters$unstandardized[ml.pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & ml.pooled.result$results$parameters$unstandardized$param == ‘F1′,3]
ml.pooled.se.f1<-ml.pooled.result$results$parameters$unstandardized[ml.pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & ml.pooled.result$results$parameters$unstandardized$param == ‘F1′,4]
ml.pooled.coef.f2<-ml.pooled.result$results$parameters$unstandardized[ml.pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & ml.pooled.result$results$parameters$unstandardized$param == ‘F2′,3]
ml.pooled.se.f2<-ml.pooled.result$results$parameters$unstandardized[ml.pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & ml.pooled.result$results$parameters$unstandardized$param == ‘F2′,4]
ml.pooled.coef.f1
ml.pooled.se.f1
ml.pooled.coef.f2
ml.pooled.se.f2
ml.pooled.low.f1<-ml.pooled.result$results$parameters$ci.unstandardized[ml.pooled.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & ml.pooled.result$results$parameters$ci.unstandardized$param==‘F1′,4]
ml.pooled.high.f1<-ml.pooled.result$results$parameters$ci.unstandardized[ml.pooled.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & ml.pooled.result$results$parameters$ci.unstandardized$param==‘F1′,8]
ml.pooled.cover.f1<-ifelse(coef1>=ml.pooled.low.f1 & coef1<=ml.pooled.high.f1,1,0)
ml.pooled.low.f2<-ml.pooled.result$results$parameters$ci.unstandardized[ml.pooled.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & ml.pooled.result$results$parameters$ci.unstandardized$param==‘F2′,4]
ml.pooled.high.f2<-ml.pooled.result$results$parameters$ci.unstandardized[ml.pooled.result$results$parameters$ci.unstandardized$paramHeader==‘F3.ON’ & ml.pooled.result$results$parameters$ci.unstandardized$param==‘F2′,8]
ml.pooled.cover.f2<-ifelse(coef2>=ml.pooled.low.f2 & coef2<=ml.pooled.high.f2,1,0)
####### Create Bayes Syntax for Pooled Data Set #######
pooled.script <- mplusObject(
TITLE = “SEM for pooled data;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = “
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 f2;
“,
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=pooled.data)
##RUN BAYES2 IN MPLUS##
pooled.result = mplusModeler(pooled.script, “pooled.data”, modelout = “pooled.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
pooled.coef.f1<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & pooled.result$results$parameters$unstandardized$param == ‘F1′,3]
pooled.se.f1<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & pooled.result$results$parameters$unstandardized$param == ‘F1′,4]
pooled.coef.f2<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & pooled.result$results$parameters$unstandardized$param == ‘F2′,3]
pooled.se.f2<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & pooled.result$results$parameters$unstandardized$param == ‘F2′,4]
pooled.coef.f1
pooled.se.f1
pooled.coef.f2
pooled.se.f2
pooled.low.f1<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & pooled.result$results$parameters$unstandardized$param==‘F1′,6]
pooled.high.f1<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & pooled.result$results$parameters$unstandardized$param==‘F1′,7]
pooled.cover.f1<-ifelse(coef1>=pooled.low.f1 & coef1<=pooled.high.f1,1,0)
pooled.low.f2<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & pooled.result$results$parameters$unstandardized$param==‘F2′,6]
pooled.high.f2<-pooled.result$results$parameters$unstandardized[pooled.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & pooled.result$results$parameters$unstandardized$param==‘F2′,7]
pooled.cover.f2<-ifelse(coef2>=pooled.low.f2 & coef2<=pooled.high.f2,1,0)
###AGDP USE OF AVERAGE PRIOR ESTIMATES###
bayes.coef.f1.mean<-mean(bayes1.coef.f1,bayes2.coef.f1,bayes3.coef.f1,bayes4.coef.f1,bayes5.coef.f1)
bayes.coef.f2.mean<-mean(bayes1.coef.f2,bayes2.coef.f2,bayes3.coef.f2,bayes4.coef.f2,bayes5.coef.f2)
bayes.var.f1.mean<-mean(bayes1.var.f1,bayes2.var.f1,bayes3.var.f1,bayes4.var.f1,bayes5.var.f1)
bayes.var.f2.mean<-mean(bayes1.var.f2,bayes2.var.f2,bayes3.var.f2,bayes4.var.f2,bayes5.var.f2)
df.f<-data.frame(bayes.coef.f1.mean,bayes.coef.f2.mean,bayes.var.f1.mean,bayes.var.f2.mean)
agdp.script <- mplusObject(
TITLE = “SEM for AGDP;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = loopReplace(“
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 (beta1);
f3 on f2 (beta2);
MODEL PRIORS:
beta1~N([[bayes.coef.f1.mean]], [[bayes.var.f1.mean]]);
beta2~N([[bayes.coef.f2.mean]], [[bayes.var.f2.mean]]);
“, df.f),
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data6)
##RUN AGDP BAYES IN MPLUS##
agdp.result = mplusModeler(agdp.script, “sim.data6”, modelout = “agdp.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES##
agdp.coef.f1<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & agdp.result$results$parameters$unstandardized$param == ‘F1′,3]
agdp.se.f1<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & agdp.result$results$parameters$unstandardized$param == ‘F1′,4]
agdp.coef.f2<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & agdp.result$results$parameters$unstandardized$param == ‘F2′,3]
agdp.se.f2<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & agdp.result$results$parameters$unstandardized$param == ‘F2′,4]
agdp.coef.f1
agdp.se.f1
agdp.coef.f2
agdp.se.f2
agdp.low.f1<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & agdp.result$results$parameters$unstandardized$param==‘F1′,6]
agdp.high.f1<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & agdp.result$results$parameters$unstandardized$param==‘F1′,7]
agdp.cover.f1<-ifelse(coef1>=agdp.low.f1 & coef1<=agdp.high.f1,1,0)
agdp.low.f2<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & agdp.result$results$parameters$unstandardized$param==‘F2′,6]
agdp.high.f2<-agdp.result$results$parameters$unstandardized[agdp.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & agdp.result$results$parameters$unstandardized$param==‘F2′,7]
agdp.cover.f2<-ifelse(coef2>=agdp.low.f2 & coef2<=agdp.high.f2,1,0)
###PARTIAL DISCOUNTING PRIORS USING bayesDP###
##BASE THE HISTORICAL DATA ON THE AGDP VALUES##
####### Create Bayes Syntax for Sixth Data Set using Naive priors #######
bayes6b.script <- mplusObject(
TITLE = “SEM for data6 naive priors;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = “
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 f2;
“,
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data6)
##RUN BAYES1 IN MPLUS##
bayes6b.result = mplusModeler(bayes6b.script, “sim.data6”, modelout = “bayes6b.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
bayes6b.coef.f1<-bayes6b.result$results$parameters$unstandardized[bayes6b.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6b.result$results$parameters$unstandardized$param == ‘F1′,3]
bayes6b.se.f1<-bayes6b.result$results$parameters$unstandardized[bayes6b.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6b.result$results$parameters$unstandardized$param == ‘F1′,4]
bayes6b.coef.f2<-bayes6b.result$results$parameters$unstandardized[bayes6b.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6b.result$results$parameters$unstandardized$param == ‘F2′,3]
bayes6b.se.f2<-bayes1.result$results$parameters$unstandardized[bayes6b.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & bayes6b.result$results$parameters$unstandardized$param == ‘F2′,4]
bayes_dp.random.coef.f1.fit<-bdpnormal( mu_t=bayes6b.coef.f1,sigma_t=bayes6b.se.f1,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f1,sigma0_t=sqrt(pooled.se.f1),N0_t=nrow(pooled.data), method=“mc” )
bayes_dp.random.coef.f1 <- round(median(bayes_dp.random.coef.f1.fit$posterior_treatment$posterior_mu),4)
bayes_dp.random.se.f1 <- round(mean(bayes_dp.random.coef.f1.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.random.coef.f1.fit)
bayes_dp.random.coef.f2.fit<-bdpnormal( mu_t=bayes6b.coef.f2,sigma_t=bayes6b.se.f2,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f2,sigma0_t=sqrt(pooled.se.f2),N0_t=nrow(pooled.data), method=“mc” )
bayes_dp.random.coef.f2 <- round(median(bayes_dp.random.coef.f2.fit$posterior_treatment$posterior_mu),4)
bayes_dp.random.se.f2 <- (mean(bayes_dp.random.coef.f2.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.random.coef.f2.fit)
bayes_dp.75.coef.f1.fit<-bdpnormal( mu_t=bayes6b.coef.f1,sigma_t=bayes6b.se.f1,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f1,sigma0_t=sqrt(pooled.se.f1),N0_t=nrow(pooled.data), method=“mc”, alpha_max=0.75, fix_alpha=TRUE)
bayes_dp.75.coef.f1 <- round(median(bayes_dp.75.coef.f1.fit$posterior_treatment$posterior_mu),4)
bayes_dp.75.se.f1 <- (mean(bayes_dp.75.coef.f1.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.75.coef.f1.fit)
bayes_dp.75.coef.f2.fit<-bdpnormal( mu_t=bayes6b.coef.f2,sigma_t=bayes6b.se.f2,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f2,sigma0_t=sqrt(pooled.se.f2),N0_t=nrow(pooled.data), method=“mc”, alpha_max=0.75, fix_alpha=TRUE )
bayes_dp.75.coef.f2 <- round(median(bayes_dp.75.coef.f2.fit$posterior_treatment$posterior_mu),4)
bayes_dp.75.se.f2 <- (mean(bayes_dp.75.coef.f2.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.75.coef.f2.fit)
bayes_dp.50.coef.f1.fit<-bdpnormal( mu_t=bayes6b.coef.f1,sigma_t=bayes6b.se.f1,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f1,sigma0_t=sqrt(pooled.se.f1),N0_t=nrow(pooled.data), method=“mc”, alpha_max=0.50, fix_alpha=TRUE)
bayes_dp.50.coef.f1 <- round(median(bayes_dp.50.coef.f1.fit$posterior_treatment$posterior_mu),4)
bayes_dp.50.se.f1 <- (mean(bayes_dp.50.coef.f1.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.50.coef.f1.fit)
bayes_dp.50.coef.f2.fit<-bdpnormal( mu_t=bayes6b.coef.f2,sigma_t=bayes6b.se.f2,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f2,sigma0_t=sqrt(pooled.se.f2),N0_t=nrow(pooled.data), method=“mc”, alpha_max=0.50, fix_alpha=TRUE )
bayes_dp.50.coef.f2 <- round(median(bayes_dp.50.coef.f2.fit$posterior_treatment$posterior_mu),4)
bayes_dp.50.se.f2 <- (mean(bayes_dp.50.coef.f2.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.50.coef.f2.fit)
bayes_dp.25.coef.f1.fit<-bdpnormal( mu_t=bayes6b.coef.f1,sigma_t=bayes6b.se.f1,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f1,sigma0_t=sqrt(pooled.se.f1),N0_t=nrow(pooled.data), method=“mc”, alpha_max=0.25, fix_alpha=TRUE)
bayes_dp.25.coef.f1 <- round(median(bayes_dp.25.coef.f1.fit$posterior_treatment$posterior_mu),4)
bayes_dp.25.se.f1 <- (mean(bayes_dp.25.coef.f1.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.25.coef.f1.fit)
bayes_dp.25.coef.f2.fit<-bdpnormal( mu_t=bayes6b.coef.f2,sigma_t=bayes6b.se.f2,N_t=nrow(sim.data6),
mu0_t=pooled.coef.f2,sigma0_t=sqrt(pooled.se.f2),N0_t=nrow(pooled.data), method=“mc”, alpha_max=0.25, fix_alpha=TRUE )
bayes_dp.25.coef.f2 <- round(median(bayes_dp.25.coef.f2.fit$posterior_treatment$posterior_mu),4)
bayes_dp.25.se.f2 <- (mean(bayes_dp.25.coef.f2.fit$posterior_treatment$posterior_sigma2))
summary(bayes_dp.25.coef.f2.fit)
bayes_dp.random.low.f1<-quantile(bayes_dp.random.coef.f1.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.random.high.f1<-quantile(bayes_dp.random.coef.f1.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.random.cover.f1<-ifelse(coef1>=bayes_dp.random.low.f1 & coef1<=bayes_dp.random.high.f1,1,0)
bayes_dp.random.low.f2<-quantile(bayes_dp.random.coef.f2.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.random.high.f2<-quantile(bayes_dp.random.coef.f2.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.random.cover.f2<-ifelse(coef2>=bayes_dp.random.low.f2 & coef2<=bayes_dp.random.high.f2,1,0)
bayes_dp.75.low.f1<-quantile(bayes_dp.75.coef.f1.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.75.high.f1<-quantile(bayes_dp.75.coef.f1.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.75.cover.f1<-ifelse(coef1>=bayes_dp.75.low.f1 & coef1<=bayes_dp.75.high.f1,1,0)
bayes_dp.75.low.f2<-quantile(bayes_dp.75.coef.f2.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.75.high.f2<-quantile(bayes_dp.75.coef.f2.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.75.cover.f2<-ifelse(coef2>=bayes_dp.75.low.f2 & coef2<=bayes_dp.75.high.f2,1,0)
bayes_dp.50.low.f1<-quantile(bayes_dp.50.coef.f1.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.50.high.f1<-quantile(bayes_dp.50.coef.f1.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.50.cover.f1<-ifelse(coef1>=bayes_dp.50.low.f1 & coef1<=bayes_dp.50.high.f1,1,0)
bayes_dp.50.low.f2<-quantile(bayes_dp.50.coef.f2.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.50.high.f2<-quantile(bayes_dp.50.coef.f2.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.50.cover.f2<-ifelse(coef2>=bayes_dp.50.low.f2 & coef2<=bayes_dp.50.high.f2,1,0)
bayes_dp.25.low.f1<-quantile(bayes_dp.25.coef.f1.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.25.high.f1<-quantile(bayes_dp.25.coef.f1.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.25.cover.f1<-ifelse(coef1>=bayes_dp.25.low.f1 & coef1<=bayes_dp.25.high.f1,1,0)
bayes_dp.25.low.f2<-quantile(bayes_dp.25.coef.f2.fit$posterior_treatment$posterior_mu, .025)
bayes_dp.25.high.f2<-quantile(bayes_dp.25.coef.f2.fit$posterior_treatment$posterior_mu, .975)
bayes_dp.25.cover.f2<-ifelse(coef2>=bayes_dp.25.low.f2 & coef2<=bayes_dp.25.high.f2,1,0)
###META ANALYSIS FOR BAYES PRIORS###
#CREATE FILE FOR PRIOR COEFFICIENT AND STANDARD ERROR ESTIMATES#
bayes.coef.f1.meta<-c(bayes1.coef.f1,bayes2.coef.f1,bayes3.coef.f1,bayes4.coef.f1,bayes5.coef.f1)
bayes.coef.f2.meta<-c(bayes1.coef.f2,bayes2.coef.f2,bayes3.coef.f2,bayes4.coef.f2,bayes5.coef.f2)
bayes.se.f1.meta<-c(bayes1.se.f1,bayes2.se.f1,bayes3.se.f1,bayes4.se.f1,bayes5.se.f1)
bayes.se.f2.meta<-c(bayes1.se.f2,bayes2.se.f2,bayes3.se.f2,bayes4.se.f2,bayes5.se.f2)
df.meta<-data.frame(bayes.coef.f1.meta,bayes.coef.f2.meta,bayes.se.f1.meta,bayes.se.f2.meta)
#RUN META ANALYSIS FOR PREVIOUS ESTIMATES#
f1.meta<-rma(yi=bayes.coef.f1.meta, sei=bayes.se.f1.meta, data=df.meta, method=“DL”)
coef.f1.meta<-as.numeric(f1.meta$beta)
se.f1.meta<-as.numeric(f1.meta$se)
f2.meta<-rma(yi=bayes.coef.f2.meta, sei=bayes.se.f2.meta, data=df.meta, method=“DL”)
coef.f2.meta<-as.numeric(f2.meta$beta)
se.f2.meta<-as.numeric(f2.meta$se)
meta.var.f1<-se.f1.meta^2
meta.var.f2<-se.f2.meta^2
df.f<-data.frame(coef.f1.meta,coef.f2.meta,meta.var.f1,meta.var.f2)
####### Create Bayes Syntax for Sixth Data Set for META ANALYSIS #######
meta.script <- mplusObject(
TITLE = “SEM for META;”,
VARIABLE = “USEVARIABLES = x1-x15;”,
ANALYSIS = “ESTIMATOR = bayes;”,
MODEL = loopReplace(“
f1 by x1-x5;
f2 by x6-x10;
f3 by x11-x15;
f3 on f1 (beta1);
f3 on f2 (beta2);
MODEL PRIORS:
beta1~N([[coef.f1.meta]], [[meta.var.f1]]);
beta2~N([[coef.f2.meta]], [[meta.var.f2]]);
“, df.f),
usevariables=c(“x1”,”x2”,”x3”,”x4”,”x5”,”x6”,”x7”,”x8”,”x9”,”x10”,
“x11”,”x12”,”x13”,”x14”,”x15”),
rdata=sim.data6)
##RUN META BAYES IN MPLUS##
meta.result = mplusModeler(meta.script, “sim.data6”, modelout = “meta.inp”,run = 1L)
##EXTRACT STRUCTURE COEFFICIENTS AND STANDARD ERRORS FOR BAYES1##
meta.coef.f1<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & meta.result$results$parameters$unstandardized$param == ‘F1′,3]
meta.se.f1<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & meta.result$results$parameters$unstandardized$param == ‘F1′,4]
meta.coef.f2<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & meta.result$results$parameters$unstandardized$param == ‘F2′,3]
meta.se.f2<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader == ‘F3.ON’ & meta.result$results$parameters$unstandardized$param == ‘F2′,4]
meta.coef.f1
meta.se.f1
meta.coef.f2
meta.se.f2
meta.low.f1<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & meta.result$results$parameters$unstandardized$param==‘F1′,6]
meta.high.f1<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & meta.result$results$parameters$unstandardized$param==‘F1′,7]
meta.cover.f1<-ifelse(coef1>=meta.low.f1 & coef1<=meta.high.f1,1,0)
meta.low.f2<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & meta.result$results$parameters$unstandardized$param==‘F2′,6]
meta.high.f2<-meta.result$results$parameters$unstandardized[meta.result$results$parameters$unstandardized$paramHeader==‘F3.ON’ & meta.result$results$parameters$unstandardized$param==‘F2′,7]
meta.cover.f2<-ifelse(coef2>=meta.low.f2 & coef2<=meta.high.f2,1,0)
###RBEST ANALYSIS###
rbest.study<-c(1,2,3,4,5)
rbest.n<-c(nrow(sim.data1),nrow(sim.data2),nrow(sim.data3),nrow(sim.data4),nrow(sim.data5))
rbest.coef.f1<-c(bayes1.coef.f1,bayes2.coef.f1,bayes3.coef.f1,bayes4.coef.f1,bayes5.coef.f1)
rbest.se.f1<-c(bayes1.se.f1,bayes2.se.f1,bayes3.se.f1,bayes4.se.f1,bayes5.se.f1)
rbest.coef.f2<-c(bayes1.coef.f2,bayes2.coef.f2,bayes3.coef.f2,bayes4.coef.f2,bayes5.coef.f2)
rbest.se.f2<-c(bayes1.se.f2,bayes2.se.f2,bayes3.se.f2,bayes4.se.f2,bayes5.se.f2)
df.rbest.f1<-data.frame(rbest.coef.f1,rbest.se.f1,rbest.n,rbest.study)
df.rbest.f2<-data.frame(rbest.coef.f2,rbest.se.f2,rbest.n,rbest.study)
map_f1<-gMAP(cbind(rbest.coef.f1,rbest.se.f1)~1|rbest.study, family=gaussian,
data=df.rbest.f1, weights=rbest.n,
tau.dist=“HalfNormal”,tau.prior=44, beta.prior=cbind(0,88))
map_f2<-gMAP(cbind(rbest.coef.f2,rbest.se.f2)~1|rbest.study, family=gaussian,
data=df.rbest.f2, weights=rbest.n,
tau.dist=“HalfNormal”,tau.prior=44, beta.prior=cbind(0,88))
map_f1.auto <- automixfit(map_f1, Nc=3)
map_f1.auto
map_f1.robust <- robustify(map_f1.auto, weight = .1, mean = 0, sigma=10)
map_f1.robust
map_f2.auto <- automixfit(map_f2, Nc=3)
map_f2.auto
map_f2.robust <- robustify(map_f2.auto, weight = .1, mean = 0, sigma=10)
map_f2.robust
rbest.auto.f1<-postmix(map_f1.auto,m=bayes6.coef.f1,n=nrow(sim.data6))
rbest.robust.f1<-postmix(map_f1.robust,m=bayes6.coef.f1,n=nrow(sim.data6))
rbest.auto.f2<-postmix(map_f2.auto,m=bayes6.coef.f2,n=nrow(sim.data6))
rbest.robust.f2<-postmix(map_f2.robust,m=bayes6.coef.f2,n=nrow(sim.data6))
rbest.auto.coef.f1<-rbest.auto.f1[1,1]*rbest.auto.f1[2,1]+rbest.auto.f1[1,2]*rbest.auto.f1[2,2]+rbest.auto.f1[1,3]*rbest.auto.f1[2,3]
rbest.auto.se.f1<-rbest.auto.f1[1,1]*rbest.auto.f1[3,1]+rbest.auto.f1[1,2]*rbest.auto.f1[3,2]+rbest.auto.f1[1,3]*rbest.auto.f1[3,3]
rbest.auto.coef.f2<-rbest.auto.f2[1,1]*rbest.auto.f2[2,1]+rbest.auto.f2[1,2]*rbest.auto.f2[2,2]+rbest.auto.f2[1,3]*rbest.auto.f2[2,3]
rbest.auto.se.f2<-rbest.auto.f2[1,1]*rbest.auto.f2[3,1]+rbest.auto.f2[1,2]*rbest.auto.f2[3,2]+rbest.auto.f2[1,3]*rbest.auto.f2[3,3]
rbest.robust.coef.f1<-rbest.robust.f1[1,1]*rbest.robust.f1[2,1]+rbest.robust.f1[1,2]*rbest.robust.f1[2,2]+rbest.robust.f1[1,3]*rbest.robust.f1[2,3]
rbest.robust.se.f1<-rbest.robust.f1[1,1]*rbest.robust.f1[3,1]+rbest.robust.f1[1,2]*rbest.robust.f1[3,2]+rbest.robust.f1[1,3]*rbest.robust.f1[3,3]
rbest.robust.coef.f2<-rbest.robust.f2[1,1]*rbest.robust.f2[2,1]+rbest.robust.f2[1,2]*rbest.robust.f2[2,2]+rbest.robust.f2[1,3]*rbest.robust.f2[2,3]
rbest.robust.se.f2<-rbest.robust.f2[1,1]*rbest.robust.f2[3,1]+rbest.robust.f2[1,2]*rbest.robust.f2[3,2]+rbest.robust.f2[1,3]*rbest.robust.f2[3,3]
rbest.auto.f1.hi<-rbest.auto.coef.f1+2*rbest.auto.se.f1
rbest.auto.f1.lo<-rbest.auto.coef.f1-2*rbest.auto.se.f1
rbest.robust.f1.hi<-rbest.robust.coef.f1+2*rbest.robust.se.f1
rbest.robust.f1.lo<-rbest.robust.coef.f1-2*rbest.robust.se.f1
rbest.auto.f2.hi<-rbest.auto.coef.f2+2*rbest.auto.se.f2
rbest.auto.f2.lo<-rbest.auto.coef.f2-2*rbest.auto.se.f2
rbest.robust.f2.hi<-rbest.robust.coef.f2+2*rbest.robust.se.f2
rbest.robust.f2.lo<-rbest.robust.coef.f2-2*rbest.robust.se.f2
rbest.auto.f1.cover<-ifelse((coef1>=rbest.auto.f1.lo & coef1<=rbest.auto.f1.hi),1,0)
rbest.auto.f2.cover<-ifelse((coef2>=rbest.auto.f2.lo & coef2<=rbest.auto.f2.hi),1,0)
rbest.robust.f1.cover<-ifelse((coef1>=rbest.robust.f1.lo & coef1<=rbest.robust.f1.hi),1,0)
rbest.robust.f2.cover<-ifelse((coef2>=rbest.robust.f2.lo & coef2<=rbest.robust.f2.hi),1,0)
###CALCULATE BIAS###
mle.bias.f1<-(mle.coef.f1-coef1)/coef1
mle.bias.f2<-(mle.coef.f2-coef2)/coef2
bayes1.bias.f1<-(bayes1.coef.f1-coef1)/coef1
bayes1.bias.f2<-(bayes1.coef.f2-coef2)/coef2
bayes6.bias.f1<-(bayes6.coef.f1-coef1)/coef1
bayes6.bias.f2<-(bayes6.coef.f2-coef2)/coef2
ml.pooled.bias.f1<-(ml.pooled.coef.f1-coef1)/coef1
ml.pooled.bias.f2<-(ml.pooled.coef.f2-coef2)/coef2
pooled.bias.f1<-(pooled.coef.f1-coef1)/coef1
pooled.bias.f2<-(pooled.coef.f2-coef2)/coef2
agdp.bias.f1<-(agdp.coef.f1-coef1)/coef1
agdp.bias.f2<-(agdp.coef.f2-coef2)/coef2
meta.bias.f1<-(meta.coef.f1-coef1)/coef1
meta.bias.f2<-(meta.coef.f2-coef2)/coef2
bayes_dp.random.bias.f1<-(bayes_dp.random.coef.f1-coef1)/coef1
bayes_dp.random.bias.f2<-(bayes_dp.random.coef.f2-coef2)/coef2
bayes_dp.75.bias.f1<-(bayes_dp.75.coef.f1-coef1)/coef1
bayes_dp.75.bias.f2<-(bayes_dp.75.coef.f2-coef2)/coef2
bayes_dp.50.bias.f1<-(bayes_dp.50.coef.f1-coef1)/coef1
bayes_dp.50.bias.f2<-(bayes_dp.50.coef.f2-coef2)/coef2
bayes_dp.25.bias.f1<-(bayes_dp.25.coef.f1-coef1)/coef1
bayes_dp.25.bias.f2<-(bayes_dp.25.coef.f2-coef2)/coef2
rbest.auto.bias.f1<-(rbest.auto.coef.f1-coef1)/coef1
rbest.auto.bias.f2<-(rbest.auto.coef.f2-coef2)/coef2
rbest.robust.bias.f1<-(rbest.robust.coef.f1-coef1)/coef1
rbest.robust.bias.f2<-(rbest.robust.coef.f2-coef2)/coef2
####COMBINE RESULTS####
coef.results<-cbind(mle.coef.f1,bayes1.coef.f1,bayes6.coef.f1,ml.pooled.coef.f1,pooled.coef.f1,agdp.coef.f1,meta.coef.f1,
bayes_dp.random.coef.f1,bayes_dp.75.coef.f1,bayes_dp.50.coef.f1,bayes_dp.25.coef.f1,rbest.auto.coef.f1,rbest.robust.coef.f1,
mle.coef.f2,bayes1.coef.f2,bayes6.coef.f2,ml.pooled.coef.f2,pooled.coef.f2,agdp.coef.f2,meta.coef.f2,
bayes_dp.random.coef.f2,bayes_dp.75.coef.f2,bayes_dp.50.coef.f2,bayes_dp.25.coef.f2,rbest.auto.coef.f2,rbest.robust.coef.f2)
se.results<-cbind(mle.se.f1,bayes1.se.f1,bayes6.se.f1,ml.pooled.se.f1,pooled.se.f1,agdp.se.f1,meta.se.f1,
bayes_dp.random.se.f1,bayes_dp.75.se.f1,bayes_dp.50.se.f1,bayes_dp.25.se.f1,rbest.auto.se.f1,rbest.robust.se.f1,
mle.se.f2,bayes1.se.f2,bayes6.se.f2,ml.pooled.se.f2,pooled.se.f2,agdp.se.f2,meta.se.f2,
bayes_dp.random.se.f2,bayes_dp.75.se.f2,bayes_dp.50.se.f2,bayes_dp.25.se.f2,rbest.auto.se.f2,rbest.robust.se.f2)
cover.results<-cbind(mle.cover.f1,mle.cover.f2,bayes1.cover.f1,bayes1.cover.f2,bayes6.cover.f1,bayes6.cover.f2,
ml.pooled.cover.f1,ml.pooled.cover.f2,pooled.cover.f1,pooled.cover.f2,agdp.cover.f1,agdp.cover.f2,meta.cover.f1,meta.cover.f2,
bayes_dp.random.cover.f1,bayes_dp.random.cover.f2,bayes_dp.75.cover.f1,bayes_dp.75.cover.f2,
bayes_dp.50.cover.f1,bayes_dp.50.cover.f2,bayes_dp.25.cover.f1,bayes_dp.25.cover.f2,rbest.auto.f1.cover,rbest.auto.f2.cover,
rbest.robust.f1.cover,rbest.robust.f2.cover)
bias.results<-cbind(mle.bias.f1,bayes1.bias.f1,bayes6.bias.f1,ml.pooled.bias.f1,pooled.bias.f1,agdp.bias.f1,meta.bias.f1,
bayes_dp.random.bias.f1,bayes_dp.75.bias.f1,bayes_dp.50.bias.f1,bayes_dp.25.bias.f1,rbest.auto.bias.f1,rbest.robust.bias.f1,
mle.bias.f2,bayes1.bias.f2,bayes6.bias.f2,ml.pooled.bias.f2,pooled.bias.f2,agdp.bias.f2,meta.bias.f2,
bayes_dp.random.bias.f2,bayes_dp.75.bias.f2,bayes_dp.50.bias.f2,bayes_dp.25.bias.f2,rbest.auto.bias.f2,rbest.robust.bias.f2)
simulation1.coef<-rbind(simulation1.coef, coef.results)
simulation1.se<-rbind(simulation1.se, se.results)
simulation1.cover<-rbind(simulation1.cover, cover.results)
simulation1.bias<-rbind(simulation1.bias, bias.results)
}
#CALCULATE MSE#
actual1<-rep(coef1,z)
actual2<-rep(coef2,z)
ml.mse1<-MSE(actual1,simulation1.coef[,1])
bayes1.mse1<-MSE(actual1,simulation1.coef[,2])
bayes6.mse1<-MSE(actual1,simulation1.coef[,3])
ml_pool.mse1<-MSE(actual1,simulation1.coef[,4])
bayes_pool.mse1<-MSE(actual1,simulation1.coef[,5])
agdp.mse1<-MSE(actual1,simulation1.coef[,6])
meta.mse1<-MSE(actual1,simulation1.coef[,7])
dp_random.mse1<-MSE(actual1,simulation1.coef[,8])
dp_75.mse1<-MSE(actual1,simulation1.coef[,9])
dp_50.mse1<-MSE(actual1,simulation1.coef[,10])
dp_25.mse1<-MSE(actual1,simulation1.coef[,11])
rbest_auto.mse1<-MSE(actual1,simulation1.coef[,12])
rbest_robust.mse1<-MSE(actual1,simulation1.coef[,13])
ml.mse2<-MSE(actual2,simulation1.coef[,14])
bayes1.mse2<-MSE(actual2,simulation1.coef[,15])
bayes6.mse2<-MSE(actual2,simulation1.coef[,16])
ml_pool.mse2<-MSE(actual2,simulation1.coef[,17])
bayes_pool.mse2<-MSE(actual2,simulation1.coef[,18])
agdp.mse2<-MSE(actual2,simulation1.coef[,19])
meta.mse2<-MSE(actual2,simulation1.coef[,20])
dp_random.mse2<-MSE(actual2,simulation1.coef[,21])
dp_75.mse2<-MSE(actual2,simulation1.coef[,22])
dp_50.mse2<-MSE(actual2,simulation1.coef[,23])
dp_25.mse2<-MSE(actual2,simulation1.coef[,24])
rbest_auto.mse2<-MSE(actual2,simulation1.coef[,25])
rbest_robust.mse2<-MSE(actual2,simulation1.coef[,26])
simulation1.mse<-cbind(ml.mse1,bayes1.mse1,bayes6.mse1,ml_pool.mse1,bayes_pool.mse1,agdp.mse1,meta.mse1,
dp_random.mse1,dp_75.mse1,dp_50.mse1,dp_25.mse1,rbest_auto.mse1,rbest_robust.mse1,
ml.mse2,bayes1.mse2,bayes6.mse2,ml_pool.mse2,bayes_pool.mse2,agdp.mse2,meta.mse2,
dp_random.mse2,dp_75.mse2,dp_50.mse2,dp_25.mse2,rbest_auto.mse2,rbest_robust.mse2)
simulation1.coef.mean<-colMeans(simulation1.coef)
simulation1.se.mean<-colMeans(simulation1.se)
simulation1.cover.mean<-colMeans(simulation1.cover)
simulation1.bias.mean<-colMeans(simulation1.bias)
simulation1.mse.mean<-colMeans(simulation1.mse)