rm(list=ls())

library(metaSEM)
library(lavaan)
library(semPlot)

# TSSEM1 Pool correlations together -------------------------------------
testdata <- readLowTriMat(file = "testdata.dat", no.var = 10)
testdata

nsize <- c(447,491,1076,600,414,1850,307)

sum(nsize)

#In the first run without constraining the heterogeneity indices,
# We observed 6 estimated of heterogeneity indices (Tau^2) which had hit the lower bound 1e-10; these have been constrained to 0
RE <- Diag(rep(0,45))
for (i in 1:45) {
  
  if (i %in% c(6,15,18,28,31,37)) {
    
    RE[i,i] = 0
  }
  else
  {
    RE[i,i] = paste("0.01*Tau2_",i,sep = "")
  }
  
}

View(RE)


tssem1pooled <- tssem1(Cov = testdata, n = nsize, method="REM", RE.type = "User", RE.constraints = RE)
summary(tssem1pooled)

# TSSEM2 ------------------------

mylavaan <- "

  F1 =~ C + I + J
  F2 =~ A + B
  F3 =~ D + E + G
  F4 =~ F + H
  HOF =~ F1 + F2 + F3 +F4
  
  HOF ~~ 1*HOF
  F1 ~~ F1
  F2 ~~ F2
  F3 ~~ F3
  F4 ~~ F4
  

  A ~~ A
  B ~~ B
  C ~~ C
  D ~~ D
  E ~~ E
  F ~~ F
  G ~~ G
  H ~~ H
  I ~~ I
  J ~~ J
  
"

plot(mylavaan, layout = "tree2")


names <- c("A","B","C","D","E","F","G","H","I","J")
myRAM <- lavaan2RAM(model = mylavaan, obs.variables = names, std.lv = T)


# Modify the S matrix to allow variance of lower order factors to be freely estimated
myRAM$S[11,11] <- "0.1*F1WITHF1"
myRAM$S[12,12] <- "0.1*F2WITHF2"
myRAM$S[13,13] <- "0.1*F3WITHF3"
myRAM$S[14,14] <- "0.1*F4WITHF4"


myRAM

tssem2model <- tssem2(tssem1.obj = tssem1pooled, RAM = myRAM , diag.constraints=T, intervals="LB")
summary(tssem2model)
tssem2model_rerun <- metaSEM::rerun(tssem2model)
summary(tssem2model_rerun)