# multilevel path/structural equation models in OpenMx

30 posts / 0 new
Offline
Joined: 05/08/2017 - 11:45
multilevel path/structural equation models in OpenMx
AttachmentSize
1.26 MB
158.41 KB

Hi,

Is it possible to estimate multilevel path/structural equation models with OpenMx?
If it is possible, is there an example code for multilevel mediation (path) model?

Attached files are the example of multilevel path model using Mplus.
Source: Heck, R. H., & Thomas, S. L. (2015). An introduction to multilevel modeling techniques: MLM and SEM approaches using Mplus. Routledge.

Hope I can do the same analysis using OpenMx.

Soyoung.

Offline
Joined: 03/01/2013 - 11:03
Yep

It is possible to estimate multilevel structural equations in OpenMx. This paper by Pritikin et al. discusses many of the details and provides some example scripts.

Many-Level Multilevel Structural Equation Modeling: An Efficient Evaluation Strategy

Offline
Joined: 12/12/2012 - 12:15
voila

Here's a similar model (see attachment). You probably need to add/remove a few paths to make it identical to what you posted. I find it tedious to read Mplus model diagrams or tables.

File attachments:
Offline
Joined: 05/08/2017 - 11:45
fabulous

Thank you all.

Two more questions.

[1]
I revised the code, and i though I made the same model.
BUT I found these results are not identical in terms of estimates.
Why did these have different results?
Which Estimation method the OpenMx use?
Mplus uses Robust-ML.

[2]
I want to estimate indirect effects in Multilevel SEM.
How can I set a product term (indirect effect parameter) in OpenMx?

File attachments:
Offline
Joined: 12/12/2012 - 12:15
questions
1. See imxRobustSE

2. Create an mxAlgebra and use mxCI or mxSE to request confidence intervals on it. See attached code.

File attachments:
Offline
Joined: 05/08/2017 - 11:45
Thank you

Thank you so much for your kind response.

All the best,
Soyoung.

Offline
Joined: 03/01/2013 - 14:09
Agreement

Hi Soyoung

Could you confirm that you were able to get the same estimates and goodness-of-fit statistics with the multilevel model in OpenMx when compared to Mplus?

Thanks
Mike

Offline
Joined: 05/08/2017 - 11:45
could not get the identical results

Hi, Mike.
I still could not get the identical results.
I am curious why these are different.

First, in terms of model fit statistics, two results were almost same.

 Mplus OpenMx RMSEA 0 0 CFI 1 NA TLI 1.001 1

BUT I don't think this example is suitable here, because these fit values are too perfect. It should be tested in the other data.

Second, in terms of estimates, I could not get the same estimates.
I'm trying to figure out why there are different results.

File attachments:
Offline
Joined: 01/24/2014 - 12:15
imxRobustSE()

I am pretty sure that imxRobustSE() does not correctly handle multilevel data; the matter of whether or not it does is currently a matter of disagreement among the OpenMx development team. I suggest you not rely too much on the robust SEs you get out of OpenMx.

Soyoung, could you post the OpenMx script you're currently working from?

Offline
Joined: 05/08/2017 - 11:45
sure

Hi Rob

Following code is based on what Dr.Pritikin provided.
You can find the data in post #1 .

library(OpenMx)
colnames(flatData) <- c("orgid",'female','white',
'satpay','morale','org1','org2','benefit','cond',
'resour','zproduct','lev1wt','lev2wt')

level1 <- flatData[,c('benefit','cond','female','white','orgid')]
level2 <- flatData[!duplicated(flatData$orgid),c('orgid', 'org1','org2','zproduct')] orgModel <- mxModel( "org", type="RAM", mxData(level2, 'raw', primaryKey = "orgid"), manifestVars = c('zproduct'), latentVars = c('Benefits','Conditions', 'org1','org2'), mxPath(c('zproduct','Benefits','Conditions'), arrows=2, values=1), mxPath('Benefits','Conditions', arrows=2), mxPath('one', 'org1', labels='data.org1', free=FALSE), mxPath('one', 'org2', labels='data.org2', free=FALSE), mxPath('org1', c('Benefits', 'Conditions')), mxPath('org2', 'zproduct', labels="a"), mxPath('one', 'zproduct'), mxPath('zproduct', c('Benefits','Conditions'), labels=paste0('b',1:2)), mxMatrix(nrow=2, ncol=1, free=TRUE, labels=paste0('b',1:2), name='B'), mxAlgebra(B * a, name="indirect") ) empModel <- mxModel( "emp", type="RAM", orgModel, mxData(level1, 'raw'), manifestVars = c('benefit','cond'), latentVars = c('female','white'), mxPath('one','female', labels='data.female', free=FALSE), mxPath('one','white', labels='data.white', free=FALSE), mxPath(c('female','white'), c('benefit','cond'), connect = "all.bivariate"), mxPath(c('benefit','cond'), arrows=2, connect = "unique.pairs", values=c(1,0,1)), mxPath('org.Benefits', 'benefit', values=1, free=FALSE, joinKey = "orgid"), mxPath('org.Conditions', 'cond', values=1, free=FALSE, joinKey = "orgid"), mxPath('one', c('benefit','cond')), mxCI("org.indirect") ) ch4M3 <- mxRun(empModel, intervals=TRUE) summary(ch4M3) imxRobustSE(ch4M3) Offline Joined: 01/24/2014 - 12:15 Under x86_64 Linux/GNU, I Under x86_64 Linux/GNU, I substantively reproduce the OpenMx results you report with all 3 gradient-descent optimizers (script and output are attached). Could you post the full output from MPlus? In particular, I'm curious if the -2logL matches. File attachments: Offline Joined: 05/08/2017 - 11:45 full output Yes. Attachment is full output. Mplus reports LL -40521.828. (-2) * (-40521.828) = 81043.656 Thank you. Soyoung. File attachments: Offline Joined: 03/01/2013 - 14:09 ML vs MLR? Hi Soyoung, I notice that the estimator in MPlus is MLR, whereas in OpenMx it is FIML. Could you try it with ML (FIML) in Mplus to see if that changes the results? Cheers Mike Offline Joined: 05/08/2017 - 11:45 different Hi Mike. Attachment is the result of FIML, but this is almost same as the result of MLR of Mplus, and the estimates does not match the result of OpenMx. I wonder why the two results are different...Where does the difference come from? Is there a way to check local-maxima in OpenMx, by any chance? Thank you. Soyoung. File attachments: Offline Joined: 05/24/2012 - 00:35 fit value One way to debug this is to load the parameter vectors into the OpenMx model using omxSetParameters and then evaluate the fit at that point using mxRun(mxModel(model, mxComputeOnce('fitfunction','fit'))). If the Mplus fit matches at a specific parameter vector then the models are likely the same. If the fit doesn't match then the models are different. Offline Joined: 01/24/2014 - 12:15 already tried it One way to debug this is to load the parameter vectors into the OpenMx model using omxSetParameters and then evaluate the fit at that point using mxRun(mxModel(model, mxComputeOnce('fitfunction','fit'))). I already tried that (see my post, presently #17, below). The fit values don't match. Offline Joined: 03/01/2013 - 14:09 Local maxima possibility Hi Soyoung OpenMx does have facilities for testing for local minima, such as mxTryHard(). However, in this case it seems to me that Mplus is the software that might need it, because the fit function from OpenMx's estimation is indicating a much higher likelihood (41000.78 vs 40522.071) than is Mplus. You could perhaps try putting the OpenMx estimates into Mplus and fixing all the parameters and seeing what it thinks the answer is with those values. If the answers still don't agree, it indicates that the models are not the same, and we should try to figure out why they aren't. Cheers Mike Offline Joined: 05/24/2012 - 00:35 robust Good point Rob. I am not sure what Robust-ML means for multilevel models. Offline Joined: 01/24/2014 - 12:15 It appears that imxRobustSE() It appears that imxRobustSE() does not calculate robust SEs in a manner theoretically correct for multilevel data (e.g., Freedman 2006). It also appears that what OpenMx and MPlus are doing in this case are different enough that their results shouldn't be compared to each other. See the attached script and output. When OpenMx calculates the fit at MPlus' solution, OpenMx returns a different fit from MPlus. Further, I can start OpenMx at MPlus' solution, and OpenMx still returns the same solution it has before (with a fit value of 82001.55). I don't know OpenMx's multilevel feature very well, and I don't know MPlus at all, so I'm not sure what might explain the difference (or even if Soyoung's syntax is even fitting equivalent models in the two software packages). Perhaps a difference in how the two programs internally handle zero-variance exogenous regressors? Does MPlus drop normalizing constants from its loglikelihood? File attachments: Offline Joined: 05/08/2017 - 11:45 I have been told that the I have been told that the Mplus software excludes any variable with no variance from analysis. and Mplus does not drop normalizing constants. In my eyes, OpenMx and Mplus are different in the manner of providing results. In OpenMx, when there are level-1 manifest variables causes level-2 latent variables, intercepts are showed at the lower level results. In Mplus, when there are level-1 manifest variables causes level-2 latent variables, intercepts are showed at the between level (upper level) results. When I look at the mathematical formulas of the model, I feel easy to understand to treat the intercept at the upper level. (probably because I got used to it) Will the outcome differ depending on the level of handling the intercept? Thank you again, Soyoung. Offline Joined: 05/24/2012 - 00:35 means > Will the outcome differ depending on the level of handling the intercept? No and you can prove it by trying it both ways. Offline Joined: 05/24/2012 - 00:35 bug OK, there is an OpenMx bug. Try the model again with, empModel$expectation$.useSufficientSets <- FALSE Offline Joined: 05/24/2012 - 00:35 looks better now Summary of emp free parameters: name matrix row col Estimate Std.Error A 1 emp.A[1,3] A benefit female 0.69700745 0.02567046 2 emp.A[2,3] A cond female 0.70677248 0.02514620 3 emp.A[1,4] A benefit white 0.22406279 0.02694068 4 emp.A[2,4] A cond white 0.24780137 0.02630581 5 emp.S[1,1] S benefit benefit 2.02929959 0.02590178 6 emp.S[1,2] S benefit cond 1.34397122 0.02165825 7 emp.S[2,2] S cond cond 1.94820623 0.02485886 8 emp.M[1,1] M 1 benefit 4.41229957 0.04570310 9 emp.M[1,2] M 1 cond 4.67685717 0.04025743 10 b1 org.A Benefits zproduct 0.31132159 0.03144856 11 b2 org.A Conditions zproduct 0.29100892 0.02721743 12 org.A[2,4] org.A Benefits org1 0.28727342 0.07941363 13 org.A[3,4] org.A Conditions org1 0.31525677 0.06869425 14 a org.A zproduct org2 -0.45125514 0.20589852 15 org.A[2,5] org.A Benefits org2 0.03564603 0.08489777 16 org.A[3,5] org.A Conditions org2 0.02533898 0.07302397 17 org.S[1,1] org.S zproduct zproduct 1.07059845 0.11934389 ! 18 org.S[2,2] org.S Benefits Benefits 0.13261944 0.01900351 19 org.S[2,3] org.S Benefits Conditions 0.09517323 0.01495373 20 org.S[3,3] org.S Conditions Conditions 0.09256918 0.01401026 21 org.M[1,1] org.M 1 zproduct 0.08692317 0.08979952 Model Statistics: | Parameters | Degrees of Freedom | Fit (-2lnL units) Model: 21 25029 81043.97 Saturated: NA NA NA Independence: NA NA NA Number of observations/statistics: 12605/25050 Offline Joined: 05/08/2017 - 11:45 Estimates are almost same/Where is the result of indirect effect Hi Joshua I am so pleased to see these estimates are almost same as the results of Mplus :) Where are the results of indirect effects such as a*b1 and a*b2? I am curious whether they are also same as the results of Mplus. Mplus (FIML) Loglikelihood H0 Value -40521.984 H1 Value -40521.828 Thank you so much, Soyoung. Offline Joined: 05/08/2017 - 11:45 Please never mind. Please never mind. I am going to try! :) Offline Joined: 05/08/2017 - 11:45 almost same Hello, l'd like to confirm that I can get the same estimates. and I also want to ask a question if I can ignore small difference of lbound/ubound. When I ran the following code, library(OpenMx) flatData <- read.table("ch4mv.dat") colnames(flatData) <- c("orgid",'female','white', 'satpay','morale','org1','org2','benefit','cond', 'resour','zproduct','lev1wt','lev2wt') level1 <- flatData[,c('benefit','cond','female','white','orgid')] level2 <- flatData[!duplicated(flatData$orgid),c('orgid', 'org1','org2','zproduct')]

orgModel <- mxModel(
"org", type="RAM",
mxData(level2, 'raw', primaryKey = "orgid"),
manifestVars = c('zproduct'),
latentVars = c('Benefits','Conditions', 'org1','org2'),
mxPath(c('zproduct','Benefits','Conditions'), arrows=2, values=1),
mxPath('Benefits','Conditions', arrows=2),
mxPath('one', 'org1', labels='data.org1', free=FALSE),
mxPath('one', 'org2', labels='data.org2', free=FALSE),
mxPath('org1', c('Benefits', 'Conditions')),
mxPath('org2', 'zproduct', labels="a"),
mxPath('one', 'zproduct'),
mxPath('zproduct', c('Benefits','Conditions'), labels=paste0('b',1:2)),
mxMatrix(nrow=2, ncol=1, free=TRUE, labels=paste0('b',1:2), name='B'),
mxAlgebra(B * a, name="indirect")
)

empModel <- mxModel(
"emp", type="RAM", orgModel,
mxData(level1, 'raw'),
manifestVars = c('benefit','cond'),
latentVars = c('female','white'),
mxPath('one','female', labels='data.female', free=FALSE),
mxPath('one','white', labels='data.white', free=FALSE),
mxPath(c('female','white'), c('benefit','cond'), connect = "all.bivariate"),
mxPath(c('benefit','cond'), arrows=2, connect = "unique.pairs", values=c(1,0,1)),
mxPath('org.Benefits', 'benefit', values=1, free=FALSE, joinKey = "orgid"),
mxPath('org.Conditions', 'cond', values=1, free=FALSE, joinKey = "orgid"),
mxPath('one', c('benefit','cond')),
mxCI("org.indirect")
)
empModel$expectation$.useSufficientSets <- FALSE
ch4M3 <- mxRun(empModel, intervals=TRUE)
summary(ch4M3)

I got the following result.

Summary of emp

free parameters:
name matrix        row        col    Estimate  Std.Error A
1  emp.A[1,3]      A    benefit     female  0.69705464 0.02567670
2  emp.A[2,3]      A       cond     female  0.70680474 0.02515014
3  emp.A[1,4]      A    benefit      white  0.22419999 0.02695238
4  emp.A[2,4]      A       cond      white  0.24789392 0.02630992
5  emp.S[1,1]      S    benefit    benefit  2.02928474 0.02590308
6  emp.S[1,2]      S    benefit       cond  1.34396875 0.02165891
7  emp.S[2,2]      S       cond       cond  1.94820919 0.02486024
8  emp.M[1,1]      M          1    benefit  4.42133171 0.04035614
9  emp.M[1,2]      M          1       cond  4.68326778 0.03568958
10         b1  org.A   Benefits   zproduct  0.30887125 0.03092515
11         b2  org.A Conditions   zproduct  0.28927449 0.02673664
12 org.A[2,4]  org.A   Benefits       org1  0.27821530 0.07658434
13 org.A[3,4]  org.A Conditions       org1  0.30882772 0.06618441
14          a  org.A   zproduct       org2 -0.45125100 0.21303763
15 org.S[1,1]  org.S   zproduct   zproduct  1.07059364 0.11968246
16 org.S[2,2]  org.S   Benefits   Benefits  0.13288795 0.01902547
17 org.S[2,3]  org.S   Benefits Conditions  0.09533848 0.01496872
18 org.S[3,3]  org.S Conditions Conditions  0.09267003 0.01402068
19 org.M[1,1]  org.M          1   zproduct  0.08692014 0.09164268

confidence intervals:
lbound   estimate      ubound note
org.indirect[1,1] -0.2768346 -0.1393785 -0.01158206
org.indirect[2,1] -0.2580154 -0.1305354 -0.01085831

Model Statistics:
|  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
Model:             19                  25031              81044.14
Saturated:             NA                     NA                    NA
Independence:             NA                     NA                    NA
Number of observations/statistics: 12605/25050

Information Criteria:
|  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       30982.14               81082.14                       NA
BIC:     -155294.78               81223.54                 81163.16
CFI: NA
TLI: 1   (also known as NNFI)
RMSEA:  0  [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2017-06-11 16:38:45
Wall clock time (HH:MM:SS.hh): 00:02:07.38
optimizer:  CSOLNP
OpenMx version number: 2.7.11
Need help?  See help(mxSummary) 

|

## Comparing the results of Mplus(Estimator=ML),

(attachment is the full output of Mplus).
Fit valus are same,
path coefficents are same,
BUT, lower bound / upper bound of indirect effects are not perfectly identical.

followings are from OpenMx,

confidence intervals:
lbound   estimate      ubound
org.indirect[1,1] -0.2768346 -0.1393785 -0.01158206
org.indirect[2,1] -0.2580154 -0.1305354 -0.01085831     

and followings are from Mplus.

CONFIDENCE INTERVALS OF MODEL RESULTS
Lower 2.5%      Estimate       Upper 2.5%
IND1       -0.269          -0.139          -0.010
IND2       -0.252          -0.131          -0.009


Can I ignore this difference?

Thank you,
Soyoung.

File attachments:
Offline
Joined: 05/24/2012 - 00:35
CI details

Use summary(model, verbose=TRUE) to get detailed information about the CIs. I'm not sure if Mplus makes this information available for comparison, but that's not a fault of OpenMx.

Offline
Joined: 05/08/2017 - 11:45
this is another example (2-1-1 model)

Hi,
I would like to ask about the other example of Multilevel SEM.

Now, I am trying 2-1-1 mediation model (using the attachment ex9.5.dat)
I checked the values of -2LnL are almost same; 4303.159 (OpenMx); 4303.158 (Mplus),
but some estimates are not identical.

If you look at the diagram file (attachment), a path and a intercept are different between two programs,
others are almost identical.

Could you please check if there is a way to get the same result?

I want to make this 2-1-1 model based on the article:
Preacher, K. J., Zhang, Z., & Zyphur, M. J. (2011). Alternative methods for assessing mediation in multilevel data: The advantages of multilevel SEM. Structural Equation Modeling, 18, 161-182.

my 2-1-1 model is based on E. 2-1-1 model (MSEM) in the Appendix above.

Soyoung.

P.S. following code is what I have run in OpenMx.

library(OpenMx)
colnames(exdata) <- c('Y', 'M', 'V1', 'V2', 'X', 'CLUS')
Data.L1 <- exdata[,c('CLUS','M','Y')]
Data.L2 <- exdata[!duplicated(exdata$CLUS),c('CLUS', 'X')] upperModel <- mxModel( "upper", type="RAM", mxData(Data.L2, 'raw', primaryKey = "CLUS"), manifestVars = c('X'), latentVars = c('Mb','Yb'), mxPath(c('X','Mb'), 'Yb', values=1), mxPath('X', 'Mb', values=1), mxPath('one', 'X'), # means mxPath(c('X','Mb','Yb'), arrows=2, values=1) ) lowerModel <- mxModel( "lower", type="RAM", upperModel, mxData(Data.L1, 'raw'), manifestVars = c('M', 'Y'), mxPath('M', 'Y', values=1), mxPath('one', c('M', 'Y')), # means mxPath(c('M', 'Y'), arrows=2), # variances mxPath('upper.Mb', 'M', values=1, free=FALSE, joinKey = "CLUS"), mxPath('upper.Yb', 'Y', values=1, free=FALSE, joinKey = "CLUS") ) lowerModel$expectation\$.useSufficientSets <- FALSE
model211 <- mxRun(lowerModel, intervals=TRUE)
summary(model211)

and i got the following result.

free parameters:
name  matrix row col    Estimate  Std.Error A
1  lower.A[2,1]       A   Y   M  0.33571331 0.02560867
2  lower.S[1,1]       S   M   M  5.39126176 0.37565626
3  lower.S[2,2]       S   Y   Y  1.45223197 0.10133471
4  lower.M[1,1]       M   1   M  0.70880568 0.20585133
5  lower.M[1,2]       M   1   Y -0.25604332 0.09926723
6  upper.A[2,1] upper.A  Mb   X  1.59487423 0.22672706
7  upper.A[3,1] upper.A  Yb   X  0.82937245 0.16255263
8  upper.A[3,2] upper.A  Yb  Mb -0.08507563 0.08169814
9  upper.S[1,1] upper.S   X   X  0.81028050 0.12078958
10 upper.S[2,2] upper.S  Mb  Mb  2.53495826 0.56401216
11 upper.S[3,3] upper.S  Yb  Yb  0.50668743 0.12456607
12 upper.M[1,1] upper.M   1   X -0.13019981 0.09488538
Model Statistics:
|  Parameters  |  Degrees of Freedom  |  Fit (-2lnL units)
Model:             12                   1078              4303.159
Saturated:             NA                     NA                    NA
Independence:             NA                     NA                    NA
Number of observations/statistics: 590/1090
Information Criteria:
|  df Penalty  |  Parameters Penalty  |  Sample-Size Adjusted
AIC:       2147.159               4327.159                       NA
BIC:      -2574.613               4379.720                 4341.624
CFI: NA
TLI: 1   (also known as NNFI)
RMSEA:  0  [95% CI (NA, NA)]
Prob(RMSEA <= 0.05): NA
To get additional fit indices, see help(mxRefModels)
timestamp: 2017-06-12 15:34:32
Wall clock time (HH:MM:SS.hh): 00:00:01.09
optimizer:  CSOLNP
OpenMx version number: 2.7.11
Need help?  See help(mxSummary)
File attachments:
Offline
Joined: 03/01/2013 - 14:09
Thank you for clarity

First, thank you for providing an example that runs nicely, and for clearly documenting the differences! This helps tremendously. Second, I have been able to reproduce the OpenMx side of the analysis and recover your estimates. The Mplus output is a bit brief and it doesn't seem to contain expected means or covariances - though these would seem likely to equal those of OpenMx because the model fit is the same. Sometimes expectations for the model implied means differ, and this may be the case here. I note the warning in the Mplus output:

*** WARNING in MODEL command
In the MODEL command, the following variable is a y-variable on the BETWEEN
level and an x-variable on the WITHIN level. This variable will be treated
as a y-variable on both levels: M

and I'm wondering if some of the components of the expected mean are omitted because of this 'treatment.' The difference between the two expected means of Y under OpenMx vs Mplus is .251*1.595*-.13-(-.085*1.595*-.13) = -0.0696696, which is the same as the expected mean 'due' to the path from X to Mb to M to Y, .336*(1.595*-.13) = -0.0696696. So it seems that the 'treated as y-variable' translates into 'the expected mean of Y lacks the effect of the mean of X transmitted through Mb to M to Y, but includes the effect of the mean of X transmitted through Yb.' Just how reasonable this is I will defer to others - maybe the Muthens will chime in here ;).

Hope this helps!

Offline
Joined: 05/08/2017 - 11:45
WARNING in Mplus is not error

Hi, Mike.
I have emailed the Mplus support team, and I have got the mail that the WARNING in Mplus is not error message. I think we can ignore the warning.

I have one more question.
I would like to ask about the way to handle level-1 variables.
Does OpenMx do latent variable decomposition?
If I want to see the latent variable decomposition results using OpenMx, what should I do?

If you look at the diagram file (attachment) variable M and Y are decomposed into two parts.
Mplus does a latent variable decomposition into a Within and Between part of M and Y.

There is some explanation about latent variable decomposition: