##Model 2: PSY --> MOV;PAIN --> PSY; PAIN-->MOV

## Load the library
library(metaSEM)

## Create the matrix ##

Studiestest10 <-

structure(list(data = structure(list(`Trost2012` = structure(c(1,-0.38,-0.11,-0.38,1,0.14,-0.11,0.14,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Thomas2008A` = structure(c(1,-0.54,-0.24,-0.54,1,0.28,-0.24,0.28,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Thomas2008B` = structure(c(1,-0.41,-0.24,-0.41,1,0.28,-0.24,0.28,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Thomas2008C` = structure(c(1,-0.35,-0.24,-0.35,1,0.28,-0.24,0.28,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Watson1997A` = structure(c(1,-0.28,-0.2,-0.28,1,0.28,-0.2,0.28,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Watson1997B` = structure(c(1,-0.27,-0.2,-0.27,1,0.28,-0.2,0.28,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Alschuler2009` = structure(c(1,-0.45,-0.11,-0.45,1,0.28,-0.11,0.28,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Grotle2004A` = structure(c(1,-0.08,-0.19,-0.08,1,0.13,-0.19,0.13,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Grotle2004B` = structure(c(1,-0.11,-0.23,-0.11,1,0.21,-0.23,0.21,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity"))), `Kim2017` = structure(c(1,-0.03,-0.02,-0.03,1,0.25,-0.02,0.25,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", 

"Psychological_factor", "Pain_intensity"), c("Movement", 

"Psychological_factor", "Pain_intensity")))), .Names = c("Trost2012", "Thomas2008A", "Thomas2008B", "Thomas2008C", "Watson1997A", "Watson1997B", "Alschuler2009", "Grotle2004A", "Grotle2004B", "Kim2017")), n = c(51, 36, 36, 36, 36, 36, 76, 123, 233, 30), Questionnaire = c("PASS", "PASS", "TSK", "PCS", "FABQ", "self_efficacy", "TSK", "FABQ", "FABQ", "FABQ")), .Names = c("data", "n", 

"Questionnaire"))

## show the studies ##

head(Studiestest10$data)

## Display the sample sizes
Studiestest10$n

## Variables used in the analysis
var.names <- c("Mov", "Psy", "Pai")

#######################################################################


### Stage 1 from Hagger (give the same results as version from Cheung)
random1 <- tssem1(Studiestest10$data, Studiestest10$n, method="REM", RE.type="Diag", acov="weighted")
summary(random1)

## Extract the fixed-effects estimates
(est_fixed <- coef(random1, select="fixed"))

## Convert the estimated vector to a symmetrical matrix
## where the diagonals are fixed at 1 (for a correlation matrix)
vec2symMat(est_fixed, diag=FALSE)


## stage 2
## Regression coefficients --> why 0.2 (same as Hagger) ???
A1 <- create.mxMatrix(c(0, "0.2*PsytoMov", "0.2*PaitoMov",
0, 0, "0.2*PaitoPsy",
0, 0, 0),
type="Full", byrow=TRUE, ncol=3, nrow=3,
as.mxMatrix=FALSE)

## This step is not necessary but it is useful for inspecting the model.
dimnames(A1)[[1]] <- dimnames(A1)[[2]] <- c("Mov","Psy","Pai")
A1

## Covariance matrix among the variables (CAUTION: NOT SURE FOR THIS STEP --> why 0.2 ???)
S1 <- create.mxMatrix(c("0.2*e_Mov",0,0,
0, "0.2*e_Psy", 0,
0, 0, "0.2*e_Pai"),
type="Full", byrow=TRUE, ncol=3, nrow=3, as.mxMatrix=FALSE)

## This step is not necessary but it is useful for inspecting the model.
dimnames(S1)[[1]] <- dimnames(S1)[[2]] <- c("Mov","Psy","Pai")
S1


## Stage 2 analysis: different option are in the website from Cheung and in its paper


## option 1 (works, from Cheung)
random2 <- tssem2(random1, Amatrix=A1, Smatrix=S1, intervals.type="z")
summary(random2)


## option 2 (not working, from Cheung)
random2 <- tssem2(random1, Amatrix=A1, Smatrix=S1, intervals.type="LB",
diag.constraints=TRUE)

## option 3 (Cheung/Hagger) --> works like option 1; but cannot have the indirect effect

random2 <- tssem2(random1, Amatrix=A1, Smatrix=S1, diag.constraints=TRUE,
intervals.type="LB", mx.algebras( Ind=mxAlgebra(PaitoPsy*PaitoMov, name="Ind"))


## Refit again to have all the LBCI
random2 <- rerun(random2)
summary(random2)

## Display the A matrix
mxEval(Amatrix, random2$mx.fit)

## Compute R2 intention
mxEval(1-Smatrix, random2$mx.fit)[1,1]

## Draw the output
library(semPlot)

## Convert the model to semPlotModel object
my.plot <- meta2semPlot(random2)

## Plot the parameter estimates
semPaths(my.plot, whatLabels="est", nCharNodes=10, color="yellow", 
         layout="circle2", edge.label.cex = 0.8)

#######################################################################################################

## Conduct stage 2 structural equation model (Hagger 2018) -  not working
##Indirect & Direct effects inventory:
##Ind1= Indirect effects of Psy on Mov through pain
random2 <- tssem2(random1, Amatrix=A1 Smatrix=S1, diag.constraints=FALSE,
                  intervals.type="LB",
                  mx.algebras=list
                 (Ind1 =mxAlgebra(PsytoPai*PaitoMov,name="Ind1"))