covariance matrix of coefficient estimates
Posted on

Forums
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)
vcov
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"))
}
}
Log in or register to post comments
In other words
2 * solve(fit$output$hessian)
or if you want to wait for the next release it will bevcov(fit)
.Log in or register to post comments
Thanks!
Log in or register to post comments
October
Log in or register to post comments