You are here

Converting a Stata -sem- path model

4 posts / 0 new
Last post
Robert Long's picture
Offline
Joined: 02/22/2013 - 10:08
Converting a Stata -sem- path model

Hi all

I am trying to convert a simple path model from Stata. In Stata I use

clear all
ssd init EE DR RO SE 
ssd set observations 321
 
ssd set sd 15.23 10.15 16.27 1.01 
 
ssd set correlations ///
 1 \ ///
 0.25 1 \ ///
 0.78 0.21 1 \ ///
 -0.47 -0.27 -0.35 1
 
 sem (SE -> EE) (EE -> DR) (EE -> RO) (DR -> RO), nocapslatent

This generates the following output

------------------------------------------------------------------------------
             |                 OIM
             |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
Structural   |
  EE <-      |
          SE |  -7.087228   .7428868    -9.54   0.000    -8.543259   -5.631196
       _cons |   59.18695   3.629161    16.31   0.000     52.07392    66.29997
  -----------+----------------------------------------------------------------
  DR <-      |
          EE |    .166612   .0360163     4.63   0.000     .0960213    .2372026
       _cons |   16.53305   1.063443    15.55   0.000     14.44874    18.61736
  -----------+----------------------------------------------------------------
  RO <-      |
          EE |   .8289902   .0385246    21.52   0.000     .7534834    .9044969
          DR |   .0256473   .0578058     0.44   0.657      -.08765    .1389446
       _cons |   51.46608   1.458225    35.29   0.000     48.60801    54.32415
-------------+----------------------------------------------------------------
Variance     |
        e.EE |   180.1515   14.22003                      154.3297    210.2938
        e.DR |   96.28271   7.599951                      82.48212    112.3924
        e.RO |   103.2753   8.151902                      88.47244    120.5549
------------------------------------------------------------------------------
LR test of model vs. saturated: chi2(2)   =     10.78, Prob > chi2 = 0.0046

My attempt to reproduce this in OpenMX is

vars <- c("SE","EE","DR","RO")
 
cormat <- matrix(c(1,0.25, 0.78, -0.47,
           0.25, 1, 0.21, -0.27,
           0.78, 0.21 , 1, -0.35,
            -0.47, -0.27, -0.35, 1
        ), nrow=4)
 
 
sds <- c(15.23, 10.15, 16.27, 1.01) 
 
rownames(cormat) <- colnames(cormat) <- vars
covmat <- round(cor2cov(cormat,sds),2)
 
m.cleary <- mxModel("mcleary",
    manifestVars = vars,
    type="RAM",
    mxPath(
        from=vars[c(1,2,3,2)], to=vars[c(2,3,4,4)], 
        values=1, 
        labels=c("p1","p2","p3","p4")
    ),
    mxPath(
        from=vars,        # “to=vars” is implicit 
        arrows=2,        # the default is 1
        values=.9,
        labels=paste("var", vars, sep="")
    ),
    mxData(
        observed=covmat, 
        type="cov",  # data is a covariance matrix
        numObs=321 
    )
)
 
f.cleary <- mxRun(m.cleary)
summary(f.cleary)

and this produced:

data:
$mcleary.data
$mcleary.data$cov
       SE     EE     DR    RO
SE 231.95  38.65 193.28 -7.23
EE  38.65 103.02  34.68 -2.77
DR 193.28  34.68 264.71 -5.75
RO  -7.23  -2.77  -5.75  1.02
 
 
free parameters:
   name matrix row col     Estimate    Std.Error Std.Estimate     Std.SE lbound ubound
1    p1      A  EE  SE   0.16663073  0.036072013    0.2500296 0.05412611              
2    p2      A  DR  EE   0.33663366  0.087610155    0.2100067 0.05465502              
3    p4      A  RO  EE  -0.02047884  0.005205191   -0.2058098 0.05231152              
4    p3      A  RO  DR  -0.01903892  0.003247224   -0.3067097 0.05231152              
5 varSE      S  SE  SE 231.94999676 18.337149808    1.0000000 0.07905648              
6 varEE      S  EE  EE  96.57972313  7.635160662    0.9374852 0.07411338              
7 varDR      S  DR  DR 253.03554502 20.004099673    0.9558972 0.07556987              
8 varRO      S  RO  RO   0.85379980  0.067498669    0.8370586 0.06617517              
 
observed statistics:  10 
estimated parameters:  8 
degrees of freedom:  2 
-2 log likelihood:  6205.556 
saturated -2 log likelihood:  5886.007 
number of observations:  321 
chi-square:  319.5493 
p:  4.080808e-70 

Could anyone advise where I am going wrong ?

Thanks
RL

mhunter's picture
Offline
Joined: 07/31/2009 - 15:26
At first glance everything

At first glance everything looks like it should be correct. The only thing I can see that might be odd is the "_cons" term that goes with all the coefficients in the Stata output you provided. I'm guessing it's a kind of constant/intercept term. This was not in the output for the previous model you asked about. http://openmx.psyc.virginia.edu/thread/1934.

The "_cons" term makes me think the Stata output was generated off raw data, and/or observed means. If so, then the corresponding OpenMx model would need the same information. Could that be the case?

Robert Long's picture
Offline
Joined: 02/22/2013 - 10:08
I think I may have

I think I may have inadvertently posted the output for the model which ran from raw data instead of the covariance matrix, although the path coefficients and variances are the same, so that isn't part of the problem...

However, I have now solved it ! I made the mistake of changing the order of the variables (and standard deviations), but not changing the covariance matrix !! So the correct model is

vars <- c("EE","DR","RO","SE")
 
cormat <- matrix(c(1,0.25, 0.78, -0.47,
           0.25, 1, 0.21, -0.27,
           0.78, 0.21 , 1, -0.35,
          -0.47, -0.27, -0.35, 1
        ), nrow=4)
 
 
sds <- c(15.23, 10.15, 16.27, 1.01)
 
rownames(cormat) <- colnames(cormat) <- vars
covmat <- round(cor2cov(cormat,sds),2)
 
m.cleary <- mxModel("mcleary",
    manifestVars = vars,
    type="RAM",
    mxPath(
        from=vars[c(4,1,1,2)], to=vars[c(1,2,3,3)], 
        values=1, 
        labels=c("p1","p2","p3","p4")
    ),
    mxPath(
        from=vars,        
        arrows=2,        
        values=.9,
        labels=paste("var", vars, sep="")
    ),
 
    mxData(
        observed=covmat, 
        type="cov",  
        numObs=321 
    )
)
 
f.cleary <- mxRun(m.cleary)
summary(f.cleary)

And this now agrees with Stata.

Thanks again for your help
RL

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
head=T

Yet another good reason to keep data labels in the data file, not separate :-)