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

source("http://www.suzannejak.nl/subgroup.functions.R")

(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,0.24,NA,NA,-0.13,1,NA,NA,-0.34,1,NA,NA,1,NA,1)))
(x3 <- vec2symMat(c(1,0.25,NA,NA,-0.14,1,NA,NA,-0.35,1,NA,NA,1,NA,1)))
(x4 <- vec2symMat(c(1,NA,0.36,0.34,-0.18,1,NA,NA,NA,1,0.04,-0.1,1,-0.32,1)))
(x5 <- 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)))
(x6 <- vec2symMat(c(1,0.7,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x7 <- vec2symMat(c(1,NA,0.611,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x8 <- vec2symMat(c(1,NA,0.43,NA,-0.29,1,NA,NA,NA,1,NA,-0.28,1,NA,1)))
(x9 <- vec2symMat(c(1,NA,NA,0.6,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x10 <- vec2symMat(c(1,0.29,0.29,NA,NA,1,0.48,NA,NA,1,NA,NA,1,NA,1)))
(x11 <- 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)))
(x12 <- vec2symMat(c(1,0.42,NA,0.49,NA,1,NA,0.74,NA,1,NA,NA,1,NA,1)))
(x13 <- vec2symMat(c(1,0.218,NA,0.437,-0.171,1,NA,0.376,-0.644,1,NA,NA,1,-0.336,1)))
(x14 <- vec2symMat(c(1,0.561,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x15 <- 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)))
(x16 <- vec2symMat(c(1,NA,NA,0.369,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x17 <- vec2symMat(c(1,0.62,NA,NA,-0.64,1,NA,NA,-0.69,1,NA,NA,1,NA,1)))
(x18 <- vec2symMat(c(1,NA,0.444,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x19 <- vec2symMat(c(1,0.18,0.66,NA,-0.18,1,0.38,NA,-0.9,1,NA,-0.43,1,NA,1)))
(x20 <- vec2symMat(c(1,0.54,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x21 <- vec2symMat(c(1,0.309,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x22 <- vec2symMat(c(1,0.31,0.66,NA,NA,1,0.29,NA,NA,1,NA,NA,1,NA,1)))
(x23 <- vec2symMat(c(1,0.234,NA,0.289,-0.238,1,NA,0.301,-0.196,1,NA,NA,1,-0.271,1)))
(x24 <- vec2symMat(c(1,0.195656566,NA,NA,-0.231157025,1,NA,NA,-0.328838384,1,NA,NA,1,NA,1)))
(x25 <- vec2symMat(c(1,0.230833333,NA,NA,-0.173461538,1,NA,NA,-0.524405594,1,NA,NA,1,NA,1)))
(x26 <- vec2symMat(c(1,NA,0.28,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x27 <- vec2symMat(c(1,0.62,0.54,0.33,NA,1,0.54,0.65,NA,1,0.34,NA,1,NA,1)))
(x28 <- vec2symMat(c(1,NA,0.454,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x29 <- vec2symMat(c(1,NA,0.832,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x30 <- vec2symMat(c(1,NA,NA,NA,-0.303,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x31 <- 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)))
(x32 <- vec2symMat(c(1,0.338,NA,0.396,-0.119,1,NA,0.53,-0.583,1,NA,NA,1,-0.52,1)))
(x33 <- vec2symMat(c(1,NA,0.655,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x34 <- vec2symMat(c(1,0.408,0.624,NA,NA,1,0.461,NA,NA,1,NA,NA,1,NA,1)))
(x35 <- vec2symMat(c(1,0.092,NA,NA,-0.234,1,NA,NA,-0.415,1,NA,NA,1,NA,1)))
(x36 <- vec2symMat(c(1,NA,NA,0.766,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x37 <- vec2symMat(c(1,NA,0.684,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x38 <- vec2symMat(c(1,NA,NA,NA,-0.282,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x39 <- vec2symMat(c(1,NA,NA,0.445,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x40 <- vec2symMat(c(1,0.613,0.51,NA,NA,1,0.739,NA,NA,1,NA,NA,1,NA,1)))
(x41 <- vec2symMat(c(1,NA,NA,0.549,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x42 <- vec2symMat(c(1,NA,NA,0.77,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x43 <- vec2symMat(c(1,0.273,0.5148,NA,NA,1,0.4016,NA,NA,1,NA,NA,1,NA,1)))
(x44 <- vec2symMat(c(1,0.304,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x45 <- vec2symMat(c(1,0.456,NA,NA,-0.31,1,NA,NA,-0.244,1,NA,NA,1,NA,1)))
(x46 <- vec2symMat(c(1,0.537,NA,NA,-0.398,1,NA,NA,-0.36,1,NA,NA,1,NA,1)))
(x47 <- 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)))
(x48 <- vec2symMat(c(1,0.47,0.48,NA,NA,1,0.41,NA,NA,1,NA,NA,1,NA,1)))
(x49 <- vec2symMat(c(1,NA,NA,0.619,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x50 <- vec2symMat(c(1,0.15,NA,0.17,NA,1,NA,0.13,NA,1,NA,NA,1,NA,1)))
(x51 <- vec2symMat(c(1,0.15,NA,0.28,NA,1,NA,0.23,NA,1,NA,NA,1,NA,1)))
(x52 <- vec2symMat(c(1,0.57,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x53 <- vec2symMat(c(1,0.4,0.33,NA,NA,1,0.53,NA,NA,1,NA,NA,1,NA,1)))
(x54 <- vec2symMat(c(1,0.386,0.63,0.384,NA,1,0.531,0.549,NA,1,0.367,NA,1,NA,1)))
(x55 <- vec2symMat(c(1,NA,NA,0.43,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x56 <- vec2symMat(c(1,0.51,0.66,NA,-0.61,1,0.76,NA,-0.57,1,NA,-0.66,1,NA,1)))
(x57 <- vec2symMat(c(1,0.499,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x58 <- vec2symMat(c(1,0.602,0.716,NA,NA,1,0.622,NA,NA,1,NA,NA,1,NA,1)))
(x59 <- vec2symMat(c(1,NA,NA,0.315,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x60 <- vec2symMat(c(1,NA,0.6414,0.0027,-0.1386,1,NA,NA,NA,1,0.0459,-0.1925,1,-0.629,1)))
(x61 <- vec2symMat(c(1,NA,NA,0.353,-0.385,1,NA,NA,NA,1,NA,NA,1,-0.478,1)))
(x62 <- vec2symMat(c(1,0.51,0.68,0.35,NA,1,0.47,0.5,NA,1,0.48,NA,1,NA,1)))
(x63 <- vec2symMat(c(1,0.62,0.77,NA,NA,1,0.63,NA,NA,1,NA,NA,1,NA,1)))
(x64 <- vec2symMat(c(1,NA,0.53,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x65 <- vec2symMat(c(1,0.28,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x66 <- vec2symMat(c(1,0.57,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x67 <- vec2symMat(c(1,0.32,0.561,NA,NA,1,0.505,NA,NA,1,NA,NA,1,NA,1)))
(x68 <- vec2symMat(c(1,0.795,0.761,NA,-0.671,1,0.719,NA,-0.688,1,NA,-0.633,1,NA,1)))
(x69 <- vec2symMat(c(1,0.493,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x70 <- vec2symMat(c(1,0.45,0.71,NA,NA,1,0.66,NA,NA,1,NA,NA,1,NA,1)))
(x71 <- vec2symMat(c(1,NA,0.546,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x72 <- vec2symMat(c(1,0.445,NA,0.368,NA,1,NA,0.301,NA,1,NA,NA,1,NA,1)))
(x73 <- vec2symMat(c(1,0.32,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x74 <- vec2symMat(c(1,0.25,NA,0.25,-0.2,1,NA,0.45,-0.37,1,NA,NA,1,-0.49,1)))
(x75 <- vec2symMat(c(1,0.529,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x76 <- vec2symMat(c(1,0.44,0.67,NA,NA,1,0.39,NA,NA,1,NA,NA,1,NA,1)))
(x77 <- vec2symMat(c(1,NA,NA,0.582,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x78 <- vec2symMat(c(1,0.365,NA,0.606,NA,1,NA,0.474,NA,1,NA,NA,1,NA,1)))
(x79 <- vec2symMat(c(1,0.41,0.351,NA,NA,1,0.347,NA,NA,1,NA,NA,1,NA,1)))
(x80 <- vec2symMat(c(1,0.4,NA,NA,-0.39,1,NA,NA,-0.65,1,NA,NA,1,NA,1)))
(x81 <- vec2symMat(c(1,0.45,NA,NA,-0.34,1,NA,NA,-0.64,1,NA,NA,1,NA,1)))
(x82 <- vec2symMat(c(1,NA,NA,0.57,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x83 <- vec2symMat(c(1,NA,NA,NA,-0.48,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x84 <- vec2symMat(c(1,0.496,NA,0.51,NA,1,NA,0.575,NA,1,NA,NA,1,NA,1)))
(x85 <- vec2symMat(c(1,0.3257,NA,NA,-0.4036,1,NA,NA,-0.6762,1,NA,NA,1,NA,1)))
(x86 <- vec2symMat(c(1,0.684,0.346,NA,-0.115,1,0.155,NA,-0.153,1,NA,-0.63,1,NA,1)))
(x87 <- vec2symMat(c(1,NA,NA,0.594,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x88 <- vec2symMat(c(1,0.588,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x89 <- vec2symMat(c(1,0.495,NA,0.52,NA,1,NA,0.45,NA,1,NA,NA,1,NA,1)))
(x90 <- vec2symMat(c(1,0.434,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x91 <- vec2symMat(c(1,0.54,0.4,NA,NA,1,0.57,NA,NA,1,NA,NA,1,NA,1)))
(x92 <- vec2symMat(c(1,NA,0.584,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x93 <- vec2symMat(c(1,NA,0.14,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x94 <- vec2symMat(c(1,0.422,NA,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x95 <- vec2symMat(c(1,0.68,NA,NA,-0.47,1,NA,NA,-0.32,1,NA,NA,1,NA,1)))
(x96 <- vec2symMat(c(1,0.62,0.68,NA,NA,1,0.618,NA,NA,1,NA,NA,1,NA,1)))
(x97 <- vec2symMat(c(1,0.52,0.59,NA,-0.38,1,0.56,NA,-0.32,1,NA,-0.38,1,NA,1)))
(x98 <- vec2symMat(c(1,0.4102,0.4796,NA,-0.3729,1,0.6113,NA,-0.642,1,NA,-0.5103,1,NA,1)))
(x99 <- vec2symMat(c(1,NA,0.213,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x100 <- vec2symMat(c(1,NA,NA,0.592,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x101 <- vec2symMat(c(1,0.612,NA,0.725,NA,1,NA,0.659,NA,1,NA,NA,1,NA,1)))
(x102 <- vec2symMat(c(1,0.591,NA,0.69,NA,1,NA,0.647,NA,1,NA,NA,1,NA,1)))
(x103 <- vec2symMat(c(1,NA,NA,0.26,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x104 <- vec2symMat(c(1,0.471,0.539,NA,NA,1,0.673,NA,NA,1,NA,NA,1,NA,1)))
(x105 <- vec2symMat(c(1,0.489,0.449,NA,NA,1,0.545,NA,NA,1,NA,NA,1,NA,1)))
(x106 <- vec2symMat(c(1,0.648,0.686,NA,NA,1,0.768,NA,NA,1,NA,NA,1,NA,1)))
(x107 <- vec2symMat(c(1,NA,0.25,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x108 <- vec2symMat(c(1,NA,0.42,NA,-0.58,1,NA,NA,NA,1,NA,-0.57,1,NA,1)))
(x109 <- vec2symMat(c(1,0.562,0.774,NA,-0.514,1,0.688,NA,-0.578,1,NA,-0.562,1,NA,1)))
(x110 <- vec2symMat(c(1,0.22,NA,0.19,-0.1,1,NA,0.22,-0.37,1,NA,NA,1,-0.2,1)))
(x111 <- vec2symMat(c(1,NA,NA,0.85,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x112 <- vec2symMat(c(1,NA,0.532,NA,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x113 <- vec2symMat(c(1,NA,NA,0.53,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x114 <- vec2symMat(c(1,0.55,0.62,NA,NA,1,0.69,NA,NA,1,NA,NA,1,NA,1)))
(x115 <- vec2symMat(c(1,0.601,0.747,NA,NA,1,0.667,NA,NA,1,NA,NA,1,NA,1)))
(x116 <- vec2symMat(c(1,NA,NA,0.598,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x117 <- vec2symMat(c(1,0.405,0.458,NA,NA,1,0.444,NA,NA,1,NA,NA,1,NA,1)))
(x118 <- vec2symMat(c(1,NA,NA,0.578,NA,1,NA,0.301,NA,1,NA,NA,1,NA,1)))
(x119 <- vec2symMat(c(1,NA,NA,0.22,NA,1,NA,NA,NA,1,NA,NA,1,NA,1)))
(x120 <- vec2symMat(c(1,NA,NA,0.22,NA,1,NA,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,x79,x80,x81,x82,x83,x84,x85,x86,x87,x88,x89,x90,x91,x92,x93,x94,x95,x96,x97,x98,x99,x100,x101,x102,x103,x104,x105,x106,x107,x108,x109,x110,x111,x112,x113,x114,x115,x116,x117,x118,x119,x120)

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

my.n6 <- c(181,349,447,178,197,305,151,222,253,182,185,486,161,798,78,550,224,274,514,345,914,221,250,211,173,100,407,32,32,153,212,224,322,356,206,368,160,167,546,844,2111,532,188,176,239,275,80,1126,338,268,253,108,217,168,180,230,430,207,989,56,193,43,247,39,160,148,216,729,621,128,385,263,154,270,156,375,540,495,372,220,296,215,223,189,229,201,260,120,408,225,728,137,152,214,73,403,298,400,473,306,388,421,201,163,163,523,371,728,403,336,262,253,242,244,262,685,204,200,797,237)

my.n6 

model1 <- 'P ~  s2p*S + t2p*T + j2p*J + e2p*E
           E ~  j2e*J + t2e*T + s2e*S
           J ~  s2j*S
           T ~  s2t*S
	   S ~~ 1*S
           J ~~ jWITHt*T
           P ~~ Errp*P
           J ~~ Errj*J
           T ~~ Errt*T 
           E ~~ Erre*E'

RAM1 <- lavaan2RAM(model1, obs.variables=c("S", "P", "J", "T","E"))
RAM1

random1 <- tssem1(my.df6, my.n6, method="REM", RE.type="Diag")
summary(random1)

random2 <- tssem2(random1, Amatrix=RAM1$A, Smatrix=RAM1$S)
summary(random2)

my.cu <- c(1,1,1,1,0,0,0,1,0,0,1,0,1,0,1,1,1,0,0,1,1,1,0,1,1,1,1,1,1,0,0,1,1,1,0,0,1,NA,0,1,1,1,1,1,1,1,0,0,0,1,1,0,0,1,1,1,NA,0,1,1,1,1,0,0,1,1,0,0,0,1,1,0,0,1,0,0,1,1,0,1,NA,1,1,0,0,1,0,0,0,0,0,0,1,0,1,0,0,0,1,0,0,0,NA,1,1,0,1,0,1,1,0,0,1,0,1,1,0,0,1,0)

my.df6_cu_low <- my.df6[which(my.cu==0)]
my.df6_cu_high <- my.df6[which(my.cu==1)]


my.n6_cu_low <- my.n6[which(my.cu==0)]
my.n6_cu_high <- my.n6[which(my.cu==1)]


stage1_g1.fit <- tssem1(Cov = my.df6_cu_low , n = my.n6_cu_low , method = "REM", RE.type = "Diag")
stage1_g2.fit <- tssem1(Cov = my.df6_cu_high , n = my.n6_cu_high, method = "REM", RE.type = "Diag")

summary(stage1_g1.fit)
summary(stage1_g2.fit)

stage2_g1.fit <- tssem2(stage1_g1.fit, Amatrix=RAM1$A, Smatrix=RAM1$S)
stage2_g2.fit <- tssem2(stage1_g2.fit, Amatrix=RAM1$A, Smatrix=RAM1$S)

summary(stage2_g1.fit)
summary(stage2_g2.fit)

## Varying parameters in S matrix
model2 <- 'P ~  s2p*S + t2p*T + j2p*J + e2p*E
           E ~  j2e*J + t2e*T + s2e*S
           J ~  s2j*S
           T ~  s2t*S
	   S ~~ 1*S
           J ~~ g2jWITHt*T
           P ~~ g2Errp*P
           J ~~ g2Errj*J
           T ~~ g2Errt*T 
           E ~~ g2Erre*E'
     
RAM2 <- lavaan2RAM(model2, obs.variables=c("S", "P", "J", "T","E"))
RAM2


# Create the models for the two groups, make sure to set the argument run=FALSE
stage2_g1 <- tssem2(stage1_g1.fit, Amatrix=RAM1$A, Smatrix=RAM1$S, run=FALSE, model.name="g1")

stage2_g2 <- tssem2(stage1_g2.fit, Amatrix=RAM1$A, Smatrix=RAM2$S, run=FALSE, model.name="g2")

# Create the multigroup model
stage2_constrained <- mxModel(model="same_regression_coef", stage2_g1, stage2_g2, mxFitFunctionMultigroup(c("g1", "g2")))

# Fit multigroup model with equality constraints
Stage2_constrained.fit <- mxRun(stage2_constrained, intervals=TRUE)

# first make a list of the fitted models in the separate groups
submodels.fit <- list(stage2_g1.fit,stage2_g2.fit)

subgroup.summary(submodels.fit,Stage2_constrained.fit, print.est=TRUE)

## Varying parameters in A matrix
model2 <- 'P ~  g2s2p*S + g2t2p*T + g2j2p*J + g2e2p*E
           E ~  g2j2e*J + g2t2e*T + g2s2e*S
           J ~  g2s2j*S
           T ~  s2t*S
	   S ~~ 1*S
           J ~~ g2jWITHt*T
           P ~~ g2Errp*P
           J ~~ g2Errj*J
           T ~~ g2Errt*T 
           E ~~ g2Erre*E'

RAM2 <- lavaan2RAM(model2, obs.variables=c("S", "P", "J", "T","E"))
RAM2

stage2_g1 <- tssem2(stage1_g1.fit, Amatrix=RAM1$A, Smatrix=RAM1$S, run=FALSE, model.name="g1")
stage2_g2 <- tssem2(stage1_g2.fit, Amatrix=RAM2$A, Smatrix=RAM1$S, run=FALSE, model.name="g2")

# Create the multigroup model
stage2_constrained <- mxModel(model="same_regression_coef", stage2_g1, stage2_g2, mxFitFunctionMultigroup(c("g1", "g2")))

# Fit multigroup model with equality constraints
Stage2_constrained.fit <- mxRun(stage2_constrained, intervals=TRUE)

# first make a list of the fitted models in the separate groups
submodels.fit <- list(stage2_g1.fit,stage2_g2.fit)

subgroup.summary(submodels.fit,Stage2_constrained.fit, print.est=TRUE)