# 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)
## Warning: package 'OpenMx' was built under R version 4.3.3
##
## Attaching package: 'OpenMx'
## The following object is masked from 'package:psych':
##
## tr
## -----------------------------
## Step 2: Load and prepare data
## -----------------------------
library(mokken)
## Warning: package 'mokken' was built under R version 4.3.3
## Loading required package: poLCA
## Warning: package 'poLCA' was built under R version 4.3.3
## 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
mxdata <- mxData(observed=data, type="raw")
## Indicate names and number of indicators
manVars <- colnames(data)[9:15]
nv <- length(manVars) # Number of manifest variables
## ------------------------------------
## 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), 7), rep(c(FALSE, TRUE), 7)),
values = c(rep(c(1, 0), 7), rep(c(0, 1), 7)), byrow = TRUE, name = "matL0", dimnames = list(manVars, c("Factor1", "Factor2")))
## Warning in matrix(values, nrow, ncol, byrow = byrow): data length differs from
## size of matrix: [28 != 7 x 2]
## Warning in matrix(free, nrow, ncol, byrow = byrow): data length differs from
## size of matrix: [28 != 7 x 2]
matC1 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = TRUE, values = 0, name = "matC1", dimnames = list(manVars, c("Factor1", "Factor2")))
matC2 <- mxMatrix(type = "Full", nrow = nv, ncol = 2, free = TRUE, values = 0, name = "matC2", dimnames = list(manVars, c("Factor1", "Factor2")))
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")
matH1 <- mxMatrix(type = "Symm", nrow = 2, ncol = 2, free = c(FALSE, TRUE, TRUE, FALSE), values = 0, name = "matH1")
matH2 <- mxMatrix(type = "Symm", nrow = 2, ncol = 2, free = c(FALSE, TRUE, TRUE, FALSE), values = 0, name = "matH2")
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 Std.Error
## 1 Configural.matT0[1,1] matT0 1 rmet_total 1.079486e-01 64.16183
## 2 Configural.matT0[1,2] matT0 1 er40_total 3.459066e-02 55.02358
## 3 Configural.matT0[1,3] matT0 1 tasit2_ssar 3.846244e-02 73.90512
## 4 Configural.matT0[1,4] matT0 1 tasit2_psar 5.991764e-02 71.60465
## 5 Configural.matT0[1,5] matT0 1 tasit3_lies 1.323200e-01 54.18878
## 6 Configural.matT0[1,6] matT0 1 tasit3_sar 2.772923e-02 73.21621
## 7 Configural.matT0[1,7] matT0 1 mean_ea 2.213360e-02 53.04572
## 8 Configural.matB1[1,1] matB1 1 rmet_total -1.962331e-01 50.60431
## 9 Configural.matB1[1,2] matB1 1 er40_total -7.278194e-02 43.24425
## 10 Configural.matB1[1,3] matB1 1 tasit2_ssar -1.036236e-01 59.04315
## 11 Configural.matB1[1,4] matB1 1 tasit2_psar -1.533808e-01 58.23348
## 12 Configural.matB1[1,5] matB1 1 tasit3_lies -2.076722e-01 42.18739
## 13 Configural.matB1[1,6] matB1 1 tasit3_sar -9.156552e-02 56.92911
## 14 Configural.matB1[1,7] matB1 1 mean_ea -1.081664e-01 40.52460
## 15 Configural.matB2[1,1] matB2 1 rmet_total -1.589890e-01 64.41853
## 16 Configural.matB2[1,2] matB2 1 er40_total -2.045204e-01 52.82008
## 17 Configural.matB2[1,3] matB2 1 tasit2_ssar -1.812061e-01 66.71674
## 18 Configural.matB2[1,4] matB2 1 tasit2_psar -1.470093e-01 66.29069
## 19 Configural.matB2[1,5] matB2 1 tasit3_lies -2.356401e-02 52.88743
## 20 Configural.matB2[1,6] matB2 1 tasit3_sar -2.102784e-01 66.23852
## 21 Configural.matB2[1,7] matB2 1 mean_ea -2.921706e-01 49.16317
## 22 Configural.matL0[1,1] matL0 rmet_total Factor1 5.120849e-01 57.86400
## 23 Configural.matL0[2,1] matL0 er40_total Factor1 3.688710e-01 48.82776
## 24 Configural.matL0[3,1] matL0 tasit2_ssar Factor1 6.488319e-01 67.15229
## 25 Configural.matL0[4,1] matL0 tasit2_psar Factor1 6.682211e-01 64.60179
## 26 Configural.matL0[5,1] matL0 tasit3_lies Factor1 2.940117e-01 49.30575
## 27 Configural.matL0[6,1] matL0 tasit3_sar Factor1 8.087997e-01 68.94065
## 28 Configural.matL0[7,1] matL0 mean_ea Factor1 2.659317e-01 46.35587
## 29 Configural.matC1[1,1] matC1 rmet_total Factor1 4.017986e-01 46.10074
## 30 Configural.matC1[2,1] matC1 er40_total Factor1 1.881635e-01 38.40574
## 31 Configural.matC1[3,1] matC1 tasit2_ssar Factor1 2.877973e-01 54.49916
## 32 Configural.matC1[4,1] matC1 tasit2_psar Factor1 2.916382e-01 52.82586
## 33 Configural.matC1[5,1] matC1 tasit3_lies Factor1 2.176438e-01 38.90872
## 34 Configural.matC1[6,1] matC1 tasit3_sar Factor1 1.038803e-01 52.46455
## 35 Configural.matC1[7,1] matC1 mean_ea Factor1 3.702052e-01 35.97013
## 36 Configural.matC1[1,2] matC1 rmet_total Factor2 -5.900348e-01 37.77054
## 37 Configural.matC1[2,2] matC1 er40_total Factor2 -4.688142e-01 31.46777
## 38 Configural.matC1[3,2] matC1 tasit2_ssar Factor2 -1.127780e-01 42.46118
## 39 Configural.matC1[4,2] matC1 tasit2_psar Factor2 -1.414132e-01 42.29799
## 40 Configural.matC1[5,2] matC1 tasit3_lies Factor2 -3.576397e-01 30.32104
## 41 Configural.matC1[6,2] matC1 tasit3_sar Factor2 -1.321486e-01 41.38810
## 42 Configural.matC1[7,2] matC1 mean_ea Factor2 -4.050947e-01 28.04381
## 43 Configural.matC2[1,1] matC2 rmet_total Factor1 -8.752343e-02 56.96404
## 44 Configural.matC2[2,1] matC2 er40_total Factor1 6.026424e-02 46.53067
## 45 Configural.matC2[3,1] matC2 tasit2_ssar Factor1 2.074382e-02 58.73894
## 46 Configural.matC2[4,1] matC2 tasit2_psar Factor1 -7.713099e-02 59.23647
## 47 Configural.matC2[5,1] matC2 tasit3_lies Factor1 -9.805239e-02 47.85003
## 48 Configural.matC2[6,1] matC2 tasit3_sar Factor1 7.712301e-02 63.15947
## 49 Configural.matC2[7,1] matC2 mean_ea Factor1 4.043018e-02 43.34004
## 50 Configural.matC2[1,2] matC2 rmet_total Factor2 4.038428e-01 46.81585
## 51 Configural.matC2[2,2] matC2 er40_total Factor2 3.503076e-01 37.25430
## 52 Configural.matC2[3,2] matC2 tasit2_ssar Factor2 2.662220e-01 46.52048
## 53 Configural.matC2[4,2] matC2 tasit2_psar Factor2 2.916632e-01 47.17285
## 54 Configural.matC2[5,2] matC2 tasit3_lies Factor2 4.823935e-01 41.18934
## 55 Configural.matC2[6,2] matC2 tasit3_sar Factor2 1.279379e-01 46.84120
## 56 Configural.matC2[7,2] matC2 mean_ea Factor2 1.895716e-01 35.16100
## 57 Configural.matE0[1,1] matE0 rmet_total rmet_total 5.076286e-01 52.74036
## 58 Configural.matE0[2,2] matE0 er40_total er40_total 7.266759e-01 41.50079
## 59 Configural.matE0[3,3] matE0 tasit2_ssar tasit2_ssar 4.010445e-01 65.87372
## 60 Configural.matE0[4,4] matE0 tasit2_psar tasit2_psar 4.462095e-01 57.89952
## 61 Configural.matE0[5,5] matE0 tasit3_lies tasit3_lies 6.815166e-01 44.00782
## 62 Configural.matE0[6,6] matE0 tasit3_sar tasit3_sar 2.499707e-01 93.26847
## 63 Configural.matE0[7,7] matE0 mean_ea mean_ea 6.806681e-01 43.68466
## 64 Configural.matD1[1,1] matD1 rmet_total rmet_total -2.400114e-01 19.28681
## 65 Configural.matD1[2,2] matD1 er40_total er40_total -1.329155e-01 22.68127
## 66 Configural.matD1[3,3] matD1 tasit2_ssar tasit2_ssar -2.768163e-01 19.33463
## 67 Configural.matD1[4,4] matD1 tasit2_psar tasit2_psar -4.087895e-01 18.77637
## 68 Configural.matD1[5,5] matD1 tasit3_lies tasit3_lies 2.295421e-02 23.33538
## 69 Configural.matD1[6,6] matD1 tasit3_sar tasit3_sar 2.805691e-01 20.04982
## 70 Configural.matD1[7,7] matD1 mean_ea mean_ea 6.470562e-04 22.46691
## 71 Configural.matD2[1,1] matD2 rmet_total rmet_total 7.605247e-03 27.66076
## 72 Configural.matD2[2,2] matD2 er40_total er40_total 8.707774e-02 30.30900
## 73 Configural.matD2[3,3] matD2 tasit2_ssar tasit2_ssar 2.959964e-01 27.41822
## 74 Configural.matD2[4,4] matD2 tasit2_psar tasit2_psar 2.170820e-01 27.54402
## 75 Configural.matD2[5,5] matD2 tasit3_lies tasit3_lies -1.112921e-01 26.74054
## 76 Configural.matD2[6,6] matD2 tasit3_sar tasit3_sar 1.413972e-01 22.02515
## 77 Configural.matD2[7,7] matD2 mean_ea mean_ea 1.884901e-01 29.02595
## 78 Configural.matP0[1,2] matP0 1 2 2.651196e-16 NA
## 79 Configural.matH1[1,2] matH1 1 2 5.188666e-01 28.91563
## 80 Configural.matH2[1,2] matH2 1 2 -2.567503e-01 37.92812
## 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 !
## 67 !
## 68 !
## 69 !
## 70 !
## 71 !
## 72 !
## 73 !
## 74 !
## 75 !
## 76 !
## 77 !
## 78 !
## 79 !
## 80 !
##
## Model Statistics:
## | Parameters | Degrees of Freedom | Fit (-2lnL units)
## Model: 80 3874 9736.655
## 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: 1988.655 9896.655 9922.369
## BIC: -14946.970 10246.384 9992.412
## 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-08-15 12:53:53
## Wall clock time: 210.1073 secs
## optimizer: SLSQP
## OpenMx version number: 2.21.11
## 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
#
#