# covariance matrix of coefficient estimates

5 posts / 0 new Offline
Joined: 10/12/2009 - 14:55
covariance matrix of coefficient estimates

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) Offline
Joined: 05/24/2012 - 00:35
vcov

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"))
}
} Offline
Joined: 07/31/2009 - 15:26
In other words

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). Offline
Joined: 10/12/2009 - 14:55
Thanks!

Thanks for your fast response. Truly appreciated. Any time table for the next release? Offline
Joined: 07/31/2009 - 15:26
October

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