You are here

Expected covariance matrix is non-positive-definite

6 posts / 0 new
Last post
Robert Long's picture
Offline
Joined: 02/22/2013 - 10:08
Expected covariance matrix is non-positive-definite

Hi all

I'm new to OpenMx. Actually I'm trying to convert some SEMs written in Stata into R for a module that I am helping to deliver, and for better or worse, we have chosen OpenMx as the R package to use.

It is a very simple path analysis. In Stata the code is just

    ssd init BMI SO RT SE
    ssd set correlations ///
    1 \ ///
    0.43 1 \ ///
    0.16 0.53 1 \ ///
    0.02 -0.18 -0.25 1 
 
ssd set sd 4.46 7.09 5.68 5.97
ssd set observations 202
 
sem ///
    ( SE <- RT ) ///
    ( RT <- SO ) ///
    ( SO <- BMI ) ///
    , standardized nocapslatent

Using OpenMX I have tried to replicate this with

require(OpenMx)
require(MBESS)
 
selVars <- c("BMI","SO","RT","SE")
 
cormat <- matrix(c(1,0.43,0.16,0.02,
        0.43,1,0.53,-0.18,
        0.16,0.53,1,-0.25,
        0.02,-0.18,-0.25, 1
        ), 
    nrow=4)
 
sds <- c(4.46, 7.09, 5.68, 5.97)
 
covmat <- round(cor2cov(cormat,sds),2)
 
rownames(covmat) <- selVars
colnames(covmat) <- selVars
 
m.tigg <- mxModel("mtigg",
    manifestVars= selVars,
    mxPath(
        from=c("BMI", "SO"), 
        arrows=1, 
        free=T, 
        values=1, 
        lbound=.01, 
        labels=c("p1",p2")
    ),
    mxPath(
        from=c("SO", "RT"), 
        arrows=1, 
        free=T, 
        values=1, 
        lbound=.01, 
        labels=c("p2","p3")
    ),
    mxPath(
        from=c("RT", "SE"), 
        arrows=1, 
        free=T, 
        values=1, 
        lbound=.01, 
        labels=c("p3","p4")
    ),
    mxData(
        observed=covmat, 
        type="cov", 
        numObs=202 
    ),
    type="RAM"
)
 
f.tigg <- mxRun(m.tigg)

However it complains that "Expected covariance matrix is non-positive-definite."

I would really appreciate any guidance about what I'm doing wrong as I have a lot of these models to convert, and very little time to do so....

FWIW, the output I get in Stata (which also agrees with MPlus) is

------------------------------------------------------------------------------
             |                 OIM
Standardized |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
Structural   |
  SE <-      |
          RT |       -.25   .0659595    -3.79   0.000    -.3792782   -.1207218
  -----------+----------------------------------------------------------------
  RT <-      |
          SO |        .53   .0504741    10.50   0.000     .4310726    .6289274
  -----------+----------------------------------------------------------------
  SO <-      |
         BMI |        .43   .0546349     7.87   0.000     .3229175    .5370825
-------------+----------------------------------------------------------------
Variance     |
        e.SE |      .9375   .0329797                      .8750389     1.00442
        e.RT |      .7191   .0535025                      .6215243    .8319945
        e.SO |      .8151   .0469861                      .7280208    .9125948
------------------------------------------------------------------------------
LR test of model vs. saturated: chi2(3)   =      4.11, Prob > chi2 = 0.2495

Ultimately, this is what I am trying to replicate.

Thanks a lot
Robert Long

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
Welcome to OpenMx!

Welcome to OpenMx!

I can attest to their being some "cost of entry", but there's some of that with any SEM program.

Here is the model I think you meant to specify.

require(OpenMx)
require(MBESS)
 
selVars <- c("BMI","SO","RT","SE")
 
cormat <- matrix(c(1,0.43,0.16,0.02,
        0.43,1,0.53,-0.18,
        0.16,0.53,1,-0.25,
        0.02,-0.18,-0.25, 1
        ), 
    nrow=4)
 
sds <- c(4.46, 7.09, 5.68, 5.97)
 
covmat <- round(cor2cov(cormat,sds),2)
 
rownames(covmat) <- selVars
colnames(covmat) <- selVars
 
m.tigg <- mxModel("mtigg",
    manifestVars= selVars,
    type="RAM",
    mxPath(
        from=selVars[1:3],
    to=selVars[2:4],
        arrows=1, 
        free=T, 
        values=1, 
        labels=c("p1","p2","p3")
    ),
    mxPath(
      from=selVars,
      arrows=2,
      values=.9,
      labels=paste("var", selVars, sep=""),
      free=TRUE
    ),
    mxData(
        observed=round(cov2cor(covmat), 3), 
        type="cov", 
        numObs=202 
    )
)
 
f.tigg <- mxRun(m.tigg)
summary(f.tigg)
Robert Long's picture
Offline
Joined: 02/22/2013 - 10:08
Dear mhunter Thank you so

Dear mhunter

Thank you so much for this ! The output agrees nicely with Stata/MPlus now. I do apreciate that there is going to be a learning curve with this - as with all packages. At this stage I am just confused about this part of the model:

from=selVars[1:3],
to=selVars[2:4],

I am a bit confused by this, as there is only one exogenous variable: the paths just go from 1->2->3->4

Thanks again
RL

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
In OpenMx, paths go FROM one

In OpenMx, paths go FROM one variable TO another. The "from" and the "to" arguments can be given vectors where the default behavior is to match elementwise. I could have specified several paths: (1) 1->2, (2) 2->3, and (3) 3->4. Instead, I specified one path:

from  to
1       2
2       3
3       4

The square brackets just select some of the elements from the vector selVars.

Hope this helps! Let us know if we can clear anything else up.

Robert Long's picture
Offline
Joined: 02/22/2013 - 10:08
Thanks again mhunter, that

Thanks again mhunter, that makes perfect sense !

When I looked at the examples in the tutorial I thought mxData was working with a covariance matrix, and I had not got far enought through the docs to see that a correlation matrix is also possible. However, I see you have used the correlation matrix - yet the type is still "cov", and changing from "cov" to "cor" doesn't seem to change anything. Am I missing something ?

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
I think it should only change

I think it should only change the degrees of freedom in this case.

Glad to help!