Converting a Stata -sem- path model
Robert Long
Joined: 02/22/2013
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
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?
Log in or register to post comments
In reply to At first glance everything by mhunter
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
Log in or register to post comments
In reply to I think I may have by Robert Long
head=T
Yet another good reason to keep data labels in the data file, not separate :-)
Log in or register to post comments