Implementation of Oberski, D. L., & Satorra, A.. (2013). Measurement error models with uncertainty about the error variance. Structural equation modeling, 20, 409-428. doi:10.1080/10705511.2013.797820. See also daob.nl/publications

oberski-satorra-lavaan.R

# This code is a lavaan implementation of the standard error correction found in 
#
# Oberski, D. L., & Satorra, A. (2013).
#   Measurement error models with uncertainty about the error variance. 
#   Structural equation modeling, 20, 409-428. 
#   DOI:10.1080/10705511.2013.797820 
#
# Author: Daniel Oberski
# Date 27 november 2018
# License: MIT (https://opensource.org/licenses/MIT)

# Core functions that implement the correction

# This function is intended only for calling from vcov_corrected:
get_additional_variance <- function(object, fixed, Sigma_fixed) {

    strict.exo <- FALSE

    if(object@[email protected]) {
        strict.exo <- TRUE
    }

    FIT <- lavaan:::lav_object_extended(object, add = fixed, remove.duplicated = FALSE)
    LIST <- FIT@ParTable

    free.idx <- LIST$free[ LIST$free > 0L & LIST$user != 10L ]
    fixed.idx <- LIST$free[ LIST$free > 0L & LIST$user == 10L ]


    Delta <- lavaan:::computeDelta(FIT@Model)[[1]] # Only single-group for now

    Delta_free <- Delta[, free.idx, drop=FALSE]
    Delta_fixed <- Delta[, fixed.idx, drop=FALSE]

    V <- lavTech(fit, "WLS.V")[[1]]

    H <- t(Delta_free) %*% V %*% Delta_fixed

    # Ji <- t(Delta_free) %*% V %*% Delta_free
    Ji <- lavTech(fit, "information.expected")

    Ji %*% H %*% Sigma_fixed %*% t(H) %*% Ji
}

# This function is intended for calling by the user:
# Input: 
#  - lavaan object, 
#  - vector of lavaan names of fixed parameters, 
#  - covariance matrix of fixed parameters
# Output:
#  estimated variance-covariance matrix of free parameters after correcting for uncertainty in fixed parameters
vcov_corrected <- function(object, fixed, Sigma_fixed) {

  stopifnot(length(fixed) == ncol(Sigma_fixed)) # Same no. parameters
  stopifnot(ncol(Sigma_fixed) == nrow(Sigma_fixed)) # Square
  stopifnot(all(eigen(Sigma_fixed, only.values = TRUE)$values > 1e-4)) #Positive definite
  stopifnot(all(t(Sigma_fixed) == Sigma_fixed)) # Symmetric

  C_psi <- get_additional_variance(object, fixed = fixed, Sigma_fixed = Sigma_fixed)
  vcov(object) + C_psi
}


## An example application of the se correction

library(lavaan)

# Errors-in-variables regression model
eiv_regression <- '
  visual  =~ 1*x1
  textual =~ 1*x4

  textual ~ visual

  x1~~1*x1
  x4~~0.3*x4
'
# First, fit the model as usual in lavaan
fit <- sem(eiv_regression, data=HolzingerSwineford1939)
summary(fit)

# Specify variance-covariance matrix of parameters fixed to estimates:
Sigma_fixed <- diag(c(0.1, 1))

# Obtain corrected covariance matrix of free parameters
vcov_corrected(fit, fixed = c("x1~~x1", "x4~~x4"), Sigma_fixed = Sigma_fixed)

vcov(fit)