##Model 1: PSY --> MOV;PSY ---> PAIN; PAIN-->MOV install.packages("OpenMx") install.packages("metaSEM") install.packages("semPlot") ## Load the library library(metaSEM) ## Loading required package: OpenMx ## To take full advantage of multiple cores, use: mxOption(NULL, 'Number of Threads', parallel::detectCores()) #now Sys.setenv(OMP_NUM_THREADS=parallel::detectCores()) #before library(OpenMx) library(OpenMx) ## "SLSQP" is set as the default optimizer in OpenMx. ## mxOption(NULL, "Gradient algorithm") is set at "central". mxOption(NULL, "Gradient algorithm") ## mxOption(NULL, "Optimality tolerance") is set at "6.3e-14". mxOption(NULL, "Optimality tolerance") ## mxOption(NULL, "Gradient iterations") is set at "2". mxOption(NULL, "Gradient iterations") ## Load the library for figures library(semPlot) ## Create the matrix ## Studiestest10 <- structure(list(data = structure(list(`Anderson2013` = structure(c(1,-0.002,-0.078,-0.002,1,0.308,-0.078,0.308,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity"))), `Basler2008` = structure(c(1,-0.2,NA,-0.2,1,NA,NA,NA,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity"))), `Christe2016` = structure(c(1,0.162,0.11,0.162,1,0.05,0.11,0.05,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity"))), `Christe2017` = structure(c(1,-0.486,0.198,-0.486,1,-0.015,0.198,-0.015,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity"))), `Courbalay2017` = structure(c(1,0.05,0.12,0.05,1,0.38,0.12,0.38,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity"))), `Demoulin2013` = structure(c(1,0.061,0.116,0.061,1,-0.1735,0.116,-0.1735,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity"))), `Dubois2016` = structure(c(1,-0.07,0.06,-0.07,1,0.34,0.06,0.34,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity"))), `George2006` = structure(c(1,-0.25,NA,-0.25,1,0.135,NA,0.135,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.22,-0.23,0.22,1), .Dim = c(3L, 3L), .Dimnames = list(c("Movement", "Psychological_factor", "Pain_intensity"), c("Movement", "Psychological_factor", "Pain_intensity")))), .Names = c("Anderson2013", "Basler2008", "Christe2016", "Christe2017", "Courbalay2017", "Demoulin2013", "Dubois2016", "George2006", "Grotle2004a", "Grotle2004b")), n = c(109, 99, 10, 10, 17, 50, 36, 63, 123, 233), 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") ######################################################################### ### Fixed-effects model ### ## Mike: There is an error message (OpenMx status1: 6). ## The reason seems to be related to Tau2_2_2 and Tau2_3_3, which ## are too small. Since the p value of the homogeneity test is .57, ## it seems reasonable to apply the fixed-effects model. ## Fixed-effects model fix1 <- tssem1(Studiestest10$data, Studiestest10$n, method="FEM") summary(fix1) ## stage 2 ## ## Regression coefficients A1 <- create.mxMatrix(c(0, "0.2*PsytoMov", "0.2*PaitoMov", 0, 0, 0, 0, "0.2*PsytoPai", 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 ## Mike: Since Psy is an independent variable, its variance should be fixed at 1.0. ## 0.2 is just the starting value. You may use other sensible values, e.g., 0.3, 0.4... S1 <- create.mxMatrix(c("0.2*e_Mov",0,0, 0, 1, 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) fix2 <- tssem2(fix1, Amatrix=A1, Smatrix=S1, intervals.type="z") summary(fix2) plot(fix2) ## option 2 ## Mike: it works fine now. fix2 <- tssem2(fix1, Amatrix=A1, Smatrix=S1, intervals.type="LB", diag.constraints=TRUE) summary(fix2) plot(fix2) ## option 3 (Cheung/Hagger) ## Mike: I think that you want PsytoPai*PaitoMov for the indirect effect. ## It works fine now. fix2 <- tssem2(fix1, Amatrix=A1, Smatrix=S1, diag.constraints=TRUE, intervals.type="LB", mx.algebras=list(Ind=mxAlgebra(PsytoPai*PaitoMov, name="Ind"))) summary(fix2) plot(fix2) ## Refit again to have all the LBCI ## Is this useful ??? fix2 <- rerun(fix2) summary(fix2) ## Display the A matrix mxEval(Amatrix, fix2$mx.fit) ## Compute R2 intention ## I don't understand this value mxEval(1-Smatrix, fix2$mx.fit)[1,1] ###################################################################################### ## if necessary for the figure ## Convert the model to semPlotModel object my.plot <- meta2semPlot(fix2) ## Plot the parameter estimates semPaths(my.plot, whatLabels="est", nCharNodes=10, color="yellow", layout="circle2", edge.label.cex = 0.8)