# final transformed scog data: spasd_spins_yj_z_df
data2 <- merge(spasd_spins_demo_df, spasd_spins_yj_z_df[, c(1, 4:10)], by = 'record_id') %>%
mutate(sex = ifelse(sex == "Female", 0, 1)) %>% # female = 0, male = 1
mutate(age = (age - mean(age))/sd(age)) # Standardize age
data <- data2 [, c ( 1:8, 10, 9, 15, 13, 11, 12, 14)]
# write.csv(data,file='spasd_spins_mnlfa_input.csv', row.names=FALSE)
# Tutorial from:
### Laura Kolbe
### Last updated: 22 February 2022
### MNLFA in OpenMx with mokken::DS14 data.
### Supplementary material
## -------------------------------
## Step 1: Install and load OpenMx
## -------------------------------
library(OpenMx)
## OpenMx may run faster if it is compiled to take advantage of multiple cores.
##
## Attaching package: 'OpenMx'
## The following object is masked from 'package:psych':
##
## tr
## -----------------------------
## Step 2: Load and prepare data
## -----------------------------
library(mokken)
## Loading required package: poLCA
## Loading required package: scatterplot3d
##
## Attaching package: 'mokken'
## The following object is masked from 'package:car':
##
## recode
## The following object is masked from 'package:psych':
##
## ICC
## The following object is masked from 'package:dplyr':
##
## recode
## Create data object
## mxData() function constructs an object with additional
## information allowing it to be processed in the mxModel() function
Simulation <- c("rmet_total", "er40_total", "tasit3_lies", "EA_Score")
Mentalizing <- c("tasit2_ssar", "tasit2_psar", "tasit3_sar")
latents <- c("sim", "men")
manifests <- c(Simulation, Mentalizing)
SimMenModel <- mxModel(model="Simulation_Mentalizing_Model", type="RAM",
manifestVars = manifests,
latentVars = latents,
# factor loadings from latents to manifests
mxPath(from = "sim", to = Simulation, free = c(TRUE, TRUE, TRUE, TRUE)), # Paths from 'sim' to 'Simulation'
mxPath(from = "men", to = Mentalizing, free = c(TRUE, TRUE, TRUE)), # Paths from 'men' to 'Mentalizing'
# Allow latent variables to covary
mxPath(from="sim" , to="men", arrows=2, free=TRUE),
# Allow latent variables to have variance
mxPath(from=latents, arrows=2, free=FALSE, values=1),
# Manifest have residual variance
mxPath(from=manifests, arrows=2),
# the data to be analysed
mxData(cov(data[,manifests], use="complete.obs"), type = "cov", numObs = 585))
mxdata1 <- mxData(observed=data, type="raw")
## Define the latent factors and manifest variables
# latents <- c("sim", "men")
manVars <- c(Simulation, Mentalizing)
nv <- length(manVars)
## Step 3: Test full scalar invariance
## ------------------------------------
## Specify matrices for configural model
# Define matrices for configural model
# create three MxMatrix objects for the indicator intercepts of the configural mode
# baselines intercepts
matT0 <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE, values = 1, name = "matT0", dimnames = list(NULL,
manVars))
# effect of background variables (B1 & B2) on intercepts
matB1 <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE, values = 0, name = "matB1", dimnames = list(NULL,
manVars))
matB2 <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE, values = 0, name = "matB2", dimnames = list(NULL,
manVars))
# create three MxMatrix objects for the factor loadings of the configural mode
# baseline factor loadings
matL0 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = c(rep(c(TRUE, FALSE), 4), rep(c(FALSE, TRUE),
3)),
values = c(rep(c(1, 0), 4), rep(c(0, 1), 3)), byrow = TRUE, name = "matL0", dimnames =
list(manVars, c("Factor1", "Factor2")))
# effect of background variables (C1 & C2) on loadings
matC1 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = c(rep(c(TRUE, FALSE), 4), rep(c(FALSE, TRUE),
3)), values = 0,
byrow=TRUE, name = "matC1", dimnames = list(manVars, c("Factor1", "Factor2")))
matC2 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = c(rep(c(TRUE, FALSE), 4), rep(c(FALSE, TRUE),
3)), values = 0,
byrow=TRUE, name = "matC2", dimnames = list(manVars, c("Factor1", "Factor2")))
# create three MxMatrix objects for the residual variance of the indicators in the configural mode
# diagonal matrix containing the baseline residual variances h0
matE0 <- mxMatrix(type = "Diag", nrow = nv, ncol = nv, free = TRUE, values = 1, name = "matE0", dimnames = list(manVars,
manVars))
# diagonal matrix containing the effects of background variables (D1 & D2) on the residual variances
matD1 <- mxMatrix(type = "Diag", nrow = nv, ncol = nv, free = TRUE, values = 0, name = "matD1", dimnames = list(manVars,
manVars))
matD2 <- mxMatrix(type = "Diag", nrow = nv, ncol = nv, free = TRUE, values = 0, name = "matD2", dimnames = list(manVars,
manVars))
# matrices specified for the common-factor variances and correlation
# symmetric matrix containing the baseline common-factor variances and the baseline correlation between the two common factors
matP0 <- mxMatrix(type = "Symm", nrow = 2, ncol = 2, free = c(FALSE, TRUE, TRUE, FALSE), values = c(1, 0, 0, 1), name =
"matP0", dimnames = list(latents, latents))
# direct effects of background variables (H1 & H2) on the common-factor variances and correlation
matH1 <- mxMatrix(type = "Symm", nrow = 2, ncol = 2, free = c(FALSE, TRUE, TRUE, FALSE), values = 0, name = "matH1", dimnames
= list(latents, latents))
matH2 <- mxMatrix(type = "Symm", nrow = 2, ncol = 2, free = c(FALSE, TRUE, TRUE, FALSE), values = 0, name = "matH2", dimnames
= list(latents, latents))
matA0 <- mxMatrix(type = "Full", nrow = 2, ncol = 1, free = FALSE, values = 0, name = "matA0")
matG1 <- mxMatrix(type = "Full", nrow = 2, ncol = 1, free = FALSE, values = 0, name = "matG1")
matG2 <- mxMatrix(type = "Full", nrow = 2, ncol = 1, free = FALSE, values = 0, name = "matG2")
# Moderators
matV1 <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = FALSE, labels = "data.sex", name = "Male")
matV2 <- mxMatrix(type = "Full", nrow = 1, ncol = 1, free = FALSE, labels = "data.age", name = "Age")
# Identity and off-diagonal matrices
matIa <- mxMatrix(type = "Diag", nrow = 2, ncol = 2, free = FALSE, values = 1, name = "matIa")
matIb <- mxMatrix(type = "Full", nrow = 2, ncol = 2, free = FALSE, values = c(0, 1, 1, 0), name = "matIb")
## Specify algebra for the dependent parameters
matT <- mxAlgebra(expression = matT0 + matB1 * Male + matB2 * Age, name = "matT")
matL <- mxAlgebra(expression = matL0 + matC1 * Male + matC2 * Age, name = "matL")
matE <- mxAlgebra(expression = matE0 * exp(matD1 * Male + matD2 * Age), name = "matE")
matA <- mxAlgebra(expression = matA0 + matG1 * Male + matG2 * Age, name = "matA")
## Specify algebra for covariance matrix of factors (transformed to ensure positive definite matrices)
matVar <- mxAlgebra(expression = matP0 * exp(matH1 * Male + matH2 * Age), name = "matVar")
matR <- mxAlgebra(expression = (exp(2 * (matP0 + matH1 * Male + matH2 * Age)) - 1) /
(exp(2 * (matP0 + matH1 * Male + matH2 * Age)) + 1),
name = "matR")
matCov <- mxAlgebra(expression = (matIa * sqrt(matVar)) %*% matR %*% (matIa * sqrt(matVar)), name = "matCov")
matP <- mxAlgebra(expression = matIa * matVar + matIb * matCov, name = "matP")
## Specify model-implied matrices
matC <- mxAlgebra(expression = matL %*% matP %*% t(matL) + matE, name = "matC", dimnames = list(manVars, manVars))
matM <- mxAlgebra(expression = matT + t(matL %*% matA), name = "matM", dimnames = list(NULL, manVars))
## Specify expectation and fit function
expF <- mxExpectationNormal(covariance = "matC", means = "matM", dimnames = manVars)
fitF <- mxFitFunctionML()
## Make mxModel object and run the model
modConfig <- mxModel(model = "Configural",
matT, matT0, matB1, matB2,
matL, matL0, matC1, matC2,
matE, matE0, matD1, matD2,
matP, matP0, matH1, matH2,
matA, matA0, matG1, matG2,
matIa, matIb, matV1, matV2,
matVar, matR, matCov, matM, matC,
expF, fitF, mxdata1)
fitConfig <- mxRun(modConfig)
## Running Configural with 66 parameters
summary(fitConfig)
## Summary of Configural
##
## free parameters:
## name matrix row col Estimate Std.Error
## 1 Configural.matT0[1,1] matT0 1 rmet_total 0.098359404 0.06454101
## 2 Configural.matT0[1,2] matT0 1 er40_total 0.021759316 0.06427288
## 3 Configural.matT0[1,3] matT0 1 tasit3_lies 0.121467583 0.06447539
## 4 Configural.matT0[1,4] matT0 1 EA_Score -0.009188188 0.06400731
## 5 Configural.matT0[1,5] matT0 1 tasit2_ssar 0.029994365 0.06212116
## 6 Configural.matT0[1,6] matT0 1 tasit2_psar 0.060177254 0.06317331
## 7 Configural.matT0[1,7] matT0 1 tasit3_sar 0.024567724 0.06138927
## 8 Configural.matB1[1,1] matB1 1 rmet_total -0.181538933 0.08296560
## 9 Configural.matB1[1,2] matB1 1 er40_total -0.057351457 0.08312885
## 10 Configural.matB1[1,3] matB1 1 tasit3_lies -0.207716255 0.08316157
## 11 Configural.matB1[1,4] matB1 1 EA_Score -0.006495984 0.08265192
## 12 Configural.matB1[1,5] matB1 1 tasit2_ssar -0.073275944 0.08156382
## 13 Configural.matB1[1,6] matB1 1 tasit2_psar -0.135900532 0.08274477
## 14 Configural.matB1[1,7] matB1 1 tasit3_sar -0.064660898 0.08110740
## 15 Configural.matB2[1,1] matB2 1 rmet_total -0.154031364 0.04084488
## 16 Configural.matB2[1,2] matB2 1 er40_total -0.176964303 0.04185281
## 17 Configural.matB2[1,3] matB2 1 tasit3_lies -0.034026544 0.04095666
## 18 Configural.matB2[1,4] matB2 1 EA_Score -0.265357164 0.03946412
## 19 Configural.matB2[1,5] matB2 1 tasit2_ssar -0.162979377 0.04304467
## 20 Configural.matB2[1,6] matB2 1 tasit2_psar -0.158747763 0.04226170
## 21 Configural.matB2[1,7] matB2 1 tasit3_sar -0.185504114 0.04339885
## 22 Configural.matL0[1,1] matL0 rmet_total Factor1 0.753154051 0.07255808
## 23 Configural.matL0[2,1] matL0 er40_total Factor1 0.443358139 0.07188052
## 24 Configural.matL0[3,1] matL0 tasit3_lies Factor1 0.457178769 0.07204461
## 25 Configural.matL0[4,1] matL0 EA_Score Factor1 0.304757822 0.07368539
## 26 Configural.matL0[5,2] matL0 tasit2_ssar Factor2 0.650122741 0.06350297
## 27 Configural.matL0[6,2] matL0 tasit2_psar Factor2 0.698841564 0.06312853
## 28 Configural.matL0[7,2] matL0 tasit3_sar Factor2 0.678660975 0.06095835
## 29 Configural.matC1[1,1] matC1 rmet_total Factor1 0.028756673 0.08984743
## 30 Configural.matC1[2,1] matC1 er40_total Factor1 0.060184523 0.09209467
## 31 Configural.matC1[3,1] matC1 tasit3_lies Factor1 0.045983558 0.09170438
## 32 Configural.matC1[4,1] matC1 EA_Score Factor1 0.167043305 0.09376471
## 33 Configural.matC1[5,2] matC1 tasit2_ssar Factor2 0.151334903 0.07962086
## 34 Configural.matC1[6,2] matC1 tasit2_psar Factor2 0.122931114 0.07979084
## 35 Configural.matC1[7,2] matC1 tasit3_sar Factor2 0.093307267 0.07825309
## 36 Configural.matC2[1,1] matC2 rmet_total Factor1 -0.002560447 0.04461248
## 37 Configural.matC2[2,1] matC2 er40_total Factor1 0.062842396 0.04654638
## 38 Configural.matC2[3,1] matC2 tasit3_lies Factor1 0.020540389 0.04555188
## 39 Configural.matC2[4,1] matC2 EA_Score Factor1 0.055685086 0.04391836
## 40 Configural.matC2[5,2] matC2 tasit2_ssar Factor2 0.035331683 0.04112133
## 41 Configural.matC2[6,2] matC2 tasit2_psar Factor2 0.030787836 0.04023883
## 42 Configural.matC2[7,2] matC2 tasit3_sar Factor2 0.044731904 0.04236563
## 43 Configural.matE0[1,1] matE0 rmet_total rmet_total 0.414581968 0.08059981
## 44 Configural.matE0[2,2] matE0 er40_total er40_total 0.757843549 0.07792913
## 45 Configural.matE0[3,3] matE0 tasit3_lies tasit3_lies 0.776147709 0.07927273
## 46 Configural.matE0[4,4] matE0 EA_Score EA_Score 0.816200984 0.08124228
## 47 Configural.matE0[5,5] matE0 tasit2_ssar tasit2_ssar 0.486076582 0.06069528
## 48 Configural.matE0[6,6] matE0 tasit2_psar tasit2_psar 0.452206287 0.06135000
## 49 Configural.matE0[7,7] matE0 tasit3_sar tasit3_sar 0.430420829 0.05817839
## 50 Configural.matD1[1,1] matD1 rmet_total rmet_total -0.135209020 0.24651765
## 51 Configural.matD1[2,2] matD1 er40_total er40_total -0.046164489 0.13333475
## 52 Configural.matD1[3,3] matD1 tasit3_lies tasit3_lies -0.051968805 0.13306936
## 53 Configural.matD1[4,4] matD1 EA_Score EA_Score -0.136296015 0.13309246
## 54 Configural.matD1[5,5] matD1 tasit2_ssar tasit2_ssar -0.268905758 0.16708600
## 55 Configural.matD1[6,6] matD1 tasit2_psar tasit2_psar -0.270780179 0.18039935
## 56 Configural.matD1[7,7] matD1 tasit3_sar tasit3_sar -0.027649279 0.17009166
## 57 Configural.matD2[1,1] matD2 rmet_total rmet_total 0.018841101 0.10972931
## 58 Configural.matD2[2,2] matD2 er40_total er40_total -0.031621950 0.06653249
## 59 Configural.matD2[3,3] matD2 tasit3_lies tasit3_lies -0.029630409 0.06317208
## 60 Configural.matD2[4,4] matD2 EA_Score EA_Score -0.199370040 0.06941286
## 61 Configural.matD2[5,5] matD2 tasit2_ssar tasit2_ssar 0.176336196 0.07907991
## 62 Configural.matD2[6,6] matD2 tasit2_psar tasit2_psar 0.055324972 0.08943823
## 63 Configural.matD2[7,7] matD2 tasit3_sar tasit3_sar 0.182748257 0.07800743
## 64 Configural.matP0[1,2] matP0 sim men 1.144577211 0.20679668
## 65 Configural.matH1[1,2] matH1 sim men -0.170096707 0.21376073
## 66 Configural.matH2[1,2] matH2 sim men 0.266213811 0.12669037
## A
## 1
## 2
## 3
## 4
## 5
## 6
## 7
## 8
## 9
## 10
## 11
## 12
## 13
## 14
## 15
## 16
## 17
## 18
## 19
## 20
## 21
## 22
## 23
## 24
## 25
## 26
## 27
## 28
## 29
## 30
## 31
## 32
## 33
## 34
## 35
## 36
## 37
## 38
## 39
## 40
## 41
## 42
## 43
## 44
## 45
## 46
## 47
## 48
## 49
## 50
## 51
## 52
## 53
## 54
## 55
## 56
## 57
## 58
## 59
## 60
## 61
## 62
## 63
## 64
## 65
## 66
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 66 4029 10414.88
## Saturated: 35 4060 NA
## Independence: 14 4081 NA
## Number of observations/statistics: 600/4095
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: 2356.877 10546.88 10563.47
## BIC: -15358.352 10837.07 10627.54
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2025-01-21 23:45:07
## Wall clock time: 27.22809 secs
## optimizer: SLSQP
## OpenMx version number: 2.21.12
## Need help? See help(mxSummary)
# ## Specify matrices scalar model
matB1 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=FALSE,
values=0,
name="matB1")
matB2 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=FALSE,
values=0,
name="matB2")
matC1 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=FALSE,
values=0,
name="matC1")
matC2 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=FALSE,
values=0,
name="matC2")
matH1 <- mxMatrix(type="Symm", nrow=2, ncol=2,
free=TRUE,
values=0,
name="matH1")
matH2 <- mxMatrix(type="Symm", nrow=2, ncol=2,
free=TRUE,
values=0,
name="matH2")
matG1 <- mxMatrix(type="Full", nrow=2, ncol=1,
free=TRUE,
values=0,
name="matG1")
matG2 <- mxMatrix(type="Full", nrow=2, ncol=1,
free=TRUE,
values=0,
name="matG2")
## Make mxModel object and run the model
modScalar <- mxModel(model="Scalar",
matT, matT0, matB1, matB2,
matL, matL0, matC1, matC2,
matE, matE0, matD1, matD2,
matP, matP0, matH1, matH2,
matA, matA0, matG1, matG2,
matIa, matIb, matV1, matV2,
matVar, matR, matCov, matM, matC,
expF, fitF, mxdata1)
fitScalar <- mxRun(modScalar)
## Running Scalar with 46 parameters
summary(fitScalar)
## Summary of Scalar
##
## free parameters:
## name matrix row col Estimate Std.Error A
## 1 Scalar.matT0[1,1] matT0 1 rmet_total 0.088035925 0.05901983
## 2 Scalar.matT0[1,2] matT0 1 er40_total 0.054932681 0.05027069
## 3 Scalar.matT0[1,3] matT0 1 tasit3_lies 0.059415737 0.04920725
## 4 Scalar.matT0[1,4] matT0 1 EA_Score 0.025737244 0.05061847
## 5 Scalar.matT0[1,5] matT0 1 tasit2_ssar 0.040393231 0.05480323
## 6 Scalar.matT0[1,6] matT0 1 tasit2_psar 0.031132167 0.05606365
## 7 Scalar.matT0[1,7] matT0 1 tasit3_sar 0.042936771 0.05418376
## 8 Scalar.matL0[1,1] matL0 rmet_total Factor1 0.688389946 0.06273100
## 9 Scalar.matL0[2,1] matL0 er40_total Factor1 0.466509789 0.05038877
## 10 Scalar.matL0[3,1] matL0 tasit3_lies Factor1 0.431310346 0.04897601
## 11 Scalar.matL0[4,1] matL0 EA_Score Factor1 0.435231546 0.04882463
## 12 Scalar.matL0[5,2] matL0 tasit2_ssar Factor2 0.667894018 0.04797422
## 13 Scalar.matL0[6,2] matL0 tasit2_psar Factor2 0.690795982 0.04877925
## 14 Scalar.matL0[7,2] matL0 tasit3_sar Factor2 0.662164001 0.04805374
## 15 Scalar.matE0[1,1] matE0 rmet_total rmet_total 0.470268748 0.07215299
## 16 Scalar.matE0[2,2] matE0 er40_total er40_total 0.765150177 0.07873167
## 17 Scalar.matE0[3,3] matE0 tasit3_lies tasit3_lies 0.801173085 0.08011676
## 18 Scalar.matE0[4,4] matE0 EA_Score EA_Score 0.804079681 0.08170481
## 19 Scalar.matE0[5,5] matE0 tasit2_ssar tasit2_ssar 0.475732347 0.05741240
## 20 Scalar.matE0[6,6] matE0 tasit2_psar tasit2_psar 0.457261643 0.05658962
## 21 Scalar.matE0[7,7] matE0 tasit3_sar tasit3_sar 0.441272902 0.05494850
## 22 Scalar.matD1[1,1] matD1 rmet_total rmet_total -0.207496023 0.19053644 !
## 23 Scalar.matD1[2,2] matD1 er40_total er40_total -0.069839079 0.13157717
## 24 Scalar.matD1[3,3] matD1 tasit3_lies tasit3_lies -0.073948914 0.12895582
## 25 Scalar.matD1[4,4] matD1 EA_Score EA_Score -0.067857420 0.13259967
## 26 Scalar.matD1[5,5] matD1 tasit2_ssar tasit2_ssar -0.234764783 0.15737727
## 27 Scalar.matD1[6,6] matD1 tasit2_psar tasit2_psar -0.272791147 0.16199294
## 28 Scalar.matD1[7,7] matD1 tasit3_sar tasit3_sar -0.070258606 0.15570749
## 29 Scalar.matD2[1,1] matD2 rmet_total rmet_total -0.039268189 0.09115807
## 30 Scalar.matD2[2,2] matD2 er40_total er40_total 0.002550887 0.06320388
## 31 Scalar.matD2[3,3] matD2 tasit3_lies tasit3_lies -0.027567673 0.06103391
## 32 Scalar.matD2[4,4] matD2 EA_Score EA_Score -0.163780395 0.06958225
## 33 Scalar.matD2[5,5] matD2 tasit2_ssar tasit2_ssar 0.177363361 0.07526218
## 34 Scalar.matD2[6,6] matD2 tasit2_psar tasit2_psar 0.046028306 0.08131896
## 35 Scalar.matD2[7,7] matD2 tasit3_sar tasit3_sar 0.188040454 0.07199445
## 36 Scalar.matP0[1,2] matP0 sim men 1.171196790 0.22196606
## 37 Scalar.matH1[1,1] matH1 1 1 0.221532086 0.19455749
## 38 Scalar.matH1[1,2] matH1 1 2 -0.185406400 0.22592839
## 39 Scalar.matH1[2,2] matH1 2 2 0.336968276 0.15567521
## 40 Scalar.matH2[1,1] matH2 1 1 0.064231725 0.08917110
## 41 Scalar.matH2[1,2] matH2 1 2 0.271920075 0.13463036
## 42 Scalar.matH2[2,2] matH2 2 2 0.095482186 0.07572514
## 43 Scalar.matG1[1,1] matG1 1 1 -0.239764267 0.10856037
## 44 Scalar.matG1[2,1] matG1 2 1 -0.137602534 0.10180589
## 45 Scalar.matG2[1,1] matG2 1 1 -0.290536657 0.06027569
## 46 Scalar.matG2[2,1] matG2 2 1 -0.260288946 0.05626054
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 46 4049 10454.09
## Saturated: 35 4060 NA
## Independence: 14 4081 NA
## Number of observations/statistics: 600/4095
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: 2356.086 10546.09 10553.90
## BIC: -15447.083 10748.34 10602.31
## CFI: NA
## TLI: 1 (also known as NNFI)
## RMSEA: 0 [95% CI (NA, NA)]
## Prob(RMSEA <= 0.05): NA
## To get additional fit indices, see help(mxRefModels)
## timestamp: 2025-01-21 23:45:21
## Wall clock time: 14.05357 secs
## optimizer: SLSQP
## OpenMx version number: 2.21.12
## Need help? See help(mxSummary)
## Compare fit of unconstrained model with constrained model
miTest <- mxCompare(fitConfig, fitScalar)
miTest[2,c(7,8,9)]
## diffLL diffdf p
## 2 39.20853 20 0.006280323
miTest$p[2] < 0.001 # FALSE - H0 full scalar-and-metric invariance can not be rejected
## [1] FALSE
# ## ---------------------------------
# ## Step 4: Select anchor indicators
# ## ---------------------------------
#
# ## Run unconstrained model for each indicator
fitAbo <- list()
for (i in 1:nv){
freeparT <- matrix(FALSE, nrow=1, ncol=nv)
freeparT[i] <- TRUE
freeparL <- matrix(FALSE, nrow=nv, ncol=2)
freeparL[i,ifelse(i<=4,1,2)] <- TRUE
matB1 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=freeparT,
values=0,
name="matB1")
matB2 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=freeparT,
values=0,
name="matB2")
matC1 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=freeparL,
values=0,
byrow=TRUE,
name="matC1")
matC2 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=freeparL,
values=0,
byrow=TRUE,
name="matC2")
modAbo <- mxModel(model=paste0("All_but_", i),
matT, matT0, matB1, matB2,
matL, matL0, matC1, matC2,
matE, matE0, matD1, matD2,
matP, matP0, matH1, matH2,
matA, matA0, matG1, matG2,
matIa, matIb, matV1, matV2,
matVar, matR, matCov, matM, matC,
expF, fitF, mxdata1)
fitAbo[[i]] <- mxRun(modAbo)
}
## Running All_but_1 with 50 parameters
## Running All_but_2 with 50 parameters
## Running All_but_3 with 50 parameters
## Running All_but_4 with 50 parameters
## Running All_but_5 with 50 parameters
## Running All_but_6 with 50 parameters
## Running All_but_7 with 50 parameters
## Compare constrained model with all unconstrained models
anchorTest <- mxCompare(fitAbo, fitScalar)
anchorOut <- data.frame(Name=paste0("Indicator",1:nv),
X2=anchorTest$diffLL[seq(2,14,2)],
df=anchorTest$diffdf[seq(2,14,2)],
p=anchorTest$p[seq(2,14,2)])
anchorOut
## Name X2 df p
## 1 Indicator1 8.1290692 4 8.696301e-02
## 2 Indicator2 3.8568814 4 4.257212e-01
## 3 Indicator3 9.9052158 4 4.205494e-02
## 4 Indicator4 28.1551375 4 1.160168e-05
## 5 Indicator5 0.6969654 4 9.517026e-01
## 6 Indicator6 1.4697660 4 8.319832e-01
## 7 Indicator7 1.8663086 4 7.603311e-01
## Select two indicators per factor with smallest X2 as anchor
anchorOut[order(anchorOut$X2[1:4]),]
## Name X2 df p
## 2 Indicator2 3.856881 4 4.257212e-01
## 1 Indicator1 8.129069 4 8.696301e-02
## 3 Indicator3 9.905216 4 4.205494e-02
## 4 Indicator4 28.155137 4 1.160168e-05
anchorOut[4+order(anchorOut$X2[5:7]),]
## Name X2 df p
## 5 Indicator5 0.6969654 4 0.9517026
## 6 Indicator6 1.4697660 4 0.8319832
## 7 Indicator7 1.8663086 4 0.7603311
## Save anchors in object
anchors1 <- c(1,2)
anchors2 <- c(5,6)
## -------------------------------
## Step 5: Test partial invariance
## -------------------------------
## Specify which DIF effects can be estimated
freeparT <- matrix(data=TRUE, nrow=1, ncol=7)
freeparT[1,c(anchors1,anchors2)] <- FALSE
freeparL <- matrix(c(rep(c(TRUE,FALSE),4),rep(c(FALSE,TRUE),3)),nv,2, byrow=TRUE)
freeparL[anchors1,1] <- FALSE
freeparL[anchors2,2] <- FALSE
## Specify matrices for unconstrained model
matB1 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=freeparT,
values=0,
name="matB1")
matB2 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=freeparT,
values=0,
name="matB2")
matC1 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=freeparL,
values=0,
name="matC1")
matC2 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=freeparL,
values=0,
name="matC2")
## Make mxModel object and run the model
modAnchors <- mxModel(model="AnchorsOnly",
matT, matT0, matB1, matB2,
matL, matL0, matC1, matC2,
matE, matE0, matD1, matD2,
matP, matP0, matH1, matH2,
matA, matA0, matG1, matG2,
matIa, matIb, matV1, matV2,
matVar, matR, matCov, matM, matC,
expF, fitF, mxdata1)
fitAnchors <- mxRun(modAnchors)
## Running AnchorsOnly with 58 parameters
## Run constrained model for each indicator (except the anchors)
testIn <- c(1:7)[-c(anchors1,anchors2)]
fitApo <- list()
for (i in testIn){
freeparTa <- freeparT
freeparLa <- freeparL
freeparTa[1,i] <- FALSE
freeparLa[i,ifelse(i<4,1,2)] <- FALSE
matB1 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=freeparTa,
values=0,
name="matB1")
matB2 <- mxMatrix(type="Full", nrow=1, ncol=nv,
free=freeparTa,
values=0,
name="matB2")
matC1 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=freeparLa,
values=0,
name="matC1")
matC2 <- mxMatrix(type="Full", nrow=nv, ncol=2,
free=freeparLa,
values=0,
name="matC2")
modApo <- mxModel(model=paste0("Anchors_plus_", i),
matT, matT0, matB1, matB2,
matL, matL0, matC1, matC2,
matE, matE0, matD1, matD2,
matP, matP0, matH1, matH2,
matA, matA0, matG1, matG2,
matIa, matIb, matV1, matV2,
matVar, matR, matCov, matM, matC,
expF, fitF, mxdata1)
fitApo[[i]] <- mxRun(modApo)
}
## Running Anchors_plus_3 with 54 parameters
## Running Anchors_plus_4 with 56 parameters
## Running Anchors_plus_7 with 54 parameters
## Compare fit of unconstrained model with all constrained models
piTest <- mxCompare(fitAnchors, fitApo)
piOut <- data.frame(Name=paste0("Indicator",testIn),
X2=piTest$diffLL[2:4],
df=piTest$diffdf[2:4],
p=piTest$p[2:4],
p.bonferroni=p.adjust(p=piTest$p[2:4], method="bonferroni"),
p.BH=p.adjust(p=piTest$p[2:4], method="BH"))
piOut
## Name X2 df p p.bonferroni p.BH
## 1 Indicator3 6.942214 4 0.138973642 0.41692093 0.20846046
## 2 Indicator4 11.014348 2 0.004057557 0.01217267 0.01217267
## 3 Indicator7 1.780703 4 0.776010949 1.00000000 0.77601095
# ## -----------------------
# ## Final model with OpenMx
# ## -----------------------
#
# ## Save start time
# start_OpenMx <- Sys.time()
#
# ## Specify all non-DIF indicators
# finalIn <- c(1,2,3,4,5,6,8,9,10,11,12,14)
#
# ## Specify free parameters
# freeparTb <- freeparT
# freeparLb <- freeparL
#
# for(i in finalIn){
# freeparTb[1,i] <- FALSE
# freeparLb[i,ifelse(i<8,1,2)] <- FALSE
# }
#
# ## Re-specify matrices
# matB1 <- mxMatrix(type="Full", nrow=1, ncol=nv,
# free=freeparTb,
# values=0,
# name="matB1")
# matB2 <- mxMatrix(type="Full", nrow=1, ncol=nv,
# free=freeparTb,
# values=0,
# name="matB2")
# matC1 <- mxMatrix(type="Full", nrow=nv, ncol=2,
# free=freeparLb,
# values=0,
# name="matC1")
# matC2 <- mxMatrix(type="Full", nrow=nv, ncol=2,
# free=freeparLb,
# values=0,
# name="matC2")
#
# ## Fit final model
# modOpenmxPartial <- mxModel(model="PartialInvariance",
# matT, matT0, matB1, matB2,
# matL, matL0, matC1, matC2,
# matE, matE0, matD1, matD2,
# matP, matP0, matH1, matH2,
# matA, matA0, matG1, matG2,
# matIa, matIb, matV1, matV2,
# matVar, matR, matCov, matM, matC,
# expF, fitF, mxdata1)
# fitOpenmxPartial <- mxRun(modOpenmxPartial)
#
# ## Save end time
# end_OpenMx <- Sys.time()
#
#
# ## ---------------------------
# ## Create plots of final model
# ## ---------------------------
#
# ## Functions for expected X7 and X13 scores
# FexpX7 <- function(theta, Male, Age){
# fitOpenmxPartial$matT0$values[,7] + fitOpenmxPartial$matB1$values[,7]*Male +
# fitOpenmxPartial$matB2$values[,7]*Age + (fitOpenmxPartial$matL0$values[7,1] +
# fitOpenmxPartial$matC1$values[7,1]*Male + fitOpenmxPartial$matC2$values[7,1]*
# Age)*theta
# }
#
# FexpX13 <- function(theta, Male, Age){
# fitOpenmxPartial$matT0$values[,13] + fitOpenmxPartial$matB1$values[,13]*Male +
# fitOpenmxPartial$matB2$values[,13]*Age + (fitOpenmxPartial$matL0$values[13,2] +
# fitOpenmxPartial$matC1$values[13,2]*Male + fitOpenmxPartial$matC2$values[13,2]*
# Age)*theta
# }
#
# ## Plot tracelines for DIF indicators
# library(ggplot2)
# myTheme <- theme(axis.title = element_text(size = 12),
# axis.text.y = element_text(size = 10),
# axis.text.x = element_text(size = 10),
# strip.text = element_text(size = 10),
# legend.text = element_text(size = 12),
# legend.title = element_text(size = 12),
# legend.key.size = grid::unit(5, units = "char"),
# legend.key = element_rect(fill = "white", colour = "white"),
# panel.background = element_rect(fill = "white", color = "grey"),
# panel.grid.minor = element_line(colour = "white"))
#
# All <- expand.grid("Theta" = rep(seq(-4, 4, 0.01), 6), "Male" = 0:1, "Age" = c(-2,2),
# "expX7" = NA, "expX13" = NA)
# for (i in 1:nrow(All)){All$expX7[i] <- FexpX7(All$Theta[i], All$Male[i], All$Age[i])
# All$expX13[i] <- FexpX13(All$Theta[i], All$Male[i], All$Age[i])
# }
# All$Male <- factor(All$Male, labels = c("Female","Male"))
# All$Age <- factor(All$Age, labels = c("Age = 38", "Age = 80"))
# # Age is converted from standardized score to actual age score:
# # Mean age = 58.66913, SD age = 10.54776
# # -2 --> -2*10.54776+58.66913 = 37.57361
# # 2 --> 2*10.54776+58.66913 = 79.76465
#
# ggplot(All, aes(x = Theta, y = expX7, lty = Age, colour = Male)) +
# geom_line(size = 1) +
# scale_color_manual(values=c("grey80", "black")) +
# scale_linetype_manual(values=c(1, 5)) +
# labs(x = "SI", y = "Expected X7") +
# myTheme
#
# ggplot(All, aes(x = Theta, y = expX13, lty = Age, colour = Male)) +
# geom_line(size = 1) +
# scale_color_manual(values=c("grey80", "black")) +
# scale_linetype_manual(values=c(1, 5)) +
# labs(x = "NA", y = "Expected X13") +
# myTheme
#
# ## Functions for moderated factor parameters
# FacMean1 <- function(Male, Age){
# fitOpenmxPartial$matA0$values[1,1] + fitOpenmxPartial$matG1$values[1,1]*Male +
# fitOpenmxPartial$matG2$values[1,1]*Age
# }
#
# FacMean2 <- function(Male, Age){
# fitOpenmxPartial$matA0$values[2,1] + fitOpenmxPartial$matG1$values[2,1]*Male +
# fitOpenmxPartial$matG2$values[2,1]*Age
# }
#
# FacVar1 <- function(Male, Age){
# fitOpenmxPartial$matP0$values[1,1] * exp(fitOpenmxPartial$matH1$values[1,1]*Male +
# fitOpenmxPartial$matH2$values[1,1]*Age)
# }
#
# FacVar2 <- function(Male, Age){
# fitOpenmxPartial$matP0$values[2,2] * exp(fitOpenmxPartial$matH1$values[2,2]*Male +
# fitOpenmxPartial$matH2$values[2,2]*Age)
# }
#
# FacCorr <- function(Male, Age){
# (exp(2*(fitOpenmxPartial$matP0$values[1,2] + fitOpenmxPartial$matH1$values[1,2]*Male +
# fitOpenmxPartial$matH2$values[1,2]*Age)) - 1)/(exp(2*(fitOpenmxPartial$matP0$values[1,2] +
# fitOpenmxPartial$matH1$values[1,2]*Male + fitOpenmxPartial$matH2$values[1,2]*Age)) + 1)
# }
#
# ## Plots factor parameters
# All2 <- expand.grid("Male" = 0:1, "Age" = seq(-3, 3, 0.01), "A_1" = NA,
# "A_2" = NA, "P_1" = NA, "P_2" = NA, "Corr" = NA)
# for (i in 1:nrow(All2)){
# All2$A_1[i] <- FacMean1(All2$Male[i], All2$Age[i])
# All2$A_2[i] <- FacMean2(All2$Male[i], All2$Age[i])
# All2$P_1[i] <- FacVar1(All2$Male[i], All2$Age[i])
# All2$P_2[i] <- FacVar2(All2$Male[i], All2$Age[i])
# All2$Corr[i] <- FacCorr(All2$Male[i], All2$Age[i])
# }
# All2$Male <- factor(All2$Male, labels = c("Female","Male"))
#
# ggplot(All2, aes(x = Age, y = A_1, colour = Male)) +
# geom_line(size = 1) +
# scale_color_manual(values=c("grey80", "black")) +
# scale_linetype_manual(values=c(1, 5)) +
# labs(x = "Age", y = "SI mean") +
# scale_x_continuous(breaks=c(-2, 0, 2), labels=c("38", "59","80")) +
# myTheme
#
# ggplot(All2, aes(x = Age, y = A_2, colour = Male)) +
# geom_line(size = 1) +
# scale_color_manual(values=c("grey80", "black")) +
# scale_linetype_manual(values=c(1, 5)) +
# labs(x = "Age", y = "NA mean") +
# scale_x_continuous(breaks=c(-2, 0, 2), labels=c("38", "59","80")) +
# myTheme
#
# ggplot(All2, aes(x = Age, y = P_1, colour = Male)) +
# geom_line(size = 1) +
# scale_color_manual(values=c("grey80", "black")) +
# scale_linetype_manual(values=c(1, 5)) +
# labs(x = "Age", y = "SI variance") +
# scale_x_continuous(breaks=c(-2, 0, 2), labels=c("38", "59","80")) +
# myTheme
#
# ggplot(All2, aes(x = Age, y = P_2, colour = Male)) +
# geom_line(size = 1) +
# scale_color_manual(values=c("grey80", "black")) +
# scale_linetype_manual(values=c(1, 5)) +
# labs(x = "Age", y = "NA variance") +
# scale_x_continuous(breaks=c(-2, 0, 2), labels=c("38", "59","80")) +
# myTheme
#
# ggplot(All2, aes(x = Age, y = Corr, colour = Male)) +
# geom_line(size = 1) +
# scale_color_manual(values=c("grey80", "black")) +
# scale_linetype_manual(values=c(1, 5)) +
# labs(x = "Age", y = "Factor correlation") +
# scale_x_continuous(breaks=c(-2, 0, 2), labels=c("38", "59","80")) +
# myTheme
#
#
# ## ------------------------
# ## Compare results to Mplus
# ## ------------------------
#
# library(MplusAutomation)
#
# ## Save start time
# start_Mplus <- Sys.time()
#
# ## Create folder
# pathfix <- "~/FinalModel"
# dir.create(pathfix)
#
# ## Store data in folder
# prepareMplusData(df=DS14, filename=paste0(pathfix, "/DS14dat.dat"))
#
# ## Create and run partial-invariance model
# modMplusPartial <- '
# DATA:
# FILE = "DS14dat.dat";
#
# VARIABLE:
# NAMES = Male Age x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14;
# CONSTRAINT = Male Age;
# MISSING = .;
#
# ANALYSIS:
# estimator is ML;
# TOLERANCE = 1;
#
# MODEL:
# SI by x1-x7* (l1-l7);
# NA by x8-x14* (l8-l14);
# SI (SIvar);
# NA (NAvar);
# SI with NA (cov);
# [SI] (SImean);
# [NA] (NAmean);
# [x1-x14] (t1-t14);
# x1-x14 (e1-e14);
# Male @1;
# Age @1;
# [Male @0];
# [Age @0];
#
# MODEL CONSTRAINT:
# NEW (
# h1_SI h1_NA g1_SI g1_NA
# h2_SI h2_NA g2_SI g2_NA
# p0 h1_cov h2_cov
# t7_0 b7_1 b7_2 t13_0 b13_1 b13_2
# l7_0 c7_1 c7_2 l13_0 c13_1 c13_2
# e1_0 e2_0 e3_0 e4_0 e5_0 e6_0 e7_0 e8_0 e9_0 e10_0 e11_0 e12_0 e13_0 e14_0
# d1_1 d2_1 d3_1 d4_1 d5_1 d6_1 d7_1 d8_1 d9_1 d10_1 d11_1 d12_1 d13_1 d14_1
# d1_2 d2_2 d3_2 d4_2 d5_2 d6_2 d7_2 d8_2 d9_2 d10_2 d11_2 d12_2 d13_2 d14_2);
#
# SIvar = 1 * EXP(h1_SI*Male + h2_SI*Age);
# NAvar = 1 * EXP(h1_NA*Male + h2_NA*Age);
# SImean = 0 + g1_SI*Male + g2_SI*Age;
# NAmean = 0 + g1_NA*Male + g2_NA*Age;
#
# cov = SQRT(EXP(h1_SI*Male + h2_SI*Age))*
# SQRT(EXP(h1_NA*Male + h2_NA*Age))*
# (EXP(2*(p0 + h1_cov*Male + h2_cov*Age))-1)/
# (EXP(2*(p0 + h1_cov*Male + h2_cov*Age))+1);
#
# e1 = e1_0 * EXP(d1_1*Male + d1_2*Age);
# e2 = e2_0 * EXP(d2_1*Male + d2_2*Age);
# e3 = e3_0 * EXP(d3_1*Male + d3_2*Age);
# e4 = e4_0 * EXP(d4_1*Male + d4_2*Age);
# e5 = e5_0 * EXP(d5_1*Male + d5_2*Age);
# e6 = e6_0 * EXP(d6_1*Male + d6_2*Age);
# e7 = e7_0 * EXP(d7_1*Male + d7_2*Age);
# e8 = e8_0 * EXP(d8_1*Male + d8_2*Age);
# e9 = e9_0 * EXP(d9_1*Male + d9_2*Age);
# e10 = e10_0 * EXP(d10_1*Male + d10_2*Age);
# e11 = e11_0 * EXP(d11_1*Male + d11_2*Age);
# e12 = e12_0 * EXP(d12_1*Male + d12_2*Age);
# e13 = e13_0 * EXP(d13_1*Male + d13_2*Age);
# e14 = e14_0 * EXP(d14_1*Male + d14_2*Age);
# t7 = t7_0 + b7_1*Male + b7_2*Age;
# l7 = l7_0 + c7_1*Male + c7_2*Age;
# t13 = t13_0 + b13_1*Male + b13_2*Age;
# l13 = l13_0 + c13_1*Male + c13_2*Age;'
#
# cat(modMplusPartial, file = paste0(pathfix, "/modPartial.inp", sep = ""))
#
# ## Run model and save output
# runModels(pathfix)
# fitMplusPartial <- readModels(pathfix)
#
# ## Save end time
# end_Mplus <- Sys.time()
#
#
# ## Compare parameters estimates
# summary(fitOpenmxPartial)$parameters
# fitMplusPartial$parameters
# # NOTE: in the Mplus parameter list, the parameters being a function of
# # the background variables are also listed (can be ignored) and equal 999.000
# # because they differ per individual
#
# ## Compare times
# end_OpenMx - start_OpenMx
# end_Mplus - start_Mplus
#
#