# final transformed scog data: spasd_spins_yj_z_df
data <- 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
# 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", "mean_ea")
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, # list the measured variables (boxes)
latentVars = latents, # list the latent variables (circles)
# factor loadings from latents to manifests
mxPath(from="sim", to=Simulation),# factor loadings
mxPath(from="men", to=Mentalizing), # factor loadings
# 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.0),
# 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))
mxdata <- 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
matT0 <- mxMatrix(type = "Full", nrow = 1, ncol = nv, free = TRUE, values = 1, name = "matT0", dimnames = list(NULL, manVars))
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))
matL0 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = c(rep(c(TRUE, FALSE), length(Simulation)),
rep(c(FALSE, TRUE), length(Mentalizing))),
values = c(rep(c(1, 0), length(Simulation)),
rep(c(0, 1), length(Mentalizing))), byrow = TRUE,
name = "matL0", dimnames = list(manVars, latents))
matC1 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = TRUE, values = 0, name = "matC1", dimnames = list(manVars, latents))
matC2 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = TRUE, values = 0, name = "matC2", dimnames = list(manVars, latents))
matE0 <- mxMatrix(type = "Diag", nrow = nv, ncol = nv, free = TRUE, values = 1, name = "matE0", dimnames = list(manVars, manVars))
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))
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))
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, mxdata)
fitConfig <- mxRun(modConfig)
## Running Configural with 80 parameters
## Warning: In model 'Configural' Optimizer returned a non-zero status code 5. The
## Hessian at the solution does not appear to be convex. See
## ?mxCheckIdentification for possible diagnosis (Mx status RED).
summary(fitConfig)
## Summary of Configural
##
## The Hessian at the solution does not appear to be convex. See ?mxCheckIdentification for possible diagnosis (Mx status RED).
##
## free parameters:
## name matrix row col Estimate
## 1 Configural.matT0[1,1] matT0 1 rmet_total 0.09477599
## 2 Configural.matT0[1,2] matT0 1 er40_total 0.02488555
## 3 Configural.matT0[1,3] matT0 1 tasit3_lies 0.12017705
## 4 Configural.matT0[1,4] matT0 1 mean_ea 0.02284209
## 5 Configural.matT0[1,5] matT0 1 tasit2_ssar 0.05302106
## 6 Configural.matT0[1,6] matT0 1 tasit2_psar 0.07332273
## 7 Configural.matT0[1,7] matT0 1 tasit3_sar 0.04860970
## 8 Configural.matB1[1,1] matB1 1 rmet_total -0.17821003
## 9 Configural.matB1[1,2] matB1 1 er40_total -0.06037405
## 10 Configural.matB1[1,3] matB1 1 tasit3_lies -0.19851614
## 11 Configural.matB1[1,4] matB1 1 mean_ea -0.10333209
## 12 Configural.matB1[1,5] matB1 1 tasit2_ssar -0.10812469
## 13 Configural.matB1[1,6] matB1 1 tasit2_psar -0.16071379
## 14 Configural.matB1[1,7] matB1 1 tasit3_sar -0.09617108
## 15 Configural.matB2[1,1] matB2 1 rmet_total -0.14373863
## 16 Configural.matB2[1,2] matB2 1 er40_total -0.18582660
## 17 Configural.matB2[1,3] matB2 1 tasit3_lies -0.01304927
## 18 Configural.matB2[1,4] matB2 1 mean_ea -0.27472782
## 19 Configural.matB2[1,5] matB2 1 tasit2_ssar -0.17008524
## 20 Configural.matB2[1,6] matB2 1 tasit2_psar -0.13906830
## 21 Configural.matB2[1,7] matB2 1 tasit3_sar -0.19848482
## 22 Configural.matL0[1,1] matL0 rmet_total sim 0.72780296
## 23 Configural.matL0[2,1] matL0 er40_total sim 0.49774635
## 24 Configural.matL0[3,1] matL0 tasit3_lies sim 0.43454322
## 25 Configural.matL0[4,1] matL0 mean_ea sim 0.28643262
## 26 Configural.matL0[5,2] matL0 tasit2_ssar men 0.61897698
## 27 Configural.matL0[6,2] matL0 tasit2_psar men 0.63136627
## 28 Configural.matL0[7,2] matL0 tasit3_sar men 0.82840583
## 29 Configural.matC1[1,1] matC1 rmet_total sim 385.96291915
## 30 Configural.matC1[2,1] matC1 er40_total sim 337.45570642
## 31 Configural.matC1[3,1] matC1 tasit3_lies sim 273.96538053
## 32 Configural.matC1[4,1] matC1 mean_ea sim 305.74395951
## 33 Configural.matC1[5,1] matC1 tasit2_ssar sim 311.05368732
## 34 Configural.matC1[6,1] matC1 tasit2_psar sim 319.65003276
## 35 Configural.matC1[7,1] matC1 tasit3_sar sim 286.48080500
## 36 Configural.matC1[1,2] matC1 rmet_total men -386.45790538
## 37 Configural.matC1[2,2] matC1 er40_total men -338.06202164
## 38 Configural.matC1[3,2] matC1 tasit3_lies men -274.40309080
## 39 Configural.matC1[4,2] matC1 mean_ea men -305.93523638
## 40 Configural.matC1[5,2] matC1 tasit2_ssar men -311.06441612
## 41 Configural.matC1[6,2] matC1 tasit2_psar men -319.65949846
## 42 Configural.matC1[7,2] matC1 tasit3_sar men -286.64597420
## 43 Configural.matC2[1,1] matC2 rmet_total sim 0.29318374
## 44 Configural.matC2[2,1] matC2 er40_total sim 0.10597091
## 45 Configural.matC2[3,1] matC2 tasit3_lies sim 0.25412162
## 46 Configural.matC2[4,1] matC2 mean_ea sim 0.12578445
## 47 Configural.matC2[5,1] matC2 tasit2_ssar sim 0.27688028
## 48 Configural.matC2[6,1] matC2 tasit2_psar sim 0.28879002
## 49 Configural.matC2[7,1] matC2 tasit3_sar sim 0.04278989
## 50 Configural.matC2[1,2] matC2 rmet_total men -0.26937598
## 51 Configural.matC2[2,2] matC2 er40_total men 0.03934818
## 52 Configural.matC2[3,2] matC2 tasit3_lies men -0.23091715
## 53 Configural.matC2[4,2] matC2 mean_ea men -0.06474535
## 54 Configural.matC2[5,2] matC2 tasit2_ssar men -0.21323824
## 55 Configural.matC2[6,2] matC2 tasit2_psar men -0.33727265
## 56 Configural.matC2[7,2] matC2 tasit3_sar men 0.13290205
## 57 Configural.matE0[1,1] matE0 rmet_total rmet_total 0.29951682
## 58 Configural.matE0[2,2] matE0 er40_total er40_total 0.69630791
## 59 Configural.matE0[3,3] matE0 tasit3_lies tasit3_lies 0.71470365
## 60 Configural.matE0[4,4] matE0 mean_ea mean_ea 0.67530856
## 61 Configural.matE0[5,5] matE0 tasit2_ssar tasit2_ssar 0.39492901
## 62 Configural.matE0[6,6] matE0 tasit2_psar tasit2_psar 0.46138792
## 63 Configural.matE0[7,7] matE0 tasit3_sar tasit3_sar 0.16368265
## 64 Configural.matD1[1,1] matD1 rmet_total rmet_total 0.38605071
## 65 Configural.matD1[2,2] matD1 er40_total er40_total -0.27992178
## 66 Configural.matD1[3,3] matD1 tasit3_lies tasit3_lies 0.03328612
## 67 Configural.matD1[4,4] matD1 mean_ea mean_ea -0.02996164
## 68 Configural.matD1[5,5] matD1 tasit2_ssar tasit2_ssar -0.22068136
## 69 Configural.matD1[6,6] matD1 tasit2_psar tasit2_psar -0.45349645
## 70 Configural.matD1[7,7] matD1 tasit3_sar tasit3_sar 0.45112842
## 71 Configural.matD2[1,1] matD2 rmet_total rmet_total -0.16231154
## 72 Configural.matD2[2,2] matD2 er40_total er40_total 0.09426167
## 73 Configural.matD2[3,3] matD2 tasit3_lies tasit3_lies -0.04917150
## 74 Configural.matD2[4,4] matD2 mean_ea mean_ea 0.20026231
## 75 Configural.matD2[5,5] matD2 tasit2_ssar tasit2_ssar 0.34746054
## 76 Configural.matD2[6,6] matD2 tasit2_psar tasit2_psar 0.30111120
## 77 Configural.matD2[7,7] matD2 tasit3_sar tasit3_sar -0.40575354
## 78 Configural.matP0[1,2] matP0 sim men 0.83085189
## 79 Configural.matH1[1,2] matH1 sim men 6.13436037
## 80 Configural.matH2[1,2] matH2 sim men -0.08269571
## Std.Error A
## 1 0.06203192 !
## 2 0.06492392
## 3 0.06376565
## 4 0.06166015 !
## 5 0.05910186 !
## 6 0.06356951
## 7 0.06053335 !
## 8 0.08212223 !
## 9 0.08310480 !
## 10 0.08318752
## 11 0.08317433 !
## 12 0.07877222 !
## 13 0.08378139
## 14 0.08021623 !
## 15 0.04264368 !
## 16 0.04374573
## 17 0.04213831 !
## 18 0.04577094
## 19 0.04370902
## 20 0.04325578 !
## 21 0.04402799
## 22 0.07491487 !
## 23 0.07508034 !
## 24 0.07296849 !
## 25 0.07449528 !
## 26 0.05858897 !
## 27 0.06309764
## 28 0.05324090 !
## 29 NA
## 30 NA
## 31 NA
## 32 NA
## 33 NA
## 34 NA
## 35 NA
## 36 NA
## 37 NA
## 38 NA
## 39 NA
## 40 NA
## 41 NA
## 42 NA
## 43 0.11166227 !
## 44 0.12169964
## 45 0.11198070
## 46 0.11485422
## 47 0.10621428
## 48 0.11579596 !
## 49 0.11289907
## 50 0.11649210 !
## 51 0.12053248
## 52 0.11279099 !
## 53 0.11204068
## 54 0.10428145 !
## 55 0.11208298 !
## 56 0.10164879
## 57 0.07286295 !
## 58 0.07625900 !
## 59 0.07420133 !
## 60 0.07096974 !
## 61 0.04737960
## 62 0.05980383 !
## 63 0.04810020
## 64 0.25933458 !
## 65 0.18065174 !
## 66 0.13827578 !
## 67 0.14592227 !
## 68 0.15932120
## 69 0.17387349 !
## 70 0.24037155
## 71 0.11337879 !
## 72 0.07282754 !
## 73 0.06424468 !
## 74 0.07011473
## 75 0.07350468
## 76 0.07960522 !
## 77 0.27429510 !
## 78 0.11857259 !
## 79 NA !
## 80 0.05317646 !
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 80 3874 9709.586
## Saturated: 35 3919 NA
## Independence: 14 3940 NA
## Number of observations/statistics: 585/3954
##
##
## ** Information matrix is not positive definite (not at a candidate optimum).
## Be suspicious of these results. At minimum, do not trust the standard errors.
##
## Information Criteria:
## | df Penalty | Parameters Penalty | Sample-Size Adjusted
## AIC: 1961.586 9869.586 9895.301
## BIC: -14974.038 10219.315 9965.344
## 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: 2024-10-02 16:25:11
## Wall clock time: 42.58822 secs
## optimizer: SLSQP
## OpenMx version number: 2.21.12
## Need help? See help(mxSummary)
### GOT TO HERE ###
# ## 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)
# summary(fitScalar)
#
# ## Compare fit of unconstrained model with constrained model
# miTest <- mxCompare(fitConfig, fitScalar)
# miTest[2,c(7,8,9)]
# miTest$p[2] < 0.001
#
#
# ## ---------------------------------
# ## Step 4: Select anchor indicators
# ## ---------------------------------
#
# ## Run unconstrained model for each indicator
# fitAbo <- list()
#
# for (i in 1:14){
# freeparT <- matrix(FALSE, nrow=1, ncol=14)
# freeparT[i] <- TRUE
# freeparL <- matrix(FALSE, nrow=14, ncol=2)
# freeparL[i,ifelse(i<8,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)
# }
#
# ## Compare constrained model with all unconstrained models
# anchorTest <- mxCompare(fitAbo, fitScalar)
# anchorOut <- data.frame(Name=paste0("Indicator",1:14),
# X2=anchorTest$diffLL[seq(2,28,2)],
# df=anchorTest$diffdf[seq(2,28,2)],
# p=anchorTest$p[seq(2,28,2)])
# anchorOut
#
# ## Select two indicators per factor with smallest X2 as anchor
# anchorOut[order(anchorOut$X2[1:7]),]
# anchorOut[7+order(anchorOut$X2[8:14]),]
#
# ## Save anchors in object
# anchors1 <- c(5,6)
# anchors2 <- c(8,9)
#
#
# ## -------------------------------
# ## Step 5: Test partial invariance
# ## -------------------------------
#
# ## Specify which DIF effects can be estimated
# freeparT <- matrix(data=TRUE, nrow=1, ncol=14)
# freeparT[1,c(anchors1,anchors2)] <- FALSE
#
# freeparL <- matrix(c(rep(c(TRUE,FALSE),7),rep(c(FALSE,TRUE),7)),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)
#
# ## Run constrained model for each indicator (except the anchors)
# testIn <- c(1:14)[-c(anchors1,anchors2)]
# fitApo <- list()
#
# for (i in testIn){
# freeparTa <- freeparT
# freeparLa <- freeparL
# freeparTa[1,i] <- FALSE
# freeparLa[i,ifelse(i<8,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)
# }
#
# ## Compare fit of unconstrained model with all constrained models
# piTest <- mxCompare(fitAnchors, fitApo)
# piOut <- data.frame(Name=paste0("Indicator",testIn),
# X2=piTest$diffLL[2:11],
# df=piTest$diffdf[2:11],
# p=piTest$p[2:11],
# p.bonferroni=p.adjust(p=piTest$p[2:11], method="bonferroni"),
# p.BH=p.adjust(p=piTest$p[2:11], method="BH"))
# piOut
#
#
# ## -----------------------
# ## 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
#
#