library("OpenMx")
library("metaSEM")

(x1 <- vec2symMat(c(1,0.67,0.67,0.63,-0.56,1,0.77,0.66,-0.58,1,0.62,-0.58,1,-0.57,1)))
(x2 <- vec2symMat(c(1,NA,0.522,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x3 <- vec2symMat(c(1,0.24,NA,NA,-0.13,1,NA,NA,-0.34,1,NA,NA,1,NA,1)))
(x4 <- vec2symMat(c(1,0.25,NA,NA,-0.14,1,NA,NA,-0.35,1,NA,NA,1,NA,1)))
(x5 <- vec2symMat(c(1,NA,0.36,0.34,-0.18,1,NA,NA,NA,1,0.04,-0.1,1,-0.32,1)))
(x6 <- vec2symMat(c(1,0.44,0.57,0.47,-0.25,1,0.62,0.53,-0.62,1,0.3,-0.62,1,-0.39,1)))
(x7 <- vec2symMat(c(1,0.7,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x8 <- vec2symMat(c(1,NA,0.668,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x9 <- vec2symMat(c(1,NA,0.43,NA,-0.29,1,NA,NA,NA,1,NA,-0.28,1,NA,1)))
(x10 <- vec2symMat(c(1,NA,NA,0.6,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x11 <- vec2symMat(c(1,0.29,0.29,NA,NA,1,0.48,NA,NA,1,NA,NA,1,NA,1)))
(x12 <- vec2symMat(c(1,NA,0.54,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x13 <- vec2symMat(c(1,0.453,0.7,0.72,-0.618,1,0.586,0.536,-0.51,1,0.68,-0.665,1,-0.623,1)))
(x14 <- vec2symMat(c(1,0.42,NA,0.49,NA,1,NA,0.74,NA,1,NA,NA,1,NA,1)))
(x15 <- vec2symMat(c(1,0.218,NA,0.437,-0.171,1,NA,0.376,-0.644,1,NA,NA,1,-0.336,1)))
(x16 <- vec2symMat(c(1,0.561,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x17 <- vec2symMat(c(1,0.564,0.456,0.413,-0.259,1,0.539,0.472,-0.451,1,0.315,-0.274,1,-0.351,1)))
(x18 <- vec2symMat(c(1,NA,NA,0.369,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x19 <- vec2symMat(c(1,0.62,NA,NA,-0.64,1,NA,NA,-0.69,1,NA,NA,1,NA,1)))
(x20 <- vec2symMat(c(1,NA,0.413,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x21 <- vec2symMat(c(1,0.18,0.66,NA,-0.18,1,0.38,NA,-0.9,1,NA,-0.43,1,NA,1)))
(x22 <- vec2symMat(c(1,NA,0.529,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x23 <- vec2symMat(c(1,0.54,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x24 <- vec2symMat(c(1,0.309,0.202,NA,NA,1,0.508,NA,NA,1,NA,NA,1,NA,1)))
(x25 <- vec2symMat(c(1,0.31,0.4,NA,NA,1,0.42,NA,NA,1,NA,NA,1,NA,1)))
(x26 <- vec2symMat(c(1,0.234,NA,0.289,-0.238,1,NA,0.301,-0.196,1,NA,NA,1,-0.271,1)))
(x27 <- vec2symMat(c(1,0.195656566,NA,NA,-0.231157025,1,NA,NA,-0.328838384,1,NA,NA,1,NA,1)))
(x28 <- vec2symMat(c(1,0.230833333,NA,NA,-0.173461538,1,NA,NA,-0.524405594,1,NA,NA,1,NA,1)))
(x29 <- vec2symMat(c(1,NA,0.28,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x30 <- vec2symMat(c(1,0.62,0.54,0.33,NA,1,0.54,0.65,NA,1,0.34,NA,1,NA,1)))
(x31 <- vec2symMat(c(1,NA,0.454,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x32 <- vec2symMat(c(1,NA,0.832,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x33 <- vec2symMat(c(1,NA,NA,NA,-0.303,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x34 <- vec2symMat(c(1,0.45,0.77,0.38,-0.51,1,0.49,0.74,-0.72,1,0.52,-0.54,1,-0.81,1)))
(x35 <- vec2symMat(c(1,0.338,NA,0.396,-0.119,1,NA,0.53,-0.583,1,NA,NA,1,-0.52,1)))
(x36 <- vec2symMat(c(1,NA,0.655,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x37 <- vec2symMat(c(1,0.408,0.624,NA,NA,1,0.461,NA,NA,1,NA,NA,1,NA,1)))
(x38 <- vec2symMat(c(1,0.092,NA,NA,-0.234,1,NA,NA,-0.415,1,NA,NA,1,NA,1)))
(x39 <- vec2symMat(c(1,NA,NA,0.766,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x40 <- vec2symMat(c(1,NA,0.533,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x41 <- vec2symMat(c(1,NA,NA,NA,-0.282,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x42 <- vec2symMat(c(1,NA,NA,0.445,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x43 <- vec2symMat(c(1,0.613,0.51,NA,NA,1,0.739,NA,NA,1,NA,NA,1,NA,1)))
(x44 <- vec2symMat(c(1,NA,NA,0.61,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x45 <- vec2symMat(c(1,NA,NA,0.77,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x46 <- vec2symMat(c(1,0.273,0.5148,NA,NA,1,0.4016,NA,NA,1,NA,NA,1,NA,1)))
(x47 <- vec2symMat(c(1,NA,0.743,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x48 <- vec2symMat(c(1,0.304,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x49 <- vec2symMat(c(1,0.456,NA,NA,-0.31,1,NA,NA,-0.244,1,NA,NA,1,NA,1)))
(x50 <- vec2symMat(c(1,0.537,NA,NA,-0.398,1,NA,NA,-0.36,1,NA,NA,1,NA,1)))
(x51 <- vec2symMat(c(1,0.376,0.502,0.633,-0.363,1,0.51,0.434,-0.447,1,0.559,-0.592,1,-0.535,1)))
(x52 <- vec2symMat(c(1,0.47,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x53 <- vec2symMat(c(1,NA,NA,0.619,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x54 <- vec2symMat(c(1,0.15,NA,0.17,NA,1,NA,0.13,NA,1,NA,NA,1,NA,1)))
(x55 <- vec2symMat(c(1,0.15,NA,0.28,NA,1,NA,0.23,NA,1,NA,NA,1,NA,1)))
(x56 <- vec2symMat(c(1,0.57,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x57 <- vec2symMat(c(1,NA,NA,NA,-0.18,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x58 <- vec2symMat(c(1,0.4,0.34,NA,NA,1,0.55,NA,NA,1,NA,NA,1,NA,1)))
(x59 <- vec2symMat(c(1,0.386,0.63,0.384,NA,1,0.531,0.549,NA,1,0.367,NA,1,NA,1)))
(x60 <- vec2symMat(c(1,NA,NA,0.43,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x61 <- vec2symMat(c(1,0.51,0.66,NA,-0.61,1,0.76,NA,-0.57,1,NA,-0.66,1,NA,1)))
(x62 <- vec2symMat(c(1,0.499,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x63 <- vec2symMat(c(1,0.602,0.716,NA,NA,1,0.622,NA,NA,1,NA,NA,1,NA,1)))
(x64 <- vec2symMat(c(1,NA,NA,0.18,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x65 <- vec2symMat(c(1,NA,0.6414,0.0027,-0.1386,1,NA,NA,NA,1,0.0459,-0.1925,1,-0.629,1)))
(x66 <- vec2symMat(c(1,NA,NA,0.353,-0.385,1,NA,NA,NA,1,NA,NA,1,-0.478,1)))
(x67 <- vec2symMat(c(1,0.51,0.68,0.35,NA,1,0.47,0.5,NA,1,0.48,NA,1,NA,1)))
(x68 <- vec2symMat(c(1,0.62,0.77,NA,NA,1,0.63,NA,NA,1,NA,NA,1,NA,1)))
(x69 <- vec2symMat(c(1,NA,0.53,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x70 <- vec2symMat(c(1,0.28,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x71 <- vec2symMat(c(1,0.57,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x72 <- vec2symMat(c(1,0.32,0.561,NA,NA,1,0.505,NA,NA,1,NA,NA,1,NA,1)))
(x73 <- vec2symMat(c(1,0.795,0.761,NA,-0.671,1,0.719,NA,-0.688,1,NA,-0.633,1,NA,1)))
(x74 <- vec2symMat(c(1,NA,0.532,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x75 <- vec2symMat(c(1,NA,0.522,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x76 <- vec2symMat(c(1,0.493,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x77 <- vec2symMat(c(1,NA,0.31,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x78 <- vec2symMat(c(1,0.45,0.71,NA,NA,1,0.66,NA,NA,1,NA,NA,1,NA,1)))


my.df6 <- list(x1, x2, x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22,x23,x24,x25,x26,x27,x28,x29,x30,x31,x32,x33,x34,x35,x36,x37,x38,x39,x40,x41,x42,x43,x44,x45,x46,x47,x48,x49,x50,x51,x52,x53,x54,x55,x56,x57,x58,x59,x60,x61,x62,x63,x64,x65,x66,x67,x68,x69,x70,x71,x72,x73,x74,x75,x76,x77,x78)

my.df6 <- lapply(my.df6, function(x)
{dimnames(x) <- list(c("S","P","J","T","E"),
                     c("S","P","J","T","E"))
  x})

my.df6

A3 <- create.mxMatrix(c(0, 0, 0, 0,0, 0, 0, "0.2*J2P", "0.2*T2P", "0.2*E2P", "0.2*S2J", 0, 0, 0, 0, "0.2*S2T", 0, 0, 0, 0, 0, 0, "0.2*J2E", "0.2*T2E", 0), type="Full", nrow=5, ncol=5, byrow=TRUE, name="A3")
A3


S3 <- create.mxMatrix(c(1, 0, "0.1*ErrVarP", 0 , 0, "0.1*ErrVarJ", 0, 0, "0.1*CorJT", "0.1*ErrVarT", 0, 0, 0,0, "0.1*ErrVarE"), type="Symm", byrow=TRUE, name="S3")
S3

my.n6 <- c(181,354,349,447,178,197,305,151,222,253,182,189,185,486,161,798,78,550,224,274,514,214,345,914,221,250,211,173,100,407,32,32,153,212,224,322,356,206,368,160,167,546,844,2111,532,188,157,176,239,275,80,1126,338,268,253,108,130,217,168,180,230,430,207,989,56,193,43,247,39,160,148,216,729,563,446,621,232,128)
my.n6

S3_high <- create.mxMatrix(c(1, 0, "0.1*ErrVarP_H", 0 , 0, "0.1*ErrVarJ_H", 0, 0, "0.1*CorJT_H", "0.1*ErrVarT_H", 0, 0, 0,0, "0.1*ErrVarE_H"), type="Symm", byrow=TRUE, name="S3_High")
S3_high

source("BB9.R")
source("BBA.R")

# moderator analysis 1


my.tc <- c(0,1,0,1,0,1,1,0,1,1,1,1,0,0,0,1,0,0,1,0,1,1,1,1,0,0,0,0,1,0,0,0,1,1,0,0,0,0,0,0,1,0,1,1,1,0,1,0,0,0,0,1,0,0,0,1,0,1,0,1,0,1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,0)
my.df6_tc_low <- my.df6[which(my.tc==0)]
my.df6_tc_high <- my.df6[which(my.tc==1)]
my.n6_tc_low <- my.n6[which(my.tc==0)]
my.n6_tc_high <- my.n6[which(my.tc==1)]
stage1_tc_low.fit <- tssem1(my.df = my.df6_tc_low, n = my.n6_tc_low , method = "REM", RE.type = "Diag")
stage1_tc_high.fit <- tssem1(my.df = my.df6_tc_high, n = my.n6_tc_high , method = "REM", RE.type = "Diag")
stage2_tc_low.fit <- tssem2(stage1_tc_low.fit, Amatrix=A3, Smatrix=S3, diag.constraints=FALSE, intervals= “z”)
stage2_tc_high.fit <- tssem2(stage1_tc_high.fit, Amatrix=A3, Smatrix=S3, diag.constraints= FALSE, intervals= “z” )
stage2_tc_low <- tssem2(stage1_tc_low.fit, Amatrix=A3, Smatrix=S3, diag.constraints=FALSE, intervals=“z”, run=FALSE, model.name="low")
stage2_tc_high <- tssem2(stage1_tc_high.fit, Amatrix=A3, Smatrix=S3_high, diag.constraints=FALSE, intervals=“z”, run=FALSE, model.name="high")
stage2_constrained <- mxModel(model="same_regression_coef", stage2_tc_low, stage2_tc_high, mxFitFunctionMultigroup(c("low", "high")))
Stage2_constrained.fit <- mxRun(stage2_constrained, intervals=TRUE)
submodels.fit <- list(stage2_tc_low.fit,stage2_tc_high.fit)
subgroup.summary(submodels.fit,Stage2_constrained.fit)



# Moderator Analysis 2

my.fam <- c(0,1,1,1,0,1,1,0,1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,1,NA,0,0,0,0,0,NA,NA,1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,1,1,0,1,0,0,0,1,0,1,1,1,1,0,1,0,0,0,NA,1,0,0,1,0,1,1,0,1,1,0)
my.df6_fam_low <- my.df6[which(my.fam==0)]
my.df6_fam_high <- my.df6[which(my.fam==1)]
my.n6_fam_low <- my.n6[which(my.fam==0)]
my.n6_fam_high <- my.n6[which(my.fam==1)]
stage1_fam_low.fit <- tssem1(my.df = my.df6_fam_low, n = my.n6_fam_low , method = "REM", RE.type = "Diag”, acov = "unweighted")
stage1_fam_high.fit <- tssem1(my.df = my.df6_fam_high, n = my.n6_fam_high , method = "REM", RE.type = "Diag", acov = "unweighted")

stage1_fam_high.fit <- rerun(stage1_fam_high.fit)


stage2_fam_low.fit <- tssem2(stage1_fam_low.fit, Amatrix=A3, Smatrix=S3, diag.constraints=FALSE, intervals= “z”)
stage2_fam_high.fit <- tssem2(stage1_fam_high.fit, Amatrix=A3, Smatrix=S3, diag.constraints= FALSE, intervals= “z” )
stage2_fam_low <- tssem2(stage1_fam_low.fit, Amatrix=A3, Smatrix=S3, diag.constraints=FALSE, intervals=“z”, run=FALSE, model.name="low")
stage2_fam_high <- tssem2(stage1_fam_high.fit, Amatrix=A3, Smatrix=S3_high, diag.constraints=FALSE, intervals=“z”, run=FALSE, model.name="high")
stage2_constrained <- mxModel(model="same_regression_coef", stage2_fam_low, stage2_fam_high, mxFitFunctionMultigroup(c("low", "high")))
Stage2_constrained.fit <- mxRun(stage2_constrained, intervals=TRUE)
submodels.fit <- list(stage2_fam_low.fit,stage2_fam_high.fit)
subgroup.summary(submodels.fit,Stage2_constrained.fit)