covariance matrix of coefficient estimates

Posted on
Picture of user. dtofighi Joined: 10/12/2009
Hello,

I am wondering how to get covariance matrix of parameter estimates (e.g., inverse of Fisher's information matrix) in an SEM with OpenMx. I have included a sample code for a path model below. For example, I would like to get covariances between estimates of b1, b2, and cp that include in the covariance matrix.


library(OpenMx)
set.seed(1234)
n <- 100
x <- rnorm(n)
m <- 0.1 * x + rnorm(n)
y <- 0.3 * m + rnorm(n)
df <- data.frame(x = x, y = y, m = m)
manifests = c("x", "m", "y")
myModel <- mxModel(
  "Path Model",
  type = "RAM",
  manifestVars = manifests,
  mxPath(
    from = 'one',
    to = manifests,
    free = FALSE,
    values = 0
  ),
  mxPath(
    from = "x",
    to = "m",
    arrows = 1,
    free = TRUE,
    values = 0,
    labels = "b1"
  ),
  mxPath(
    from = c("x", "m"),
    to = "y",
    arrows = 1,
    free = TRUE,
    values = 0,
    labels = c("cp", "b2")
  ),
  mxPath(
    from = manifests,
    arrows = 2,
    free = TRUE,
    values = 1,
    labels = c("s2x", "s2m", "s2y")
  ),
  mxData(df, type = "raw")
)
fit <- mxRun(myModel)
summary(fit)
Replied on Thu, 09/21/2017 - 19:03
Picture of user. jpritikin Joined: 05/23/2012

Here's how to implement vcov. This will be included in the next release,


vcov.MxModel <- function(object, ...) {
if (!object@.wasRun) stop("This model has not been run yet. Tip: Use\n model = mxRun(model)\nto estimate a model.")
assertModelFreshlyRun(object)
fu <- object$output$fitUnits
if (fu %in% c("-2lnL", "r'Wr")) {
got <- NULL
if (!is.null(object$output[['ihessian']])) {
got <- 2 * object$output[['ihessian']]
} else if (!is.null(object$output[['hessian']])) {
got <- 2 * solve(object$output$hessian)
}
if (is.null(got)) {
stop(paste("Parameter variance covariance matrix is not available.",
"Turn on with mxOption(model, 'Calculate Hessian', 'Yes')", sep="\n"))
} else {
if (mxOption(object, "Calculate Hessian")== "No") {
warning(paste("The 'Calculate Hessian' option is disabled. This may result in a poor accuracy vcov matrix.",
"Turn on with mxOption(model, 'Calculate Hessian', 'Yes')", sep="\n"))
}
got
}
} else {
stop(paste("Don't know how to extract vcov from object",
omxQuotes(object$name), "fit in", omxQuotes(fu), "units"))
}
}

Replied on Thu, 09/21/2017 - 21:24
Picture of user. mhunter Joined: 07/31/2009

So if you need it now you can do 2 * solve(fit$output$hessian) or if you want to wait for the next release it will be vcov(fit).
Replied on Fri, 09/22/2017 - 12:19
Picture of user. dtofighi Joined: 10/12/2009

Thanks for your fast response. Truly appreciated. Any time table for the next release?
Replied on Fri, 09/22/2017 - 14:10
Picture of user. mhunter Joined: 07/31/2009

Internally we've discussed mid-October, so maybe 2-3 weeks.