You are here

covariance matrix of coefficient estimates

5 posts / 0 new
Last post
dtofighi's picture
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)
jpritikin's picture
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"))
  }
}
mhunter's picture
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).

dtofighi's picture
Offline
Joined: 10/12/2009 - 14:55
Thanks!

Thanks for your fast response. Truly appreciated. Any time table for the next release?

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
October

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