Score and criterion tests

What do you want to test?

Suppose a factor model has two loadings, a and b, and we want to test whether they are equal. We can examine the score after imposing equality, or compare the fitting criteria with and without equality. Both are tests of a restriction within a model. Testing whether the whole model fits the population moments exactly is a separate question.

semtest() builds all of these tests. Its first argument is the model fitted under the null. The second, when present, is the alternative.

pvalues_nested() also compares nested models, but it tests a different null: that both models fit exactly, using Satorra’s (2000) construction. The nested tests in this vignette allow the larger model to be misspecified, as explained below, so the two need not agree.

We first work through the equality restriction using complete continuous data.

One restriction, two tests

The same factor model before and after imposing equality. Blue paths identify the restriction. Residuals are uncorrelated; intercepts and variances are omitted from the drawing.

The first loading is fixed to 1 in both models, so their factor scales agree. The loadings of x2 and x3 are labeled a and b. Our null is a = b; the loading of x4, variances, and intercepts remain free in both models.

d <- lavaan::HolzingerSwineford1939
model <- "f =~ x1 + a*x2 + b*x3 + x4"
fit <- lavaan::cfa(model, d, meanstructure = TRUE)
restricted <- lavaan::cfa(paste(model, "a == b", sep = "\n"),
                          d, meanstructure = TRUE)
lavaan::coef(fit)[c("a", "b")]
#>         a         b 
#> 0.5466181 0.7042767

score <- semtest(restricted, fit)
score
#> score test: Common criterion pseudo-target satisfies the restricted model 
#>   type statistic df rank      reference   p_value
#>  score  2.312464  1    1 weighted_chisq 0.1119368
#>                                                     hypothesis
#>  Common criterion pseudo-target satisfies the restricted model

template <- lavaan::cfa(model, d, meanstructure = TRUE, do.fit = FALSE)
semtest(restricted, template)
#> score test: Common criterion pseudo-target satisfies the restricted model 
#>   type statistic df rank      reference   p_value
#>  score  2.312464  1    1 weighted_chisq 0.1119368
#>                                                     hypothesis
#>  Common criterion pseudo-target satisfies the restricted model

criterion <- semtest(restricted, fit, statistic = "criterion")
criterion
#> criterion difference test: Equal population minimum under the shared criterion 
#>                  type statistic df rank      reference   p_value
#>  criterion difference  2.170501  1    1 weighted_chisq 0.1235645
#>                                           hypothesis
#>  Equal population minimum under the shared criterion

What if the factor model is misspecified?

A model can fail to reproduce the population moments even at its best-fitting parameter value. That value is the pseudo-parameter: the minimizer of the population fitting criterion. Under structural misspecification, a = b refers to its two loadings, rather than to parameters of a true factor model. The target can depend on the fitting criterion.

For nested inference, we require a shared criterion and, under the null, a common population minimizer in the restricted model. We also need the usual local regularity conditions, including identification and suitable curvature. A sandwich calculation does not remove those requirements.

It helps to distinguish four issues:

• Nonnormality: the observations need not follow a normal distribution, even when the model reproduces their means and covariances.
• Structural misspecification: the model cannot reproduce the target moments; restrictions concern the criterion’s pseudo-parameter.
• Missing data: FIML uses the observed-data likelihood contributions. Its population target also depends on the observation process and assumptions.
• Estimated fitting weights: for weighted least squares, the criterion’s weight matrix is estimated from the same data as the sample moments.

Satorra (1989) gives a unified treatment of restriction tests, and Satorra (2000) develops scaled and adjusted versions under nonnormality. The extensions to missing data, structural misspecification, and estimated fitting weights used here are developed in manuscripts in preparation.

How the statistic and its reference are formed

Two separate decisions shape a score test: which matrix weights its quadratic form, and how the underlying matrices are estimated. semtest() takes the first as its weight argument. The second belongs to inference_control(), whose settings can also be passed to semtest() directly.

Write \(s\) for the score at the restricted fit, \(H\) for the curvature of the fitting criterion (how the score changes with the parameters), and \(J\) for the covariance of the score. Profiling uses \(H\) to remove the directions in which the restricted model can still move. Call the profiled score \(\tilde s\), its covariance \(\tilde J\), and the profiled weight matrix \(\tilde M\). The statistic is \[ T = \tilde s^{\mathsf T} \tilde M^{-1} \tilde s . \] Its reference distribution is a weighted sum of independent \(\chi^2_1\) variables, weighted by the eigenvalues of \(\tilde M^{-1} \tilde J\). Here curvature means \(H\) itself, not its inverse.

weight chooses \(M\):

• "curvature", the default, uses \(H\) itself.
• "expected_curvature" uses the expected curvature even when \(H\) is the observed curvature. This keeps observed curvature in the profiling, which the limit requires when the model may be misspecified, while weighting by the more stable expected matrix.
• "score_covariance" uses \(J\). Every eigenvalue is then 1, and \(T\) has a chi-square reference: the score is fully studentized.

The criterion difference always weights by curvature, so weight does not apply to it.

With a single restriction there is a single eigenvalue, and all three weightings give the same p-value, although the statistics differ. The choice matters when several restrictions are tested jointly. Here all three free loadings are held equal:

model3 <- "f =~ x1 + a*x2 + b*x3 + c*x4"
fit3 <- lavaan::cfa(model3, d, meanstructure = TRUE)
equal3 <- lavaan::cfa(paste(model3, "a == b", "b == c", sep = "\n"),
                      d, meanstructure = TRUE)
weights <- c("curvature", "expected_curvature", "score_covariance")
do.call(rbind, lapply(weights, function(w) {
  test <- semtest(equal3, fit3, weight = w)
  data.frame(weight = w, statistic = test$statistic, df = test$df,
             reference = test$reference, p_value = test$p_value)
}))
#>               weight statistic df      reference   p_value
#> 1          curvature  2.901017  2 weighted_chisq 0.3163009
#> 2 expected_curvature  2.710859  2 weighted_chisq 0.2866476
#> 3   score_covariance  3.486518  2          chisq 0.1749493

inference_control()
#> semTests inference control
#>   curvature         <by test>
#>   score_covariance  <by test>
#>   weight_effect     estimated
#>   derivative_step   1e-04
#>   tolerance         1e-10
#>   keep              FALSE

semtest(equal3, fit3, curvature = "expected")
#> score test: Common criterion pseudo-target satisfies the restricted model 
#>   type statistic df rank      reference   p_value
#>  score  2.817358  2    2 weighted_chisq 0.2684397
#>                                                     hypothesis
#>  Common criterion pseudo-target satisfies the restricted model
semtest(equal3, fit3, control = inference_control(curvature = "expected"))
#> score test: Common criterion pseudo-target satisfies the restricted model 
#>   type statistic df rank      reference   p_value
#>  score  2.817358  2    2 weighted_chisq 0.2684397
#>                                                     hypothesis
#>  Common criterion pseudo-target satisfies the restricted model

test3 <- semtest(equal3, fit3)
test3$eigenvalues
#> [1] 1.7863363 0.8110898
pvalues(test3, c("ALL", "SB", "SS", "PEBA2"))
#>       all        sb        ss     peba2 
#> 0.3163009 0.3272989 0.3297492 0.3244743 
#> test: score | reference: weighted_chisq | df: 2

Global exact fit

Return now to the separate question: does the unrestricted factor model fit the population moments exactly?

global <- semtest(fit)
global
#> score test: Exact fit against saturated moments 
#>   type statistic df rank      reference     p_value
#>  score  12.02448  2    2 weighted_chisq 0.008817575
#>                           hypothesis
#>  Exact fit against saturated moments
pvalues(global, c("SB", "ALL"))
#>          sb         all 
#> 0.006658492 0.008817575 
#> test: score | reference: weighted_chisq | df: 2

This compares the fitted model with saturated means and covariances, using likelihood scores at the fitted moments. Rejection concerns the whole model, not specifically the equality of a and b. The robust reference handles nonnormality, but the null is still exact fit.

For overall fit, curvature refers to the curvature of the saturated model. The classical discrepancy statistic for the same question is what pvalues(fit) computes, so semtest() does not duplicate it: statistic = "criterion" requires an alternative.

Missing data and estimated fitting weights

pvalues() and pvalues_nested() keep their own scope; see ?semTests-support.

For FIML, likelihood scores are observed-data scores for each missingness pattern, and the score covariance is their empirical covariance, as for complete data. See vignette("fiml-missing-data", package = "semTests") for the missing-data workflow. "STD" gives the non-robust reference, which assumes the information equality:

d_missing <- d
set.seed(1)
d_missing$x2[sample(nrow(d), 60)] <- NA
fit_fiml <- lavaan::cfa(model, d_missing, meanstructure = TRUE,
                        missing = "fiml")
pvalues(semtest(fit_fiml), c("STD", "SB", "ALL"))
#>         std          sb         all 
#> 0.003311375 0.011381091 0.013294814 
#> test: score | reference: weighted_chisq | df: 2

For least squares, the fitting weights are the matrix weighting discrepancies between observed and model-implied moments. They are not sampling weights. Their estimation can contribute to sampling variability when the structural model is misspecified. The contribution vanishes at exact fit.

Here is the equality of all three loadings, fitted by continuous DWLS:

fit_dwls <- lavaan::cfa(model3, d, meanstructure = TRUE, estimator = "DWLS",
                        ordered = FALSE)
equal_dwls <- lavaan::cfa(paste(model3, "a == b", "b == c", sep = "\n"), d,
                          meanstructure = TRUE, estimator = "DWLS",
                          ordered = FALSE)
semtest(equal_dwls, fit_dwls, weight = "score_covariance")
#> score test: Common criterion pseudo-target satisfies the restricted model 
#>   type statistic df rank reference   p_value
#>  score  3.557322  2    2     chisq 0.1688641
#>                                                     hypothesis
#>  Common criterion pseudo-target satisfies the restricted model
semtest(equal_dwls, fit_dwls, weight = "score_covariance",
        weight_effect = "fixed")
#> score test: Common criterion pseudo-target satisfies the restricted model 
#>   type statistic df rank reference   p_value
#>  score  3.080126  2    2     chisq 0.2143676
#>                                                     hypothesis
#>  Common criterion pseudo-target satisfies the restricted model

The default, weight_effect = "estimated", includes the sampling variability of both the sample moments and the estimated fitting weights. "fixed" treats the weights as known, which is correct only when the model fits exactly. Least-squares tests require the centered score covariance, which is the nested default.

These tests currently exclude sampling weights, clusters, nonlinear constraints, and parameter bounds. The ordinal path requires complete, all-ordered data with nonempty pairwise cells; it does not support mixed indicators or pairwise missingness. See vignette("categorical-data", package = "semTests") for ordinal fitting.

Inspecting a result

Pass keep = TRUE to retain the numerical ingredients in $ingredients, and inspect $diagnostics and $provenance when checking a calculation. For example, the nuisance stationarity norm checks that the restricted fit is a stationary point of its own criterion:

score$diagnostics$nuisance_stationarity_norm
#> [1] 1.156657e-07

For ordinal data, the first-stage influence of the thresholds and polychoric correlations is computed in closed form, using lavaan’s score-product convention for the first-stage weight. So is the derivative of the estimated weight, for DWLS and full WLS alike.

References

Satorra, A. (1989). Alternative test criteria in covariance structure analysis: A unified approach. Psychometrika, 54, 131–151. https://doi.org/10.1007/BF02294453

Satorra, A. (2000). Scaled and adjusted restricted tests in multi-sample analysis of moment structures. In R. D. H. Heijmans, D. S. G. Pollock, & A. Satorra (Eds.), Innovations in Multivariate Statistical Analysis (pp. 233–247). Kluwer Academic. https://doi.org/10.1007/978-1-4615-4603-0_17