Uncertainty in SEM

Introduction

Fitting a Structural Equation Model (SEM) yields point estimates (e.g., a path coefficient of \(0.45\) or a variance component of \(0.21\)). However, scientific inference requires representing the uncertainty of these estimates. Moreover, we are interested in the fit of the model itself to the data.

umx (+ OpenMx) provides three primary ways to assess uncertainty: 1. Standard Maximum Likelihood Standard Errors (Wald-based). 2. Robust Standard Errors (Huber-White / “Sandwich” estimator). 3. Profile Likelihood Confidence Intervals (likelihood-based).

This vignette explains the differences between these methods, when to use each, and how to implement them in your umx workflow.

1. Point Estimates vs. Inference (SEs vs. CIs)

Standard Errors (Wald CIs)

Standard errors represent the standard deviation of a parameter’s sampling distribution. By default, packages construct Wald confidence intervals using the standard error: \[\hat{\theta} \pm z_{\alpha/2} \times \text{SE}\]

  • Assumption: The sampling distribution of the estimator is symmetric and normal.
  • Limitations: This assumption breaks down in several common scenarios:
    • Small sample sizes.
    • Parameters near boundaries (e.g., variances, which are bounded at 0).
    • Highly non-linear parameter functions (e.g., indirect effects in mediation models).

Profile Likelihood Confidence Intervals

Instead of assuming a normal distribution, Profile Likelihood Confidence Intervals construct the interval directly from the log-likelihood surface.

To find the lower boundary of a 95% profile CI for a parameter \(\theta\): 1. The parameter \(\theta\) is fixed to a value below its maximum likelihood estimate (MLE). 2. All other parameters in the model are re-optimized to find the maximum possible log-likelihood. 3. Steps 1-2 are repeated until we find the point where the log-likelihood drops by exactly \(\chi^2_1(0.95)/2 \approx 1.92\) units compared to the MLE.

Because they re-optimize the model at each step, profile CIs: * Respect the asymmetry of the parameter space (e.g., a CI for variance might be \([0.02, 0.45]\) rather than \([-0.10, 0.35]\)). * Do not rely on local curvature (Hessian matrix) at the MLE. * Restrict values to boundary constraints.

2. ML Standard Errors vs. Robust Standard Errors

Standard ML SEs

Standard Maximum Likelihood standard errors are derived from the local curvature of the log-likelihood surface at the peak (the inverse of the Hessian matrix / information matrix). * When to use: When the data are multivariate normal, the model is correctly specified, and the sample size is large.

Robust Standard Errors (Sandwich Estimator)

When data are non-normal, standard ML standard errors are typically underestimated, leading to inflated Type-I error rates. Robust standard errors (also called Huber-White or “sandwich” SEs) correct for this.

The robust covariance matrix is calculated as: \[\text{Cov}_{\text{robust}} = H^{-1} B H^{-1}\] * \(H\) (The “Bread”): The standard Hessian matrix (curvature of the likelihood surface). * \(B\) (The “Meat”): The variability of the individual case-wise gradients.

When to prefer Robust SEs:

  • Mild-to-moderate violation of multivariate normality.
  • Minor model miss-specifications.
  • Clustered or nested data (if using a cluster-robust correction).

Limitations:

  • Robust standard errors do not fix biased point estimates resulting from severe misspecification.
  • They can be unstable in small samples (\(N < 100\)).
  • They still rely on the Wald symmetry assumption for confidence intervals.

3. Practical Workflow in umx

Let’s illustrate these options using the built-in demoOneFactor dataset.

Step 1: Fit a Base Model

We fit a simple one-factor model using umxRAM(). By default, umxSummary() reports standard ML standard errors.

# Load data
data(demoOneFactor)
manifests = names(demoOneFactor)

# Fit model
m1 = umxRAM("One_Factor_Model", data = demoOneFactor, 
    umxPath("g", to = manifests),
    umxPath(v.m. = manifests),
    umxPath(v1m0 = "g")
)
#> 
#> 
#> Table: Parameter loadings for model 'One_Factor_Model'
#> 
#> |   |name       | Estimate|SE    |CI                      |type            |
#> |:--|:----------|--------:|:-----|:-----------------------|:---------------|
#> |1  |g_to_x5    |    0.795|0.027 |0.795 [0.743, 0.848]    |Factor loading  |
#> |2  |g_to_x4    |    0.702|0.024 |0.702 [0.655, 0.749]    |Factor loading  |
#> |3  |g_to_x3    |    0.577|0.02  |0.577 [0.537, 0.617]    |Factor loading  |
#> |4  |g_to_x2    |    0.503|0.018 |0.503 [0.467, 0.539]    |Factor loading  |
#> |5  |g_to_x1    |    0.397|0.016 |0.397 [0.366, 0.427]    |Factor loading  |
#> |11 |g_with_g   |    1.000|0     |1 [1, 1]                |Factor Variance |
#> |12 |one_to_x5  |   -0.076|0.037 |-0.076 [-0.147, -0.004] |Mean            |
#> |13 |one_to_x4  |   -0.056|0.033 |-0.056 [-0.12, 0.008]   |Mean            |
#> |14 |one_to_x3  |   -0.056|0.027 |-0.056 [-0.109, -0.002] |Mean            |
#> |15 |one_to_x2  |   -0.046|0.024 |-0.046 [-0.093, 0.001]  |Mean            |
#> |16 |one_to_x1  |   -0.040|0.02  |-0.04 [-0.079, -0.001]  |Mean            |
#> |6  |x5_with_x5 |    0.036|0.004 |0.036 [0.029, 0.043]    |Residual        |
#> |7  |x4_with_x4 |    0.039|0.003 |0.039 [0.033, 0.046]    |Residual        |
#> |8  |x3_with_x3 |    0.041|0.003 |0.041 [0.035, 0.047]    |Residual        |
#> |9  |x2_with_x2 |    0.038|0.003 |0.038 [0.032, 0.043]    |Residual        |
#> |10 |x1_with_x1 |    0.041|0.003 |0.041 [0.035, 0.046]    |Residual        |

We can inspect the standard ML standard errors and Wald CIs:

umxSummary(m1, std = TRUE, uncertainty = "SE")
#> 
#> 
#> Table: Parameter loadings for model 'One_Factor_Model'
#> 
#> |   |name       | Std.Estimate|Std.SE |CI                  |type            |
#> |:--|:----------|------------:|:------|:-------------------|:---------------|
#> |1  |g_to_x5    |         0.97|0      |0.97 [0.97, 0.98]   |Factor loading  |
#> |2  |g_to_x4    |         0.96|0      |0.96 [0.95, 0.97]   |Factor loading  |
#> |3  |g_to_x3    |         0.94|0.01   |0.94 [0.93, 0.95]   |Factor loading  |
#> |4  |g_to_x2    |         0.93|0.01   |0.93 [0.92, 0.95]   |Factor loading  |
#> |5  |g_to_x1    |         0.89|0.01   |0.89 [0.87, 0.91]   |Factor loading  |
#> |11 |g_with_g   |         1.00|0      |1 [1, 1]            |Factor Variance |
#> |12 |one_to_x5  |        -0.09|0.04   |-0.09 [-0.18, 0]    |Mean            |
#> |13 |one_to_x4  |        -0.08|0.04   |-0.08 [-0.16, 0.01] |Mean            |
#> |14 |one_to_x3  |        -0.09|0.04   |-0.09 [-0.18, 0]    |Mean            |
#> |15 |one_to_x2  |        -0.08|0.04   |-0.08 [-0.17, 0]    |Mean            |
#> |16 |one_to_x1  |        -0.09|0.04   |-0.09 [-0.18, 0]    |Mean            |
#> |6  |x5_with_x5 |         0.05|0.01   |0.05 [0.04, 0.07]   |Residual        |
#> |7  |x4_with_x4 |         0.07|0.01   |0.07 [0.06, 0.09]   |Residual        |
#> |8  |x3_with_x3 |         0.11|0.01   |0.11 [0.09, 0.13]   |Residual        |
#> |9  |x2_with_x2 |         0.13|0.01   |0.13 [0.11, 0.15]   |Residual        |
#> |10 |x1_with_x1 |         0.21|0.02   |0.21 [0.17, 0.24]   |Residual        |

Step 2: Request Robust Standard Errors

To obtain robust standard errors, we can use uncertainty = "RobustSE".

Compute robust standard errors

umxSummary(m1, std = TRUE, uncertainty = "RobustSE")
#> Warning: In model 'Saturated One_Factor_Model' Optimizer returned a non-zero
#> status code 6. The model does not satisfy the first-order optimality conditions
#> to the required accuracy, and no improved point for the merit function could be
#> found during the final linesearch (Mx status RED)
#> 
#> 
#> Table: Parameter loadings for model 'One_Factor_Model'
#> 
#> |   |name       | Std.Estimate|Std.SE |CI                  |type            |
#> |:--|:----------|------------:|:------|:-------------------|:---------------|
#> |1  |g_to_x5    |         0.97|0      |0.97 [0.97, 0.98]   |Factor loading  |
#> |2  |g_to_x4    |         0.96|0      |0.96 [0.95, 0.97]   |Factor loading  |
#> |3  |g_to_x3    |         0.94|0.01   |0.94 [0.93, 0.95]   |Factor loading  |
#> |4  |g_to_x2    |         0.93|0.01   |0.93 [0.92, 0.94]   |Factor loading  |
#> |5  |g_to_x1    |         0.89|0.01   |0.89 [0.87, 0.91]   |Factor loading  |
#> |11 |g_with_g   |         1.00|0      |1 [1, 1]            |Factor Variance |
#> |12 |one_to_x5  |        -0.09|0.04   |-0.09 [-0.18, 0]    |Mean            |
#> |13 |one_to_x4  |        -0.08|0.04   |-0.08 [-0.16, 0.01] |Mean            |
#> |14 |one_to_x3  |        -0.09|0.04   |-0.09 [-0.18, 0]    |Mean            |
#> |15 |one_to_x2  |        -0.08|0.04   |-0.08 [-0.17, 0]    |Mean            |
#> |16 |one_to_x1  |        -0.09|0.04   |-0.09 [-0.18, 0]    |Mean            |
#> |6  |x5_with_x5 |         0.05|0.01   |0.05 [0.04, 0.06]   |Residual        |
#> |7  |x4_with_x4 |         0.07|0.01   |0.07 [0.06, 0.09]   |Residual        |
#> |8  |x3_with_x3 |         0.11|0.01   |0.11 [0.09, 0.13]   |Residual        |
#> |9  |x2_with_x2 |         0.13|0.01   |0.13 [0.11, 0.15]   |Residual        |
#> |10 |x1_with_x1 |         0.21|0.02   |0.21 [0.17, 0.24]   |Residual        |
#> 
#> *Statistical Note*: Robust ML indices use sandwich-based scaling + Brosseau-Liard & Savalei (2012) adjustments to CFI/TLI/RMSEA (c_model = 1.005, c_null = 0.999).

We can also use imxRobustSE() from OpenMx on our fitted model, returning a vector of SEs:

m1Robust = imxRobustSE(m1)

Compare SEs

If robust SEs are larger than standard ML SEs, that often signals non-normality or other sandwich inflation.

stdSE = m1$output$standardErrors
compareSE = data.frame(
    Parameter = row.names(stdSE),
    Standard_ML_SE = round(as.numeric(stdSE), 4),
    Robust_SE = round(as.numeric(m1Robust), 4)
)
print(compareSE)
#>     Parameter Standard_ML_SE Robust_SE
#> 1     g_to_x5         0.0266    0.0249
#> 2     g_to_x4         0.0240    0.0244
#> 3     g_to_x3         0.0204    0.0215
#> 4     g_to_x2         0.0182    0.0185
#> 5     g_to_x1         0.0155    0.0153
#> 6  x5_with_x5         0.0037    0.0032
#> 7  x4_with_x4         0.0034    0.0036
#> 8  x3_with_x3         0.0031    0.0031
#> 9  x2_with_x2         0.0028    0.0026
#> 10 x1_with_x1         0.0028    0.0029
#> 11  one_to_x5         0.0366    0.0366
#> 12  one_to_x4         0.0326    0.0327
#> 13  one_to_x3         0.0273    0.0274
#> 14  one_to_x2         0.0241    0.0242
#> 15  one_to_x1         0.0199    0.0199

Step 3: Profile likelihood CIs (two parameters)

Profile CIs re-optimize the model many times, so keep the set small. Here we show two loadings:

m1CI = umxConfint(m1, parm = c("g_to_x1", "g_to_x2"), run = TRUE)
#>         lbound estimate ubound lbound Code ubound Code
#> g_to_x1  0.367    0.397  0.429           0           0
#> g_to_x2  0.469    0.503  0.541           0           0
umxSummary(m1CI, uncertainty = "CI", refModels = FALSE)
#> 
#> 
#> Table: Parameter loadings for model 'One_Factor_Model'
#> 
#> |   |name       | Estimate|CI                  |type            |
#> |:--|:----------|--------:|:-------------------|:---------------|
#> |1  |g_to_x5    |     0.80|0.8 [0.74, 0.85]    |Factor loading  |
#> |2  |g_to_x4    |     0.70|0.7 [0.66, 0.75]    |Factor loading  |
#> |3  |g_to_x3    |     0.58|0.58 [0.54, 0.62]   |Factor loading  |
#> |4  |g_to_x2    |     0.50|0.5 [0.47, 0.54]    |Factor loading  |
#> |5  |g_to_x1    |     0.40|0.4 [0.37, 0.43]    |Factor loading  |
#> |11 |g_with_g   |     1.00|1 [1, 1]            |Factor Variance |
#> |12 |one_to_x5  |    -0.08|-0.08 [-0.15, 0]    |Mean            |
#> |13 |one_to_x4  |    -0.06|-0.06 [-0.12, 0.01] |Mean            |
#> |14 |one_to_x3  |    -0.06|-0.06 [-0.11, 0]    |Mean            |
#> |15 |one_to_x2  |    -0.05|-0.05 [-0.09, 0]    |Mean            |
#> |16 |one_to_x1  |    -0.04|-0.04 [-0.08, 0]    |Mean            |
#> |6  |x5_with_x5 |     0.04|0.04 [0.03, 0.04]   |Residual        |
#> |7  |x4_with_x4 |     0.04|0.04 [0.03, 0.05]   |Residual        |
#> |8  |x3_with_x3 |     0.04|0.04 [0.03, 0.05]   |Residual        |
#> |9  |x2_with_x2 |     0.04|0.04 [0.03, 0.04]   |Residual        |
#> |10 |x1_with_x1 |     0.04|0.04 [0.04, 0.05]   |Residual        |
#>            lbound  estimate    ubound
#> g_to_x1 0.3674249 0.3967551 0.4286077
#> g_to_x2 0.4690729 0.5031577 0.5406042

Profile CIs are typically asymmetric and respect the likelihood surface and bounds.

4. Summary: When to Prefer Which?

Feature Standard ML SEs Robust SEs (Sandwich) Profile Likelihood CIs
Computational Cost Extremely low (one Hessian inversion) Low (gradient calculations) High (requires re-optimization)
Normality Assumption Strict Relaxed Strict (unless robust version used)
Handles Boundaries? No (can yield negative variances) No Yes (respects constraints)
Asymmetric CIs? No (always symmetric) No Yes (reflects log-likelihood shape)
Best For… Large samples, normal data Non-normal data, rapid estimation Key parameters, small samples, bounded parameters

Recommendations for Reporting

  1. Use Robust SEs for general parameters (such as minor residual variances or covariances) if your data shows signs of non-normality.
  2. Compute Profile Likelihood CIs (using umxConfint) for your key structural parameters (e.g., hypothesized paths, mediation effects) to ensure your scientific inferences are robust to asymmetry and boundary conditions.