You are here

Comparing sem() and mxRun() output

8 posts / 0 new
Last post
Matifou's picture
Offline
Joined: 12/08/2009 - 09:29
Comparing sem() and mxRun() output
AttachmentSize
Binary Data from sem to openmx.R3.29 KB

Hi

Are the algorithms different in the sem package and OpenMx package for solving a sem model? Or the starting values? The question is actually whether one is expected to get the same results from the different packages? (as an aside: how can I obtain the value of the evaluated objective function? after runnnig mxRun?)

I'm asking as I tried to replicate example "Wheaton et al. alienation data " (see ?sem), I think I succeeded in having the same model but coef differ slighty or much...

run attached example and compare:
summary(Run)
name matrix row col Estimate Std.Error
1 A SEI SES 5.2730210 4.0774084
2 A Alienation67 SES -0.6321773 0.5495341
3 A Alienation71 SES -0.2325439 0.5333972
4 A Alienation71 Alienation67 0.6105953 0.4505951
5 the1 S Anomia67 Anomia67 4.1128138 2.3309929
6 the2 S Powerless67 Powerless67 3.1626654 1.6110255
7 S Anomia67 Anomia71 1.5514337 2.2886042
8 S Education Education 2.8734482 4.7788626
9 S SEI SEI 262.9797477 173.9318663
10 S SES SES 6.7365517 6.1315362
11 S Alienation67 Alienation67 5.7430329 3.9999990
12 S Alienation71 Alienation71 4.5054334 3.2001184

summary(sem.wh.1)
Parameter Estimates
Estimate Std Error z value Pr(>|z|)
lamb 5.36880 0.433982 12.3710 0.0000e+00 SEI <--- SES
gam1 -0.62994 0.056128 -11.2233 0.0000e+00 Alienation67 <--- SES
beta 0.59312 0.046820 12.6680 0.0000e+00 Alienation71 <--- Alienation67
gam2 -0.24086 0.055202 -4.3632 1.2817e-05 Alienation71 <--- SES
the1 3.60787 0.200589 17.9864 0.0000e+00 Anomia67 <--> Anomia67
the2 3.59494 0.165234 21.7567 0.0000e+00 Powerless67 <--> Powerless67
the3 2.99366 0.498972 5.9996 1.9774e-09 Education <--> Education
the4 259.57583 18.321121 14.1681 0.0000e+00 SEI <--> SEI
the5 0.90579 0.121710 7.4422 9.9032e-14 Anomia71 <--> Anomia67
psi1 5.67050 0.422906 13.4084 0.0000e+00 Alienation67 <--> Alienation67
psi2 4.51481 0.334993 13.4773 0.0000e+00 Alienation71 <--> Alienation71
phi 6.61632 0.639506 10.3460 0.0000e+00 SES <--> SES

Thanks a lot!!

mspiegel's picture
Offline
Joined: 07/31/2009 - 15:24
You can get the -2 log

You can get the -2 log likelihood like so:

modelOut <- mxRun(model)
summary(modelOut)$Minus2LogLikelihood

In general, you can use names(summary(modelOut)) to see what values are available from the summary.

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
Also, you set can set

Also, you set can set starting values in OpenMx. For RAM models, they default to 1.0.

Other helpful elements are:
fittedModel@output
summary(model)$SaturatedLikelihood

You can get things from a fitted model by name
mxEval(subModel.whatYouWant, model)

You might also get some value from

http://openmx.psyc.virginia.edu/wiki/faq-openmx#Where_are_the_example_scripts_and_demo_data_files

best,
t

carey's picture
Offline
Joined: 10/19/2009 - 15:38
matifou, there is something

matifou,

there is something funky going on, but not with the algorithms that arrive at a solution. it is with (1): summary(MxRunObject) and (2) with standard errors. both are known problems with OpenMx.

to check this out:
(1) add the following to the code for factorModel:
mxPath(from="Powerless67", to="Powerless71", labels="the3", arrows=2),
(2) substitute:
mxData(S.wh2, type="cov", numObs=932))
for
mxData(S.wh2, type="cov", numObs=nrow(S.wh)))

i now get the same parameter estimates from OpenMx and sem and the same chi-square. summary(Run), however, gives an incorrect number of observed statistics. it says 22 when there are really 21. hence the degrees of freedom for the chi square are off as are the all quantities based on the df.

attached are the R code that i used; a text file for the sem commands (something funky about the line endings is causing my R editor to freeze with read.moments and specify.model); and the output that i got.

best,
greg

Matifou's picture
Offline
Joined: 12/08/2009 - 09:29
Okay! So I had

Okay! So I had only:
mxPath(from="Anomia71", to="Anomia67", labels="the3", arrows=2), #

And I needed in facts:

this is added to original script

mxPath(from="Powerless67", to="Powerless71", labels="the3", arrows=2),

I thought that setting arrow=2 would do the trick!

Thanks a lot Carey for taking time to answer and find the right solution! Now I'm a little bit worried that df are not counted as they should... is this known to the developers?

Thanks a lot (and sorry for replying so late to your answers who came so fast!)

Matifou's picture
Offline
Joined: 12/08/2009 - 09:29
Oh and by te way... a second

Oh and by te way... a second question...

I did not understant why you suggest:
(2) substitute:
mxData(S.wh2, type="cov", numObs=932))
for
mxData(S.wh2, type="cov", numObs=nrow(S.wh)))

From the help file I see rather that it hsoul dbe the number of observations from the sample, right?

Thanks!

tbates's picture
Offline
Joined: 07/31/2009 - 14:25
> it should be the number of

> it should be the number of observations from the sample, right?

Exactly so
mxData(S.wh2, type="cov", numObs=932))

is correct, while
mxData(S.wh2, type="cov", numObs=nrow(S.wh)))

would be right for raw data, where rows are people, but for cov data, nrow is just nVar

Matifou's picture
Offline
Joined: 12/08/2009 - 09:29
Ok, thanks for the

Ok, thanks for the information!

About my second concern, that the DF are not accounted as they should, do you confirm it tbates? Any idea, whether this is known issue to the developers and when this will be solved?

Thanks a lot!!