modsem implements a Monte-Carlo correction for LMS and
QML models with ordinal data. Here we refer to these informally as
MC-LMS-ORD and MC-QML-ORD.
The MC-LMS-ORD and MC-QML-ORD algorithms
are based on Slupphaug, Mehmetoglu,
and Mittner (2026) For a more direct implementation of the original
algorithm, we recommend checking out the plssem
package.
Example
Here we ordinalize the data in the oneInt dataset.
ordinalize <- function(x, probs = c(0, 0.35, 0.7, 1)) {
x <- (x - mean(x)) / sd(x)
cut(
x,
breaks = stats::quantile(x, probs = probs),
include.lowest = TRUE,
ordered_result = TRUE
)
}
oneIntOrd <- as.data.frame(lapply(oneInt, ordinalize))Now we can estimate our model, indicating which variables are
ordinal, using the ordered= argument.
model <- "
X =~ x1 + x2 + x3
Z =~ z1 + z2 + z3
Y =~ y1 + y2 + y3
Y ~ X + Z + X:Z
"
fit_lms_ord <- modsem(
model,
data = oneIntOrd,
method = "lms",
ordered = colnames(oneIntOrd)
)
#> Warning: modsem->emLms():
#> The largest absolute gradient is 0.00173 (`convergence.grad = 0.001`, no
#> progress). The model might not have converged, or it might be poorly
#> identified. Consider increasing `em.control$grad.max.steps`.
summary(fit_lms_ord)
#> Warning: modsem->emLms():
#> The largest absolute gradient is 0.00118 (`convergence.grad = 0.001`, no
#> progress). The model might not have converged, or it might be poorly
#> identified. Consider increasing `em.control$grad.max.steps`.
#> Warning: modsem->emLms():
#> The largest absolute gradient is 0.00119 (`convergence.grad = 0.001`, no
#> progress). The model might not have converged, or it might be poorly
#> identified. Consider increasing `em.control$grad.max.steps`.
#> Warning: modsem->emLms():
#> The largest absolute gradient is 0.00125 (`convergence.grad = 0.001`, no
#> progress). The model might not have converged, or it might be poorly
#> identified. Consider increasing `em.control$grad.max.steps`.
#>
#> modsem (1.0.24) ended normally after 57 iterations
#>
#> Estimator MC-LMS
#> Optimization method ROBBINS-MONRO
#> Naive optimization method EMA-NLMINB
#> Number of model parameters 13
#>
#> Number of observations 2000
#>
#> Naive Loglikelihood and Information Criteria:
#> Naive Loglikelihood -19223.44
#> Naive Akaike (AIC) 38472.87
#> Naive Bayesian (BIC) 38545.68
#>
#> Numerical Integration:
#> Standard Effective
#> Points of integration (per dim) 24 112
#> Dimensions 1 1
#> Total points of integration 24 112
#>
#> Fit Measures for Baseline Model (H0):
#> Standard
#> Chi-square 30.56
#> Degrees of Freedom (Chi-square) 33
#> P-value (Chi-square) 0.589
#> RMSEA 0.000
#>
#> Naive Loglikelihood -19343.47
#> Naive Akaike (AIC) 38710.94
#> Naive Bayesian (BIC) 38778.15
#>
#> Comparative Fit to H0 (LRT test):
#> Loglikelihood change 120.04
#> Difference test (D) 240.07
#> Degrees of freedom (D) 1
#> P-value (D) 0.000
#>
#>
#> Parameter Estimates:
#> Coefficients standardized
#> Information observed
#> Standard errors mc-delta
#>
#> Latent Variables:
#> Estimate Std.Error z.value P(>|z|)
#> X =~
#> x1 0.932 0.007 140.787 0.000
#> x2 0.900 0.007 130.155 0.000
#> x3 0.919 0.007 127.371 0.000
#> Z =~
#> z1 0.919 0.008 112.824 0.000
#> z2 0.904 0.007 127.140 0.000
#> z3 0.908 0.008 108.723 0.000
#> Y =~
#> y1 0.975 0.004 264.950 0.000
#> y2 0.957 0.004 213.621 0.000
#> y3 0.964 0.004 233.808 0.000
#>
#> Regressions:
#> Estimate Std.Error z.value P(>|z|)
#> Y ~
#> X 0.434 0.020 21.320 0.000
#> Z 0.380 0.021 17.985 0.000
#> X:Z 0.485 0.027 17.996 0.000
#>
#> Covariances:
#> Estimate Std.Error z.value P(>|z|)
#> X ~~
#> Z 0.190 0.027 7.124 0.000
#>
#> Thresholds:
#> Estimate Std.Error z.value P(>|z|)
#> x1|t1 -0.387 0.027 -14.289 0.000
#> x1|t2 0.533 0.031 17.069 0.000
#> x2|t1 -0.387 0.030 -13.022 0.000
#> x2|t2 0.529 0.030 17.727 0.000
#> x3|t1 -0.394 0.029 -13.646 0.000
#> x3|t2 0.522 0.029 17.732 0.000
#> z1|t1 -0.390 0.030 -13.151 0.000
#> z1|t2 0.525 0.030 17.745 0.000
#> z2|t1 -0.381 0.030 -12.825 0.000
#> z2|t2 0.520 0.029 17.657 0.000
#> z3|t1 -0.393 0.031 -12.847 0.000
#> z3|t2 0.524 0.030 17.606 0.000
#> y1|t1 -0.322 0.025 -13.111 0.000
#> y1|t2 0.496 0.030 16.654 0.000
#> y2|t1 -0.326 0.025 -13.127 0.000
#> y2|t2 0.508 0.030 16.818 0.000
#> y3|t1 -0.323 0.026 -12.538 0.000
#> y3|t2 0.505 0.030 17.064 0.000
#>
#> Variances:
#> Estimate Std.Error z.value P(>|z|)
#> .x1 0.132
#> .x2 0.190
#> .x3 0.156
#> .z1 0.155
#> .z2 0.182
#> .z3 0.175
#> .y1 0.049
#> .y2 0.084
#> .y3 0.070
#> X 1.000
#> Z 1.000
#> .Y 0.369