## A function to convert rows into a 3x3 correlation matrix
create.matrix <- function(x, type=c(1, 2, 3)) {
  mat <- matrix(NA, ncol=3, nrow=3)
  diag(mat) <- 1
  type <- as.character(type)
  ## DC, M1, DV1
  ## 1: DC and M1
  ## 2: DC and DV1
  ## 3: M1 and DV1
  switch(type,
         "1" = mat[1, 2] <- mat[2, 1] <- unlist(x),
         "2" = mat[1, 3] <- mat[3, 1] <- unlist(x),
         "3" = mat[2, 3] <- mat[3, 2] <- unlist(x))
  mat
}

varlist <- c("DC", "M1", "DV1")

my.cor <- lapply(split(my.df, seq(nrow(my.df))),
                 function(x, y) create.matrix(x["Effect_size"], x["Type_of_Association"]))
my.cor <- lapply(my.cor, function(x) {dimnames(x) <- list(varlist, varlist); x}  )
names(my.cor) <- my.df$StudyID

## Correlation matrices in the analysis
my.cor

## Sample sizes
my.n <- my.df$ni
my.n
## Number of studies in each cell
pattern.na(my.cor.imputed, show.na = FALSE)
## Total sample sizes in each cell
pattern.n(my.cor.imputed, my.n)
## Stage 1 analysis: find an average correlation matrix
stage1 <- tssem1(my.cor.imputed, my.n)
summary(stage1)
## Average correlation matrix
meanR <- vec2symMat(coef(stage1, select = "fixed"), diag = FALSE)
dimnames(meanR) <- list(varlist, varlist)
meanR
## Absolute heterogeneity variance: tau^2
tau2 <- vec2symMat(coef(stage1, select = "random"), diag = FALSE)
dimnames(tau2) <- list(varlist, varlist)
tau2
## Relative heterogeneity index: I^2
I2 <- vec2symMat(summary(stage1)$I2.values[, "Estimate"], diag = FALSE)
dimnames(I2) <- list(varlist, varlist)
I2
## Proposed model
model <- "DV1 ~ c*DC + b*M1
          M1 ~ a*DC
          DC ~~ 1*DC"
plot(model, color="yellow")
RAM1 <- lavaan2RAM(model, obs.variables = varlist)
RAM1
## Stage 2 analysis: fit the path model
stage2 <- tssem2(stage1, RAM=RAM1, intervals.type = "LB",
                 mx.algebras = list(Indirect=mxAlgebra(a*b, name="Indirect"),
                                    Direct=mxAlgebra(c, name="Direct")))
summary(stage2)
plot(stage2, color="yellow")
# svg("Fig4.svg", width=5, height=5)
# plot(stage2, col="yellow")
# dev.off()

## Testing the hypothesis c = a*b
stage2b <- tssem2(stage1, RAM=RAM1, intervals.type = "LB",
                  mx.algebras = list(Indirect=mxAlgebra(a*b, name="Indirect"),
                                     Direct=mxAlgebra(c, name="Direct")),
                  run=FALSE)
## Add a constraint on c=a*b
## Rerun to remove errors
stage2b <- mxModel(stage2b, mxConstraint(c==a*b, name="constraint"))
stage2b <- mxRun(stage2b)
stage2b <- mxTryHard(stage2b, extraTries = 100)
stage2b <- mxRun(stage2b, intervals = TRUE)
summary(stage2b)
anova(stage2$mx.fit, stage2b)