## Load the libraries for TSSEM
library("metaSEM")
library("semPlot")
library("lavaan")

## Display some data
test <- structure(list(data=structure(list(`DualéEDONIS` = structure(c(1,-0.053, 0.068,-0.053,1,0.028,0.068,0.028,1),
                                                                     .Dim = c(3L, 3L),
                                                                     .Dimnames = list(c("Cata_t1", "Pain_t1", "Pain_t2"),
                                                                                      c("Cata_t1", "Pain_t1", "Pain_t2"))),
                                           `Opstelten` = structure(c(1,-0.067,0.019,-0.067,1,-0.078,0.019,-0.078,1),
                                                                   .Dim = c(3L, 3L),
                                                                   .Dimnames = list(c("Cata_t1", "Pain_t1", "Pain_t2"),
                                                                                    c("Cata_t1", "Pain_t1", "Pain_t2"))),
                                           `Postsurgical` = structure(c(1,-0.045,0.126,-0.045,1,-0.128,0.126,-0.128,1),
                                                                      .Dim = c(3L, 3L),
                                                                      .Dimnames = list(c("Cata_t1", "Pain_t1", "Pain_t2"),
                                                                                       c("Cata_t1", "Pain_t1", "Pain_t2")))),
                                      .Names = c("DualéEDONIS", "Opstelten", "Postsurgical")),
                       n=c(1436,124,74),
                       Questionnaire = c("PCS", "PCL", "CSQ"),
                       NumberMen = c(400,62,10)),
                  .Names = c("data", "n", "Quest", "NumberMen"))

#### Fixed-effects model
## Stage 1 analysis
fixed1 <- tssem1(test$data, test$n, method="FEM")
summary(fixed1)

coef(fixed1)

##Model 1
## Using lavaan's syntax
model2 <- "## Regression model
           Pain_t2~Pain_t1+Cata_t1
           ## Variances fixed at 1
           Pain_t1 ~~ 1*Pain_t1
           Cata_t1 ~~ 1*Cata_t1"
plot(model2)

RAM2 <- lavaan2RAM(model2)
RAM2

## Stage 2 analysis
fixed2 <- tssem2(fixed1, RAM=RAM2, diag.constraints=TRUE, 
                 intervals="LB")

summary(fixed2)

test2 <- tssem2(fixed1, RAM=RAM2)
summary(test2)

## First stage of TSSEM: Random-effects model
rand1 <- tssem1(test$data, test$n, method="REM")

## Show the summary
summary(rand1)

## Plot the effect sizes
plot.meta(rand1, axis.labels=c("Cor(LS,JS)", "Cor(LS,JA)", "Cor(JS,JA)"), 
          study.min.cex=0.5, main="Effect sizes and their CIs")

## Build a model with indirect effect
## A1: Regression coefficients; 
## parameters are labelled with "a", "b", and "c"
## S1: Variance-covariance of parameters
A1 <- create.mxMatrix(c(0,"0.2*a","0.2*b",
                        0,0,"0.2*c",
                        0,0,0), 
                      type="Full", ncol=3, nrow=3, as.mxMatrix=FALSE, byrow=TRUE)
dimnames(A1) <- list(c("Cata_t1", "Pain_t1", "Pain_t2"),
                     c("Cata_t1", "Pain_t1", "Pain_t2"))
A1

S1 <- create.mxMatrix(c("0.2*Error_y","0.2*Error_m",1), type="Diag", as.mxMatrix=FALSE)
dimnames(S1) <- list(c("Cata_t1", "Pain_t1", "Pain_t2"),
                     c("Cata_t1", "Pain_t1", "Pain_t2"))
S1

## Second stage of TSSEM: Random-effects model
## Estimate direct and indirect effects
rand2 <- tssem2(rand1, Amatrix=A1, Smatrix=S1, 
                diag.constraints=TRUE, intervals.type="LB",
                mx.algebras=list( indirect=mxAlgebra(a*b, name="indirect"), 
                                  direct=mxAlgebra(c, name="direct") ))

## Show the summary
summary(rand2)

## Plot the model
my.plot1 <- meta2semPlot(rand2)
semPaths(my.plot1, whatLabels="est", edge.label.cex=1.5, sizeMan=8, color="yellow")