#Fitting a Mediation model
##################Creating the Correlation Matrice for all the variables we have###########
##Load the Data 
my.data <- read_dta("MASEM2.dta")  
 
#Load MetaSEM package
library(metaSEM)

#load OpenMx and change optimizer to NPSOL from the default SLSQP loaded by the metaSEM. The NPSOL allows for faster optimization

library(OpenMx)

head(my.data)
colnames(my.data)
dim(my.data)

# locate studies with information on at least 1 correlation, and no missing sample size
keepstudy <- rowSums(is.na(my.data[,5:82]))!=78 & is.na(my.data$NUMBER.OF.BANKS)==FALSE
# keep only the studies with information
data <- my.data[keepstudy,]

# check data
head(data)
length(data)
summary(data)

# varnames and labels
nvar <- 13
varnames <- c("A","B","C","D","E","F","G","H","I","J","K","L","M")
labels <- list(varnames,varnames)
# create list with correlation matrices
cordat <- list()
for (i in 1:nrow(data)){	
  cordat[[i]] <- vec2symMat(as.matrix(data[i,5:82]),diag = FALSE)
  dimnames(cordat[[i]]) <- labels
}
# put NA on diagonal of correlation matrix if variable is missing
for (i in 1:length(cordat)){
  for (j in 1:nrow(cordat[[i]])){	
    if (sum(is.na(cordat[[i]][j,]))==nvar-1) 
    {cordat[[i]][j,j] <- NA}
  }}

# show number of studies per correlation coefficient
pattern.na(cordat, show.na = FALSE)

# show total N per correlation coefficient
pattern.n(cordat, data$NUMBER.OF.BANKS)

#### Running of Stage1 model
stage1random <- tssem1(Cov=cordat, n=data$NUMBER.OF.BANKS, method="REM", RE.type="Diag")
summary(stage1random)

##Model requires a rerun
# Rerun model for convergence
stage1_rerun <- metaSEM::rerun(stage1random)
summary(stage1_rerun)

####Model
model10 <- "M ~ c*E +f*F+i*G+l*H+s*J+t*K+x*L
RISK ~ a*E +d*F+g*G+j*H
PERF ~ b*I
E ~ k*A + w*B + b1*C + e*D+s1*J+t1*K+x1*L
F ~ h*A + m*B + n*C + o*D+s2*J+t2*K+x2*L
G ~ p*A + q*B + r*C + u*D+s3*J+t3*K+x3*L
H ~ v*A + w1*B + y*C + z*D+s4*J+t4*K+x4*L
T ~ a1*A +d1*B +g1*C +j1*D

#Covariance
E ~~ F
E ~~ G
E ~~ H
E ~~ J
E ~~ K
E ~~ L
F ~~ G
F ~~ H
F ~~ J
F ~~ K
F ~~ L
G ~~ H
G ~~ H
G ~~ H
G ~~ H
H ~~ J
H ~~ K
H ~~ L
J ~~ K
J ~~ L
K ~~ L

#variance
A ~~ 1*A
B ~~ 1*B
C ~~ 1*C
D ~~ 1*D
E ~~ 1*E
F ~~ 1*F
G ~~ 1*G
H ~~ 1*H
J ~~ 1*J
K ~~ 1*K
L ~~ 1*L"

plot(model10)

## Convert the lavaan syntax to RAM specification used in metaSEM
RAM1 <- lavaan2RAM(model10, obs.variables= varnames)
RAM1

## Request the likelihood-based confidence interval
## Indirect effect: ind = a*b
## Direct effect: dir = c
tssem.fit <- tssem2(stage1_rerun, RAM=RAM1, intervals.type = "LB",
               mx.algebras = list(ind_E=mxAlgebra(a*b, name="ind_E"),ind_F=mxAlgebra(d*b, name="ind_F"),ind_G=mxAlgebra(g*b, name="ind_G"),ind_H=mxAlgebra(j*b, name="ind_H"),
               Dir_BINDP=mxAlgebra(c, name="Dir_E"),Dir_DUA=mxAlgebra(f, name="Dir_F"),Dir_G=mxAlgebra(i, name="Dir_G"),Dir_H=mxAlgebra(l, name="Dir_H"),Dir_J=mxAlgebra(s, name="Dir_J"),Dir_K=mxAlgebra(t, name="Dir_K"),Dir_L=mxAlgebra(x, name="Dir_L")))
               summary(tssem.fit)