##Load MetaSEM and OpenMx libraries
library(metaSEM)
library(OpenMx)
library(semPlot)
library(purrr)

data <- readLowTriMat(file="cormat3.dat", no.var=10, na.strings = "NA")
## Add variable names into the data set
varnames <- c("Aan",	"Aav",	"Par",	"Hal",	"Emoreg",	"BaO",	"BaV",	
              "RC",	"NSE",	"Ase")


nvar <- 10

labels <- list(varnames,varnames)

##Specify sample sizes
n <-  c(60,37,55,44, 176, 54, 180, 361, 40, 503,302,113)



#label studies
Ascone <- data[1]
Castilho <- data[2]
Pilton <- data[3]
Robson <- data[4]
Wickham <- data[5]
Gumley <- data[6]
Cole <- data[7]
Hart <- data[8]
Ascone2 <- data[9]
Pickering <- data[10]
Udachina <- data[11]
Wickham2 <- data[12]

#### Anxiety 	Beliefs about Others	Hallucinations ####
Model <- rbind (Wickham, Wickham2, Pickering, Cole)
Model_n <- c(176,113,503,180)

Model_A <- list(Model[[1]][c(1,6,4),c(1,6,4)],Model[[2]][c(1,6,4),c(1,6,4)], Model[[3]][c(1,6,4),c(1,6,4)],Model[[4]][c(1,6,4),c(1,6,4)])
# put NA on diagonal if variable is missing
for (i in 1:length(Model_A)){
  for (j in 1:nrow(Model_A[[i]])){
    if (sum(is.na(Model_A[[i]][j,]))==nvar-1)
    {Model_A[[i]][j,j] <- NA}
  }}
pattern.na(Model_A, show.na=FALSE)

nvar=3
# put NA on diag for var with least present correlations
for (i in 1:length(Model_A)){
  for (j in 1:nrow(Model_A[[i]])){
    for (k in 1:nvar){
      if (is.na(Model_A[[i]][j,k])==TRUE
          &is.na(Model_A[[i]][j,j])!=TRUE
          &is.na(Model_A[[i]][k,k])!=TRUE){
        if(sum(is.na(Model_A[[i]])[j,])>sum(is.na(Model_A[[i]])[k,]))
        {Model_A[[i]][k,k] <- NA}
        if(sum(is.na(Model_A[[i]])[j,])<=sum(is.na(Model_A[[i]])[k,]))
        {Model_A[[i]][j,j] <- NA}
      }}}}

Model4 <- tssem1(Cov=Model_A, Model_n, method="REM", RE.type="Diag") # only estimate the study level variance (RE)
summary(Model4)

Model4 <- rerun(Model4, autofixtau2 = TRUE)
summary(Model4)


## outputs pooled covariance matrix (the parameter estimates) extracted
vcov(Model4)



###### Prepare models for stage 2 analysis
## A1: asymmetric matrix (regression coefficients)
A1 <-  create.mxMatrix(
  c( 0,0,0,
     "0.1*b1",0,0,
     "0.1*b2","0.1*b3",0),
  type = "Full",
  nrow = 3,
  ncol = 3,
  byrow = TRUE)

A1
# b1 is predictor/mediator, b2 is predictor/outcome, b3 mediator/outcome


##create a symmetric matrix (variance covariance matrix among variables)

S1 <- create.mxMatrix(
  c("1*p11",
    0,"1*p22",
    0,0,"1*p33"),
  type="Symm",
  byrow = TRUE, as.mxMatrix=FALSE)


## Conduct stage 2 structural equation model

random4 <- tssem2(Model4, Amatrix=A1, Smatrix=S1, diag.constraints=TRUE, 
                  intervals="LB", mx.algebras=list(Ind=mxAlgebra(b1*b3, name="Ind")))

summary(random4)

my.plot4 <- meta2semPlot(random4)
semPaths(my.plot4, whatLabels="est", layout = "tree2")