Multi-Factor Measurement Model (OpenMx)

library(lavaan)
library(R2spa)
library(OpenMx)
library(umx)
library(mirt)

With a Linear Factor Model

The example is from https://lavaan.ugent.be/tutorial/sem.html.

# CFA
my_cfa <- "
  # latent variables
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
"
(fs_dat <- get_fs(PoliticalDemocracy, model = my_cfa, std.lv = TRUE)) |> head()
#>     fs_ind60   fs_dem60 fs_ind60_se fs_dem60_se ind60_by_fs_ind60
#> 1 -0.8101568 -1.3119114   0.1859987   0.3012901         0.9553858
#> 2  0.1888466 -1.3644831   0.1859987   0.3012901         0.9553858
#> 3  1.0960931  1.3107705   0.1859987   0.3012901         0.9553858
#> 4  1.8702043  1.4849083   0.1859987   0.3012901         0.9553858
#> 5  1.2446060  0.9193277   0.1859987   0.3012901         0.9553858
#> 6  0.3348621  0.4886331   0.1859987   0.3012901         0.9553858
#>   ind60_by_fs_dem60 dem60_by_fs_ind60 dem60_by_fs_dem60 ev_fs_ind60
#> 1        0.05816994        0.01834111         0.8688889  0.03459552
#> 2        0.05816994        0.01834111         0.8688889  0.03459552
#> 3        0.05816994        0.01834111         0.8688889  0.03459552
#> 4        0.05816994        0.01834111         0.8688889  0.03459552
#> 5        0.05816994        0.01834111         0.8688889  0.03459552
#> 6        0.05816994        0.01834111         0.8688889  0.03459552
#>   ecov_fs_dem60_fs_ind60 ev_fs_dem60
#> 1            0.004017388  0.09077571
#> 2            0.004017388  0.09077571
#> 3            0.004017388  0.09077571
#> 4            0.004017388  0.09077571
#> 5            0.004017388  0.09077571
#> 6            0.004017388  0.09077571

tspa_fit <- tspa(model = "dem60 ~ ind60", data = fs_dat,
                 fsT = attr(fs_dat, "fsT"),
                 fsL = attr(fs_dat, "fsL"))
parameterestimates(tspa_fit)
#>         lhs op      rhs   est    se     z pvalue ci.lower ci.upper
#> 1     ind60 =~ fs_ind60 0.955 0.000    NA     NA    0.955    0.955
#> 2     ind60 =~ fs_dem60 0.058 0.000    NA     NA    0.058    0.058
#> 3     dem60 =~ fs_ind60 0.018 0.000    NA     NA    0.018    0.018
#> 4     dem60 =~ fs_dem60 0.869 0.000    NA     NA    0.869    0.869
#> 5  fs_ind60 ~~ fs_ind60 0.035 0.000    NA     NA    0.035    0.035
#> 6  fs_ind60 ~~ fs_dem60 0.004 0.000    NA     NA    0.004    0.004
#> 7  fs_dem60 ~~ fs_dem60 0.091 0.000    NA     NA    0.091    0.091
#> 8     dem60  ~    ind60 0.460 0.113 4.089      0    0.240    0.681
#> 9     ind60 ~~    ind60 1.000 0.169 5.900      0    0.668    1.332
#> 10    dem60 ~~    dem60 0.788 0.150 5.267      0    0.495    1.081

Using OpenMx

Create OpenMx model without latent variables. Here we also use the umx package.

fsreg_umx <- umxLav2RAM(
    "
    dem60 ~ ind60
    dem60 + ind60 ~ 1
    ",
    printTab = FALSE)
#> 
#> ?plot.MxModel options: std, means, digits, strip_zero, file, splines=T/F/ortho,..., min=, max =, same = , fixed, resid= 'circle|line|none'

Create loading and error covariance matrices (need reordering)

# Loading
matL <- mxMatrix(
    type = "Full", nrow = 2, ncol = 2,
    free = FALSE,
    values = attr(fs_dat, "fsL")[2:1, 2:1],
    name = "L"
)
# Error
matE <- mxMatrix(
    type = "Symm", nrow = 2, ncol = 2,
    free = FALSE,
    values = attr(fs_dat, "fsT")[2:1, 2:1],
    name = "E"
)

Run 2S-PA

tspa_mx <- tspa_mx_model(fsreg_umx, data = fs_dat,
                         mat_ld = matL, mat_ev = matE,
                         fs_lv_names = c(ind60 = "fs_ind60",
                                         dem60 = "fs_dem60"))
# Run OpenMx
tspa_mx_fit <- mxRun(tspa_mx)
#> Running 2SPAD with 5 parameters
# Summarize the results (takes 20 seconds or so)
summary(tspa_mx_fit)
#> Summary of 2SPAD 
#>  
#> free parameters:
#>               name matrix   row   col      Estimate Std.Error A
#> 1   ind60_to_dem60   m1.A dem60 ind60  4.604652e-01 0.1126087  
#> 2 dem60_with_dem60   m1.S dem60 dem60  7.879711e-01 0.1495940  
#> 3 ind60_with_ind60   m1.S ind60 ind60  1.000000e+00 0.1694829  
#> 4     one_to_dem60   m1.M     1 dem60  1.083067e-07 0.1105181  
#> 5     one_to_ind60   m1.M     1 ind60 -3.096769e-08 0.1176358  
#> 
#> Model Statistics: 
#>                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
#>        Model:              5                    295              393.7496
#>    Saturated:             NA                     NA                    NA
#> Independence:             NA                     NA                    NA
#> Number of observations/statistics: 75/300
#> 
#> Information Criteria: 
#>       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
#> AIC:      -196.2504               403.7496                 404.6192
#> BIC:      -879.9094               415.3370                 399.5784
#> 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-09-06 04:19:54 
#> Wall clock time: 0.05677581 secs 
#> optimizer:  SLSQP 
#> OpenMx version number: 2.21.12 
#> Need help?  See help(mxSummary)
# Standardize coefficients
mxStandardizeRAMpaths(tspa_mx_fit, SE = TRUE)$m1[1, ]
#>        name          label matrix   row   col Raw.Value    Raw.SE Std.Value
#> 1 m1.A[1,2] ind60_to_dem60      A dem60 ind60 0.4604652 0.1126087 0.4604654
#>      Std.SE
#> 1 0.1001049

# Alternative specification
tspa_mx <- tspa_mx_model(fsreg_umx, data = fs_dat,
                         mat_ld = attr(fs_dat, "fsL"),
                         mat_ev = attr(fs_dat, "fsT"))

With Mean Structure

my_cfa <- "
  # latent variables
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    x1 + y1 ~ 0
    ind60 + dem60 ~ 1
"
fs_dat <- get_fs(PoliticalDemocracy, model = my_cfa,
                 meanstructure = TRUE)

lavaan

tspa_fit <- tspa(model = "dem60 ~ ind60\ndem60 + ind60 ~ 1", data = fs_dat,
                 fsT = attr(fs_dat, "fsT"),
                 fsL = attr(fs_dat, "fsL"),
                 fsb = attr(fs_dat, "fsb"))
parameterEstimates(tspa_fit)
#>         lhs op      rhs    est    se      z pvalue ci.lower ci.upper
#> 1     ind60 =~ fs_ind60  0.955 0.000     NA     NA    0.955    0.955
#> 2     ind60 =~ fs_dem60  0.182 0.000     NA     NA    0.182    0.182
#> 3     dem60 =~ fs_ind60  0.006 0.000     NA     NA    0.006    0.006
#> 4     dem60 =~ fs_dem60  0.869 0.000     NA     NA    0.869    0.869
#> 5  fs_ind60 ~~ fs_ind60  0.016 0.000     NA     NA    0.016    0.016
#> 6  fs_ind60 ~~ fs_dem60  0.006 0.000     NA     NA    0.006    0.006
#> 7  fs_dem60 ~~ fs_dem60  0.398 0.000     NA     NA    0.398    0.398
#> 8  fs_ind60 ~1           0.193 0.000     NA     NA    0.193    0.193
#> 9  fs_dem60 ~1          -0.203 0.000     NA     NA   -0.203   -0.203
#> 10    dem60  ~    ind60  1.439 0.352  4.089  0.000    0.749    2.129
#> 11    dem60 ~1          -1.810 1.794 -1.009  0.313   -5.327    1.706
#> 12    ind60 ~1           5.054 0.079 64.155  0.000    4.900    5.209
#> 13    ind60 ~~    ind60  0.449 0.076  5.900  0.000    0.300    0.598
#> 14    dem60 ~~    dem60  3.453 0.656  5.267  0.000    2.168    4.738
standardizedSolution(tspa_fit)
#>         lhs op      rhs est.std    se       z pvalue ci.lower ci.upper
#> 1     ind60 =~ fs_ind60   0.973 0.004 267.474   0.00    0.966    0.980
#> 2     ind60 =~ fs_dem60   0.061 0.006   9.671   0.00    0.049    0.074
#> 3     dem60 =~ fs_ind60   0.019 0.002   9.092   0.00    0.015    0.023
#> 4     dem60 =~ fs_dem60   0.918 0.011  83.338   0.00    0.897    0.940
#> 5  fs_ind60 ~~ fs_ind60   0.036 0.006   6.124   0.00    0.024    0.047
#> 6  fs_ind60 ~~ fs_dem60   0.072 0.000      NA     NA    0.072    0.072
#> 7  fs_dem60 ~~ fs_dem60   0.101 0.017   6.124   0.00    0.069    0.134
#> 8  fs_ind60 ~1            0.294 0.024  12.247   0.00    0.247    0.341
#> 9  fs_dem60 ~1           -0.102 0.008 -12.247   0.00   -0.119   -0.086
#> 10    dem60  ~    ind60   0.460 0.100   4.600   0.00    0.264    0.657
#> 11    dem60 ~1           -0.865 0.817  -1.058   0.29   -2.467    0.737
#> 12    ind60 ~1            7.547 0.650  11.606   0.00    6.272    8.821
#> 13    ind60 ~~    ind60   1.000 0.000      NA     NA    1.000    1.000
#> 14    dem60 ~~    dem60   0.788 0.092   8.547   0.00    0.607    0.969

OpenMx

# Force symmetric E
matE <- attr(fs_dat, "fsT")
matE <- (matE + t(matE)) / 2
matb <- t(as.matrix(attr(fs_dat, "fsb")))
tspa_mx <- tspa_mx_model(fsreg_umx, data = fs_dat,
                         mat_ld = attr(fs_dat, "fsL"),
                         mat_ev = matE,
                         mat_int = matb)
# Run OpenMx
tspa_mx_fit <- mxRun(tspa_mx)
#> Running 2SPAD with 5 parameters
# Summarize the results (takes 20 seconds or so)
summary(tspa_mx_fit)
#> Summary of 2SPAD 
#>  
#> free parameters:
#>               name matrix   row   col   Estimate  Std.Error A
#> 1   ind60_to_dem60   m1.A dem60 ind60  1.4393156 0.35192047  
#> 2 dem60_with_dem60   m1.S dem60 dem60  3.4532716 0.65558270  
#> 3 ind60_with_ind60   m1.S ind60 ind60  0.4485416 0.07601934  
#> 4     one_to_dem60   m1.M     1 dem60 -1.8101866 1.79371569  
#> 5     one_to_ind60   m1.M     1 ind60  5.0543839 0.07878429  
#> 
#> Model Statistics: 
#>                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
#>        Model:              5                    295              444.4392
#>    Saturated:             NA                     NA                    NA
#> Independence:             NA                     NA                    NA
#> Number of observations/statistics: 75/300
#> 
#> Information Criteria: 
#>       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
#> AIC:      -145.5608               454.4392                 455.3088
#> BIC:      -829.2198               466.0266                 450.2680
#> 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-09-06 04:19:55 
#> Wall clock time: 0.04486251 secs 
#> optimizer:  SLSQP 
#> OpenMx version number: 2.21.12 
#> Need help?  See help(mxSummary)
# Standardize coefficients
mxStandardizeRAMpaths(tspa_mx_fit, SE = TRUE)$m1[1, ]
#>        name          label matrix   row   col Raw.Value    Raw.SE Std.Value
#> 1 m1.A[1,2] ind60_to_dem60      A dem60 ind60  1.439316 0.3519205 0.4604657
#>      Std.SE
#> 1 0.1000874

Compare to joint model with OpenMx

jreg_umx_fit <- umxRAM(
    "
    # latent variables
    ind60 =~ x1 + x2 + x3
    dem60 =~ y1 + y2 + y3 + y4
    # latent regression
    dem60 ~ ind60
    ",
    data = PoliticalDemocracy)
#> 2 latent variables were created:ind60, dem60.
#> You have raw data, but no means model. I added
#> mxPath('one', to = manifestVars)
#> Running latent_variables with 22 parameters
#> ?umxSummary options: std=T|F', digits=, report= 'html', filter= 'NS' & more
#> Running Saturated latent_variables with 35 parameters
#> Running Independence latent_variables with 14 parameters
#> 
#> 
#> Table: Parameter loadings for model 'latent_variables'
#> 
#> |   |name             | Estimate|SE    |type             |
#> |:--|:----------------|--------:|:-----|:----------------|
#> |1  |ind60_to_x1      |    1.000|0     |Factor loading   |
#> |2  |ind60_to_x2      |    2.180|0.14  |Factor loading   |
#> |3  |ind60_to_x3      |    1.818|0.152 |Factor loading   |
#> |5  |dem60_to_y1      |    1.000|0     |Factor loading   |
#> |6  |dem60_to_y2      |    1.419|0.212 |Factor loading   |
#> |7  |dem60_to_y3      |    1.095|0.167 |Factor loading   |
#> |8  |dem60_to_y4      |    1.428|0.191 |Factor loading   |
#> |4  |ind60_to_dem60   |    1.439|0.377 |Factor to factor |
#> |16 |ind60_with_ind60 |    0.449|0.087 |Factor Variance  |
#> |17 |dem60_with_dem60 |    3.453|0.917 |Factor Variance  |
#> |18 |one_to_y1        |    5.465|0.301 |Mean             |
#> |19 |one_to_y2        |    4.256|0.453 |Mean             |
#> |20 |one_to_y3        |    6.563|0.376 |Mean             |
#> |21 |one_to_y4        |    4.453|0.384 |Mean             |
#> |22 |one_to_x1        |    5.054|0.084 |Mean             |
#> |23 |one_to_x2        |    4.792|0.173 |Mean             |
#> |24 |one_to_x3        |    3.558|0.161 |Mean             |
#> |9  |y1_with_y1       |    2.404|0.586 |Residual         |
#> |10 |y2_with_y2       |    6.552|1.286 |Residual         |
#> |11 |y3_with_y3       |    5.363|1.053 |Residual         |
#> |12 |y4_with_y4       |    2.137|0.822 |Residual         |
#> |13 |x1_with_x1       |    0.081|0.02  |Residual         |
#> |14 |x2_with_x2       |    0.120|0.072 |Residual         |
#> |15 |x3_with_x3       |    0.467|0.089 |Residual         |
#> 
#> Model Fit: χ²(13) = 23.74, p = 0.034; CFI = 0.972; TLI = 0.955; RMSEA = 0.105
#> RMSEA is worse than desired (<.06)
mxStandardizeRAMpaths(jreg_umx_fit, SE = TRUE)[4, ]
#>                      name          label matrix   row   col Raw.Value    Raw.SE
#> 4 latent_variables.A[9,8] ind60_to_dem60      A dem60 ind60  1.439314 0.3769464
#>   Std.Value    Std.SE
#> 4 0.4604654 0.1010676

Combined with IRT

Example from Lai & Hsiao (2021, Psychological Methods)

Not accounting for error

# Simulate data with mirt
set.seed(1235)
num_obs <- 1000
# Simulate theta
eta <- MASS::mvrnorm(num_obs, mu = c(0, 0), Sigma = diag(c(1, 1 - 0.5^2)),
                     empirical = TRUE)
th1 <- eta[, 1]
th2 <- -1 + 0.5 * th1 + eta[, 2]
# items and response data
a1 <- matrix(1, 10)
d1 <- matrix(rnorm(10))
a2 <- matrix(runif(10, min = 0.5, max = 1.5))
d2 <- matrix(rnorm(10))
dat1 <- simdata(a = a1, d = d1, N = num_obs, itemtype = "2PL", Theta = th1)
dat2 <- simdata(a = a2, d = d2, N = num_obs, itemtype = "2PL", Theta = th2)
# Factor scores
mod1 <- mirt(dat1, model = 1, itemtype = "Rasch", verbose = FALSE)
mod2 <- mirt(dat2, model = 1, itemtype = "2PL", verbose = FALSE)
fs1 <- fscores(mod1, full.scores.SE = TRUE)
fs2 <- fscores(mod2, full.scores.SE = TRUE)
lm(fs2[, 1] ~ fs1[, 1])  # attenuated coefficient
#> 
#> Call:
#> lm(formula = fs2[, 1] ~ fs1[, 1])
#> 
#> Coefficients:
#> (Intercept)     fs1[, 1]  
#>  -3.730e-06    2.812e-01

Joint model

With WLS

dat <- cbind(dat1, dat2)
colnames(dat) <- paste0("i", 1:20)
wls_fit <- sem("
f1 =~ i1 + i2 + i3 + i4 + i5 + i6 + i7 + i8 + i9 + i10
f2 =~ i11 + i12 + i13 + i14 + i15 + i16 + i17 + i18 + i19 + i20
f2 ~ f1
", data = dat, ordered = TRUE, std.lv = TRUE)
coef(wls_fit)["f2~f1"]
#>     f2~f1 
#> 0.5490816

With Maximum Likelihood

Note: This is extremely computational intensive

Using `tspa_mx_model()`

Because EAP (shrinkage) scores are used, we use the following properties for these scores ( $\tilde \eta_i$ for person $i$ ):

Standard error (SE, from software output) = $\sqrt{V(\eta) (1 - \rho_{\tilde \eta, i})}$ , where $\rho_{\tilde \eta, i}$ is the reliability of $\tilde \eta_i$
Loading of $\tilde \eta_i$ on $\eta$ is equal to $\rho_{\tilde \eta, i}$ = $1 - \text{SE}^2 / V(\eta)$
Standard error of measurement for $\tilde \eta_i$ is equal to $\rho_{\tilde \eta, i} \text{SE}$
Based on the above, the total variance of $\tilde \eta_i$ is also $\rho_{\tilde \eta, i}$ .

# Combine into data set
fs_dat <- cbind(fs1, fs2) |>
    as.data.frame() |>
    setNames(c("fs1", "se_fs1", "fs2", "se_fs2")) |>
    # Compute reliability and error variances
    within(expr = {
        rel_fs1 <- 1 - se_fs1^2
        rel_fs2 <- 1 - se_fs2^2
        ev_fs1 <- se_fs1^2 * (1 - se_fs1^2)
        ev_fs2 <- se_fs2^2 * (1 - se_fs2^2)
    })

# Loading
matL <- mxMatrix(
    type = "Diag", nrow = 2, ncol = 2,
    free = FALSE,
    labels = c("data.rel_fs2", "data.rel_fs1"),
    name = "L"
)
# Error
matE <- mxMatrix(
    type = "Diag", nrow = 2, ncol = 2,
    free = FALSE,
    labels = c("data.ev_fs2", "data.ev_fs1"),
    name = "E"
)

fsreg_umx <- umxLav2RAM(
    "
      fs2 ~ fs1
      fs2 + fs1 ~ 1
    ",
    printTab = FALSE)
#> 
#> ?plot.MxModel options: std, means, digits, strip_zero, file, splines=T/F/ortho,..., min=, max =, same = , fixed, resid= 'circle|line|none'

tspa_mx <- tspa_mx_model(fsreg_umx, data = fs_dat,
                         mat_ld = matL, mat_ev = matE)
# Run OpenMx
tspa_mx_fit <- mxRun(tspa_mx)
#> Running 2SPAD with 5 parameters
# Summarize the results
summary(tspa_mx_fit)
#> Summary of 2SPAD 
#>  
#> free parameters:
#>           name matrix row col     Estimate  Std.Error A
#> 1   fs1_to_fs2   m1.A fs2 fs1  0.500417224 0.05514849  
#> 2 fs2_with_fs2   m1.S fs2 fs2  0.778233329 0.07526000  
#> 3 fs1_with_fs1   m1.S fs1 fs1  0.905084857 0.06795953  
#> 4   one_to_fs2   m1.M   1 fs2 -0.000463287 0.04094454  
#> 5   one_to_fs1   m1.M   1 fs1 -0.001261722 0.03807646  
#> 
#> Model Statistics: 
#>                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
#>        Model:              5                   3995              4644.151
#>    Saturated:             NA                     NA                    NA
#> Independence:             NA                     NA                    NA
#> Number of observations/statistics: 1000/4000
#> 
#> Information Criteria: 
#>       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
#> AIC:      -3345.849               4654.151                 4654.212
#> BIC:     -22952.331               4678.690                 4662.810
#> 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-09-06 04:19:59 
#> Wall clock time: 0.7556674 secs 
#> optimizer:  SLSQP 
#> OpenMx version number: 2.21.12 
#> Need help?  See help(mxSummary)

Alternatively, one can use named matrices (mat_ld and mat_ev) to specify the names of the columns for the loading and error covariance matrices.

cross_load <- matrix(c("rel_fs2", NA, NA, "rel_fs1"), nrow = 2) |>
    `dimnames<-`(rep(list(c("fs2", "fs1")), 2))
err_cov <- matrix(c("ev_fs2", NA, NA, "ev_fs1"), nrow = 2) |>
    `dimnames<-`(rep(list(c("fs2", "fs1")), 2))
tspa_mx <- tspa_mx_model(fsreg_umx, data = fs_dat,
                         mat_ld = cross_load, mat_ev = err_cov)
# Run OpenMx
tspa_mx_fit <- mxRun(tspa_mx)
#> Running 2SPAD with 5 parameters
# Summarize the results
summary(tspa_mx_fit)
#> Summary of 2SPAD 
#>  
#> free parameters:
#>           name matrix row col     Estimate  Std.Error A
#> 1   fs1_to_fs2   m1.A fs2 fs1  0.500417224 0.05514849  
#> 2 fs2_with_fs2   m1.S fs2 fs2  0.778233329 0.07526000  
#> 3 fs1_with_fs1   m1.S fs1 fs1  0.905084857 0.06795953  
#> 4   one_to_fs2   m1.M   1 fs2 -0.000463287 0.04094454  
#> 5   one_to_fs1   m1.M   1 fs1 -0.001261722 0.03807646  
#> 
#> Model Statistics: 
#>                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
#>        Model:              5                   3995              4644.151
#>    Saturated:             NA                     NA                    NA
#> Independence:             NA                     NA                    NA
#> Number of observations/statistics: 1000/4000
#> 
#> Information Criteria: 
#>       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
#> AIC:      -3345.849               4654.151                 4654.212
#> BIC:     -22952.331               4678.690                 4662.810
#> 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-09-06 04:20:00 
#> Wall clock time: 0.7746925 secs 
#> optimizer:  SLSQP 
#> OpenMx version number: 2.21.12 
#> Need help?  See help(mxSummary)

Compare to joint model with OpenMx

Note: FIML couldn’t finish within 12 hours

jreg_umx_fit <- umxRAM(
    "
    f1 =~ a * i1 + a * i2 + a * i3 + a * i4 + a * i5 +
          a * i6 + a * i7 + a * i8 + a* i9 + a * i10
    f2 =~ i11 + i12 + i13 + i14 + i15 + i16 + i17 + i18 + i19 + i20
    f2 ~ f1
    i1 + i2 + i3 + i4 + i5 + i6 + i7 + i8 + i9 + i10 ~ 0
    i11 + i12 + i13 + i14 + i15 + i16 + i17 + i18 + i19 + i20 ~ 0
    i1 ~~ 1 * i1
    i2 ~~ 1 * i2
    i3 ~~ 1 * i3
    i4 ~~ 1 * i4
    i5 ~~ 1 * i5
    i6 ~~ 1 * i6
    i7 ~~ 1 * i7
    i8 ~~ 1 * i8
    i9 ~~ 1 * i9
    i10 ~~ 1 * i10
    i11 ~~ 1 * i11
    i12 ~~ 1 * i12
    i13 ~~ 1 * i13
    i14 ~~ 1 * i14
    i15 ~~ 1 * i15
    i16 ~~ 1 * i16
    i17 ~~ 1 * i17
    i18 ~~ 1 * i18
    i19 ~~ 1 * i19
    i20 ~~ 1 * i20
    ",
    data = lapply(as.data.frame(dat), FUN = mxFactor,
                 levels = c(0, 1)) |>
               as.data.frame(),
    type = "DWLS")
#> 2 latent variables were created:f1, f2.
#> object 'threshMat' created to handle: 20 binary variables:'i1', 'i2', 'i3', 'i4', 'i5', 'i6', 'i7', 'i8', 'i9', 'i10', 'i11', 'i12', 'i13', 'i14', 'i15', 'i16', 'i17', 'i18', 'i19', and 'i20'
#> 20 trait(s) are binary: 'i1', 'i2', 'i3', 'i4', 'i5', 'i6', 'i7', 'i8', 'i9', 'i10', 'i11', 'i12', 'i13', 'i14', 'i15', 'i16', 'i17', 'i18', 'i19', and 'i20'
#> For these, you you MUST fix the mean and variance of the latent traits driving each variable (usually 0 & 1 respectively) .
#> See ?umxThresholdMatrix
#> Using deviation-based model: Thresholds will be in'threshMat' based on deviations in 'deviations_for_thresh'
#> Running m1 with 32 parameters
#> ?umxSummary options: std=T|F', digits=, report= 'html', filter= 'NS' & more
#> 
#> 
#> Table: Parameter loadings for model 'm1'
#> 
#> |   |name         | Estimate|SE    |type             |
#> |:--|:------------|--------:|:-----|:----------------|
#> |1  |a            |    1.000|0     |Factor loading   |
#> |12 |f2_to_i11    |    1.000|0     |Factor loading   |
#> |13 |f2_to_i12    |    0.784|0.254 |Factor loading   |
#> |14 |f2_to_i13    |    1.805|0.4   |Factor loading   |
#> |15 |f2_to_i14    |    1.268|0.277 |Factor loading   |
#> |16 |f2_to_i15    |    1.935|0.42  |Factor loading   |
#> |17 |f2_to_i16    |    2.137|0.453 |Factor loading   |
#> |18 |f2_to_i17    |    1.279|0.289 |Factor loading   |
#> |19 |f2_to_i18    |    0.824|0.198 |Factor loading   |
#> |20 |f2_to_i19    |    0.730|0.194 |Factor loading   |
#> |21 |f2_to_i20    |    0.941|0.238 |Factor loading   |
#> |11 |f1_to_f2     |    0.335|0.068 |Factor to factor |
#> |42 |f1_with_f1   |    0.346|0.026 |Factor Variance  |
#> |43 |f2_with_f2   |    0.131|0.048 |Factor Variance  |
#> |22 |i1_with_i1   |    1.000|0     |Residual         |
#> |23 |i2_with_i2   |    1.000|0     |Residual         |
#> |24 |i3_with_i3   |    1.000|0     |Residual         |
#> |25 |i4_with_i4   |    1.000|0     |Residual         |
#> |26 |i5_with_i5   |    1.000|0     |Residual         |
#> |27 |i6_with_i6   |    1.000|0     |Residual         |
#> |28 |i7_with_i7   |    1.000|0     |Residual         |
#> |29 |i8_with_i8   |    1.000|0     |Residual         |
#> |30 |i9_with_i9   |    1.000|0     |Residual         |
#> |31 |i10_with_i10 |    1.000|0     |Residual         |
#> |32 |i11_with_i11 |    1.000|0     |Residual         |
#> |33 |i12_with_i12 |    1.000|0     |Residual         |
#> |34 |i13_with_i13 |    1.000|0     |Residual         |
#> |35 |i14_with_i14 |    1.000|0     |Residual         |
#> |36 |i15_with_i15 |    1.000|0     |Residual         |
#> |37 |i16_with_i16 |    1.000|0     |Residual         |
#> |38 |i17_with_i17 |    1.000|0     |Residual         |
#> |39 |i18_with_i18 |    1.000|0     |Residual         |
#> |40 |i19_with_i19 |    1.000|0     |Residual         |
#> |41 |i20_with_i20 |    1.000|0     |Residual         |
#> 
#> Model Fit: χ²(178) = 264.48, p < 0.001; CFI = NA; TLI = NA; RMSEA = 0.022
#> Algebra'threshMat':
#> 
#> 
#> |     |    i1|     i2|     i3|    i4|     i5|    i6|     i7|     i8|    i9|   i10|  i11|   i12|   i13|    i14|   i15|   i16|   i17|    i18|    i19|   i20|
#> |:----|-----:|------:|------:|-----:|------:|-----:|------:|------:|-----:|-----:|----:|-----:|-----:|------:|-----:|-----:|-----:|------:|------:|-----:|
#> |th_1 | 0.193| -0.801| -0.231| 0.134| -1.027| 0.336| -1.159| -0.058| 0.469| 0.149| 1.13| 1.484| 0.643| -0.082| 0.953| 0.702| 0.689| -0.026| -0.065| 0.581|
jreg_std <- mxStandardizeRAMpaths(jreg_umx_fit, SE = TRUE)
jreg_std[jreg_std$label == "f1_to_f2", ]
#>           name    label matrix row col Raw.Value     Raw.SE Std.Value
#> 11 m1.A[22,21] f1_to_f2      A  f2  f1 0.3349439 0.06796137  0.478655
#>        Std.SE
#> 11 0.04374742

Multidimensional Measurement Model

When a multidimensional model is used, for EAP scores, we use the following properties generalized from the unidimensional case. Specifically, let $V(\boldsymbol \eta)$ = $\boldsymbol{\psi}$ .

Loading of $\tilde{\boldsymbol \eta}_i$ on $\boldsymbol \eta$ is $\boldsymbol{\Lambda}_{\tilde \eta, i}$ , which is also the reliability matrix (i.e., the squared covariance between $\tilde{\boldsymbol \eta}_i$ and $\boldsymbol \eta$ ).
Error matrix ( $\mathbf{E}$ , from software output) = $(\mathbf{I} - \boldsymbol{\Lambda}_{\tilde \eta, i})\boldsymbol{\psi}$ . Therefore, $\boldsymbol{\Lambda}_{\tilde \eta, i} = \mathbf{I} - \mathbf{E} \boldsymbol{\psi}^{-1}$ .
Error variance of $\tilde \eta_i$ is $\boldsymbol{\Lambda}_{\tilde \eta, i} \mathbf{E}$ .

dat <- cbind(dat1, dat2)
colnames(dat) <- paste0("Item_", 1:20)
# Multidimensional IRT
mod <- "
F1 = 1-10
F2 = 11-20
COV = F1*F2
CONSTRAIN = (1-10, a1)
"
mfit <- mirt(dat, model = mod, itemtype = "2PL", verbose = FALSE)

# Factor scores
fs <- fscores(mfit)
# Error portion of factor scores
fs_ecov <- fscores(mfit, return.acov = TRUE)
# Latent variable covariance
psi <- diag(2)
psi[1, 2] <- psi[2, 1] <- coef(mfit)$GroupPars[1, "COV_21"]
inv_psi <- solve(psi)
load <- lapply(fs_ecov, FUN = function(err) {
    diag(2) - err %*% inv_psi
})
fs_err <- Map(`%*%`, load, fs_ecov)
# Convert to data frame
load_dat <- lapply(load, FUN = c) |>
    do.call(what = rbind) |>
    as.data.frame() |>
    setNames(c("F1_by_fs1", "F1_by_fs2", "F2_by_fs1", "F2_by_fs2"))
fs_err_dat <- lapply(fs_err, FUN = `[`, c(1, 2, 4)) |>
    do.call(what = rbind) |>
    as.data.frame() |>
    setNames(c("ev_fs1", "ecov_fs1_fs2", "ev_fs2"))

# Combine into data set
fs_dat <- cbind(fs, load_dat, fs_err_dat)
# Build OpenMx model
fsreg_umx <- umxLav2RAM(
    "
      F2 ~ F1
      F2 + F1 ~ 1
    ",
    printTab = FALSE)
#> 
#> ?plot.MxModel options: std, means, digits, strip_zero, file, splines=T/F/ortho,..., min=, max =, same = , fixed, resid= 'circle|line|none'
cross_load <- matrix(c("F1_by_fs1", "F1_by_fs2",
                       "F2_by_fs1", "F2_by_fs2"), nrow = 2) |>
    `dimnames<-`(rep(list(c("F1", "F2")), 2))
err_cov <- matrix(c("ev_fs1", "ecov_fs1_fs2",
                    "ecov_fs1_fs2", "ev_fs2"), nrow = 2) |>
    `dimnames<-`(rep(list(c("F1", "F2")), 2))
tspa_mx2 <- tspa_mx_model(fsreg_umx, data = fs_dat,
                          mat_ld = cross_load, mat_ev = err_cov)
# Run OpenMx
tspa_mx2_fit <- mxRun(tspa_mx2)
#> Running 2SPAD with 5 parameters
# Summarize the results
summary(tspa_mx2_fit)
#> Summary of 2SPAD 
#>  
#> free parameters:
#>         name matrix row col      Estimate  Std.Error A
#> 1   F1_to_F2   m1.A  F2  F1  4.753222e-01 0.05207370  
#> 2 F2_with_F2   m1.S  F2  F2  7.740976e-01 0.07360716  
#> 3 F1_with_F1   m1.S  F1  F1  1.000064e+00 0.07398508  
#> 4  one_to_F2   m1.M   1  F2  5.438812e-05 0.04081612  
#> 5  one_to_F1   m1.M   1  F1 -1.312738e-06 0.03972658  
#> 
#> Model Statistics: 
#>                |  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
#>        Model:              5                   3995              4327.057
#>    Saturated:             NA                     NA                    NA
#> Independence:             NA                     NA                    NA
#> Number of observations/statistics: 1000/4000
#> 
#> Information Criteria: 
#>       |  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
#> AIC:      -3662.943               4337.057                 4337.117
#> BIC:     -23269.425               4361.596                 4345.716
#> 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-09-06 04:20:11 
#> Wall clock time: 0.9092374 secs 
#> optimizer:  SLSQP 
#> OpenMx version number: 2.21.12 
#> Need help?  See help(mxSummary)
# Standardize coefficients
mxStandardizeRAMpaths(tspa_mx2_fit, SE = TRUE)$m1[1, ]
#> Warning in mxStandardizeRAMpaths(tspa_mx2_fit, SE = TRUE): 'model' (or one of
#> its submodels) contains definition variables; interpret results of
#> mxStandardizeRAMpaths() cautiously
#>        name    label matrix row col Raw.Value    Raw.SE Std.Value     Std.SE
#> 1 m1.A[1,2] F1_to_F2      A  F2  F1 0.4753222 0.0520737 0.4753271 0.04563049