##Load the Data
my.data <- read.cvs("dataLeo.csv")

#Load MetaSEM package
library(metaSEM)

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

# locate studies with information on at least 1 correlation, and no missing sample size
keepstudy <- rowSums(is.na(my.data[,4:48]))!=45 & 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 <- 10
varnames <- c("BS","IND","WOB","DUA","BMEET","DIROWN","BNKRISK","BSIZE","LEV","PERF")
labels <- list(varnames,varnames)


# create list with correlation matrices
cordat <- list()
for (i in 1:nrow(data)){	
  cordat[[i]] <- vec2symMat(as.matrix(data[i,4:48]),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)

## Average correlation matrix under a random-effects model
averageR <- vec2symMat(coef(stage1random, select="fixed"), diag = FALSE)
dimnames(averageR) <- list(obs.vars1, obs.vars1)
averageR

## Heterogeneity variances of the random-effects
coef(stage1random, select="random")

######Full model
model10 <- "PERF ~ c*BS + b*BNKRISK+f*IND+i*WOB+l*DUA+o*BMEET+r*DIROWN+BSIZE+LEV
BNKRISK ~ a*BS
BNKRISK ~ d*IND
BNKRISK ~ g*WOB
BNKRISK ~ j*DUA
BNKRISK ~ m*BMEET
BNKRISK ~ p* DIROWN
BS ~~ 1*BS
IND ~~ 1*IND
WOB ~~ 1*WOB
DUA ~~ 1*DUA
BMEET ~~ 1*BMEET
DIROWN ~~ 1*DIROWN
BSIZE ~~ 1*BSIZE
LEV ~~ 1*LEV
## Define indirect and direct effects. 
Ind_BS := a*b
Ind_IND := d*b
Ind_WOB := g*b
Ind_DUA := j*b
Ind_BMEET := m*b
Ind_DIROWN := p*b
Dir_BS := c
Dir_IND := f
Dir_WOB := i
Dir_DUA := l
Dir_BMEET := o
Dir_DIROWN := r"

plot(model10)
## Convert the lavaan syntax to RAM specification used in metaSEM
RAM1 <- lavaan2RAM(model1, 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_BS=mxAlgebra(a*b, name="ind_BS"),ind_IND=mxAlgebra(d*b, name="ind_IND"),ind_WOB=mxAlgebra(g*b, name="ind_WOB"),ind_DUA=mxAlgebra(j*b, name="ind_DUA"),ind_BMEET=mxAlgebra(m*b, name="ind_BMEET"),ind_DIROWN=mxAlgebra(p*b, name="ind_DIROWN"),
               Dir_BS=mxAlgebra(c, name="Dir_BS"),Dir_IND=mxAlgebra(f, name="Dir_IND"),Dir_WOB=mxAlgebra(i, name="Dir_WOB"),Dir_DUA=mxAlgebra(l, name="Dir_DUA"),Dir_BMEET=mxAlgebra(o, name="Dir_BMEET"),Dir_DIROWN=mxAlgebra(r, name="Dir_DIROWN")))
               summary(tssem.fit)