Wed, 12/09/2009 - 05:12

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

2

3

4

5 the1 S Anomia67 Anomia67 4.1128138 2.3309929

6 the2 S Powerless67 Powerless67 3.1626654 1.6110255

7

8

9

10

11

12

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!!

Attachment | Size |
---|---|

from sem to openmx.R | 3.29 KB |

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!)

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!

> 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

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!!

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

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_sc...

best,

t

You can get the -2 log likelihood like so:

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