In models with ordinal indicators, the response proportions per level of the ordinal factors is often modeled as a “threshold”: a value of the standard normal distribution where the probability mass in the left tail corresponds to the response proportion of that level of the ordinal factor.
An alternative way to parametrize ordinal responses is to model response proportions per level of the ordinal factors in terms of threshold differences, which are strictly positive. This benefits model convergence, especially when for example estimating an LCA model. The thresholds are then computed by summing successive deviances.
The choice of parameterization does not change the substantive meaning of the thresholds, but it does change the model parameters. Consequently, it also changes how hypotheses must be written for a Wald test.
This vignette walks through two possible specifications using
as_ram():
model.Thresholds; andthreshold_method = "mx_deviances", in which thresholds
are constructed from deviation parameters stored in
model.mat_dev.The examples use the same simulated data and the same threshold model
under both parameterizations. The main objective is to show how to
translate a substantive hypothesis about an ordinal threshold into the
parameter expression expected by wald_test().
We begin by simulating a small data set with three ordinal variables.
library(tidySEM)
library(OpenMx)
set.seed(1)
df <- data.frame(rbind(
matrix(rnorm(300, mean = 2), ncol = 3),
matrix(rnorm(600), ncol = 3)
))
df$X1 <- ordered(cut(df$X1, 2, labels = FALSE))
df$X2 <- ordered(cut(df$X2, 2, labels = FALSE))
df$X3 <- ordered(cut(df$X3, 3, labels = FALSE))set.seed(1) makes the simulation reproducible. The first
block contributes 100 rows for three variables with mean 2, whereas the
second contributes 200 rows with mean 0. After rbind(), the
data set therefore contains 300 observations on three columns, named
X1, X2, and X3 by
data.frame(). The first two variables are converted to
ordered factors with two categories; X3 is converted to an
ordered factor with three categories. A two-category ordinal variable
requires one threshold to separate those whose score is in category 1
from those in category 2. A three-category ordinal variable requires two
thresholds. The model therefore contains four thresholds in total: one
for X1, one for X2, and two for
X3.
as_ram()The threshold model is specified with the |
operator.
Each line has the form variable | threshold_label, where
threshold labels are always t[number of the threshold]. The
model asks for the first threshold (t1) of X1
and X2, and both the first (t1) and second
(t2) thresholds of X3.
The labels t1 and t2 identify threshold
positions within a variable. They do not imply that
thresholds carrying the same label are constrained to be equal across
variables. For example, the first threshold of X1 and the
first threshold of X2 are distinct model parameters unless
an equality constraint is imposed explicitly, as in:
With the default threshold method, the parameters are thresholds,
which you can verify by running coef(mod).
The model is estimated with run_mx().
run_mx() fits the OpenMx model created by
as_ram() and returns the fitted model object. A tidy
summary can be obtained with:
table_results(res)
#> label est_sig se pval confint
#> 1 Variances.X1 1.00 <NA> <NA> <NA>
#> 2 Variances.X2 1.00 <NA> <NA> <NA>
#> 3 Variances.X3 1.00 <NA> <NA> <NA>
#> 4 model.Thresholds[1,1] 0.06 0.07 0.42 [-0.08, 0.20]
#> 5 model.Thresholds[1,2] 0.11 0.07 0.13 [-0.03, 0.25]
#> 6 model.Thresholds[1,3] -0.52*** 0.08 0.00 [-0.67, -0.38]
#> 7 model.Thresholds[2,3] 1.01*** 0.09 0.00 [0.84, 1.18]For ordinal indicators, the output includes the estimated threshold parameters along with their standard errors and other model results. This table is usually the most convenient place to inspect the estimates substantively.
For Wald tests, however, it is important to know how those estimates
are represented internally. Under the default specification, hypotheses
can refer directly to entries of model.Thresholds.
The first Wald test is:
wald_test(res, "model.Thresholds[1,1] = 0")
#> Wald tests for linear hypotheses:
#> Hypothesis df chisq p
#> model.Thresholds[1,1] = 0 1 0.65 0.42Here, the first matrix index refers to the threshold position and the
second to the observed variable. Thus,
model.Thresholds[1,1] is the first threshold of the first
ordinal indicator, X1.
wald_test() evaluates the stated equality using the
fitted estimate and its estimated sampling covariance. A small p-value
indicates that the threshold differs significantly from zero according
to the Wald approximation.
The second test is:
wald_test(res, "model.Thresholds[2,3] = 1")
#> Wald tests for linear hypotheses:
#> Hypothesis df chisq p
#> model.Thresholds[2,3] = 1 1 0.0089 0.92X3 is the third observed variable, and its second
threshold occupies row 2, column 3 of the threshold matrix. This direct
indexing is the main convenience of the default parameterization: a
hypothesis about a threshold can be written as a hypothesis about the
corresponding entry of model.Thresholds.
threshold_method = "mx_deviances"The same substantive model can be constructed with a different threshold parameterization:
The observed variables and requested thresholds are unchanged. What
changes is the way the thresholds are parameterized inside OpenMx.
Verify that these are now defined as deviances by running
coef(mod2).
With threshold_method = "mx_deviances", the model uses a
matrix of deviation parameters, model.mat_dev. Thresholds
are then constructed by cumulatively summing these deviations. This
cumulative representation is useful for encoding the ordered structure
of thresholds. It also means that the second and later thresholds are
generally functions of multiple model parameters,
rather than single free parameters.
Fit and summarize this model in the same way:
res2 <- run_mx(mod2)
table_results(res2)
#> label est_sig se pval confint
#> 1 Thresholds[1,1] 0.06 0.07 0.42 [-0.08, 0.20]
#> 2 Thresholds[1,2] 0.11 0.07 0.13 [-0.03, 0.25]
#> 3 Thresholds[1,3] -0.52*** 0.08 0.00 [-0.67, -0.38]
#> 4 Thresholds[2,3] 1.01*** 0.09 0.00 [0.84, 1.18]
#> 5 Variances.X1 1.00 <NA> <NA> <NA>
#> 6 Variances.X2 1.00 <NA> <NA> <NA>
#> 7 Variances.X3 1.00 <NA> <NA> <NA>The substantive threshold estimates reported by
table_results() are already reconstructed from the
deviances, and can be interpreted in the same way as before. The
difference becomes important when referring to the free parameters
directly.
For a variable with only one threshold, the first deviation is equal to the threshold. For a variable with multiple thresholds, later thresholds must be reconstructed by adding the relevant deviations.
Because the first threshold is equal to the first deviation, the corresponding test for deviances is:
wald_test(res2, "model.mat_dev[1,1] = 0")
#> Wald tests for linear hypotheses:
#> Hypothesis df chisq p
#> model.mat_dev[1,1] = 0 1 0.65 0.42Verify that the test result is identical to:
The difference between the parameterizations is clearest for the
second threshold of X3. It must be written as a constraint
on the sum of the first two deviations:
wald_test(res2, "model.mat_dev[1,3]+model.mat_dev[2,3] = 1")
#> Wald tests for linear hypotheses:
#> Hypothesis df chisq p
#> model.mat_dev[1,3]+model.mat_dev[2,3] = 1 1 0.0089 0.92This is the deviation-parameter counterpart of:
wald_test(res, "model.Thresholds[2,3] = 1")
#> Wald tests for linear hypotheses:
#> Hypothesis df chisq p
#> model.Thresholds[2,3] = 1 1 0.0089 0.92Do not make the mistake of testing only
model.mat_dev[2,3] = 1. That would test whether the
second deviation equals 1, not whether the second
threshold equals 1. Those are different hypotheses because the
second threshold is cumulative.
The default and deviation specifications are alternative parameterizations of the same thresholds. When both models represent the same statistical model and converge to the same solution, a Wald test formulated for the same substantive threshold should lead to equivalent conclusions, apart from small numerical differences.
The expressions passed to wald_test() differ because a
Wald test operates on the parameterization of the fitted model.
Reparameterizing a threshold as a sum of deviations therefore requires
re-expressing the null hypothesis in terms of that sum.
This distinction is especially important for second and later
thresholds. In the direct specification, the target threshold is a
single entry of model.Thresholds; in the deviation
specification, the same target is a linear combination of entries in
model.mat_dev.
While strictly positive deviances may aid in achieving model
convergence, once convergence has been reached, it is possible to
reparametrize the model to use thresholds instead. At present, this is
supported by the utility function
deviances_to_thresholds(), which takes in a model with
deviances and returns a model with thresholds. The model is
automatically re-run.
Use this functionality as follows:
coef(res2)
#> model.mat_dev[1,1] model.mat_dev[1,2] model.mat_dev[1,3] model.mat_dev[2,3]
#> 0.059 0.109 -0.524 1.533
res2_thresholds <- deviances_to_thresholds(res2)
coef(res2_thresholds)
#> model.Thresholds[1,1] model.Thresholds[1,2] model.Thresholds[1,3]
#> 0.059 0.109 -0.524
#> model.Thresholds[2,3]
#> 1.008
wald_test(res2_thresholds, "model.Thresholds[1,1] = 0")
#> Wald tests for linear hypotheses:
#> Hypothesis df chisq p
#> model.Thresholds[1,1] = 0 1 0.65 0.42