[R] Confirmatory factor analysis problems using sem package (works in Amos)
Solomon Messing
messing at stanford.edu
Tue Jul 28 02:00:47 CEST 2009
Dear John,
Would it possible to use a different optimizer with the sem package?
Perhaps optim(..., method = c("Nelder-Mead", "BFGS", "CG", "L-BFGS-B",
"SANN"),...) for example?
Thank you very much,
-Solomon
> -----Original Message-----
> From: John Fox [mailto:jfox at mcmaster.ca]
> Sent: Friday, May 22, 2009 6:25 AM
> To: 'S. Messing'
> Cc: r-help at r-project.org
> Subject: RE: [R] Confirmatory factor analysis problems using sem package
> (works in Amos)
>
> Dear Solomon,
>
> > -----Original Message-----
> > From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
> On
> > Behalf Of S. Messing
> > Sent: May-22-09 1:27 AM
> > To: r-help at r-project.org
> > Subject: [R] Confirmatory factor analysis problems using sem package
> (works
> > in Amos)
> >
> >
> > Hello all,
> >
> > I'm trying to replicate a confirmatory factor analysis done in Amos.
>
> Occasionally in an ill-conditioned problem, one program will produce a
> solution and another won't. As a general matter, I'd expect Amos to be
> more
> robust than sem() since Amos is written specifically for SEMs, while sem()
> uses nlm(), a general-purpose optimizer.
>
> > The
> > idea is to compare a one-factor and a two-factor model. I get the
> following
> > warning message when I run either model:
> >
> > "Could not compute QR decomposition of Hessian.
> > Optimization probably did not converge."
> >
> > I have no idea what to do here.
>
> A general strategy is to set debug=TRUE in the call to sem() and see what
> happens in the optimization. Then there are several things that you can do
> to try to get the optimization to converge; see ?sem. In this case,
> however,
> I wasn't able to get a solution.
>
> The one-factor model is equivalent to a one-factor exploratory FA, which
> can
> be fit by ML using factanal():
>
> > factanal(factors=1, covmat=correl, n.obs=1100)
>
> Call:
> factanal(factors = 1, covmat = correl, n.obs = 1100)
>
> Uniquenesses:
> pvote jmposaff jmnegaff boposaff bonegaff
> obama.therm mccain.therm oddcand.D evencand.D
> 0.100 0.496 0.497 0.277 0.397
> 0.129 0.312 0.466 0.585
>
> Loadings:
> Factor1
> pvote -0.949
> jmposaff 0.710
> jmnegaff -0.709
> boposaff -0.850
> bonegaff 0.777
> obama.therm -0.934
> mccain.therm 0.830
> oddcand.D 0.731
> evencand.D 0.645
>
> Factor1
> SS loadings 5.744
> Proportion Var 0.638
>
> Test of the hypothesis that 1 factor is sufficient.
> The chi square statistic is 1710.03 on 27 degrees of freedom.
> The p-value is 0
>
> As you can see, the one-factor model fits the data very poorly (as does a
> two-factor EFA); I suspect, but am not sure, that this is the source of
> the
> problem in sem(). I couldn't get a solution from sem() even when I used
> the
> factanal() solution as start values.
>
>
> > I believe posters reported the same
> > problem.
>
> In almost all cases, the models haven't been properly specified, which is
> not the case here. Here, the model just doesn't fit the data.
>
> > It seems that the ability to invert the correlation matrix may
> > have something to do with this error, but solve(correl) yields a
> solution.
>
> No, the input correlation matrix is positive-definite. sem() would have
> complained if it were not:
>
> > eigen(correl, only.values=TRUE)
> $values
> [1] 6.12561630 0.82418329 0.71616585 0.51263750 0.24467315 0.18248909
> 0.17024374
> [8] 0.13905585 0.08493524
>
>
> I'll keep your problem as a test case to see whether I can produce a
> solution, possibly using a different optimizer -- as I mentioned, sem()
> uses
> nlm().
>
> Regards,
> John
>
>
> >
> > Here are my correlation matrix and model specifications:
> >
> > --------------------------- R CODE BEGIN
> > ------------------------------------------------
> >
> > library(sem)
> > correl <- matrix(
> > c(1.0000000,-0.6657822,0.6702089,0.7997673,-0.7225454,0.8992372,
> > -0.8026491,-0.6715168,-0.5781565,-
> > 0.6657822,1.0000000,-0.5107568,
> > -0.5030886,0.6971188,-
> > 0.6306937,0.7200848,0.5121227,0.4496810,
> > 0.6702089,-0.5107568,1.0000000,0.7276558,-
> > 0.3893792,0.6043672,
> > -0.7176532,-0.5247434,-0.4670362,0.7997673,-
> > 0.5030886,0.7276558,
> > 1.0000000,-0.6251056,0.8164190,-0.6728515,-
> > 0.6371453,-0.5531964,
> > -0.7225454,0.6971188,-0.3893792,-
> > 0.6251056,1.0000000,-0.7760765,
> > 0.6175124,0.5567924,0.4914176,0.8992372,-
> > 0.6306937,0.6043672,
> > 0.8164190,-0.7760765,1.0000000,-0.7315507,-
> > 0.6625136,-0.5814590,
> > -0.8026491,0.7200848,-0.7176532,-
> > 0.6728515,0.6175124,-0.7315507,
> >
> 1.0000000,0.5910688,0.5096898,-0.6715168,0.5121227,-
> > 0.5247434,
> > -0.6371453,0.5567924,-
> > 0.6625136,0.5910688,1.0000000,0.8106496,
> > -0.5781565,0.4496810,-0.4670362,-
> > 0.5531964,0.4914176,-0.5814590,
> > 0.5096898,0.8106496,1.0000000),
> ,nrow=9,ncol=9)
> >
> > rownames(correl) = c("pvote", "jmposaff", "jmnegaff",
> > "boposaff","bonegaff",
> > "obama.therm", "mccain.therm",
> > "oddcand.D", "evencand.D" )
> >
> > colnames(correl) = c("pvote", "jmposaff", "jmnegaff",
> > "boposaff","bonegaff",
> > "obama.therm", "mccain.therm",
> > "oddcand.D", "evencand.D" )
> >
> > #One Factor Model:
> > model.all <- specify.model()
> > allmeasures -> pvote, b1, NA
> > allmeasures -> obama.therm, b2, NA
> > allmeasures -> mccain.therm, b3, NA
> > allmeasures -> jmposaff, b4, NA
> > allmeasures -> jmnegaff, b5, NA
> > allmeasures -> boposaff, b6, NA
> > allmeasures -> bonegaff, b7, NA
> > allmeasures -> oddcand.D, b8, NA
> > allmeasures -> evencand.D, b9, NA
> > allmeasures <-> allmeasures, NA,1
> > pvote <-> pvote, v1, NA
> > obama.therm <-> obama.therm, v2, NA
> > mccain.therm <-> mccain.therm, v3, NA
> > jmposaff <-> jmposaff, v4, NA
> > jmnegaff <-> jmnegaff, v5, NA
> > boposaff <-> boposaff, v6, NA
> > bonegaff <-> bonegaff, v7, NA
> > oddcand.D <-> oddcand.D, v8, NA
> > evencand.D <-> evencand.D, v9, NA
> >
> >
> > sem.all <- sem(model.all, correl, 1100)
> >
> > summary(sem.all)
> >
> > #Two Factor Model:
> > model.vi <- specify.model()
> > verbal -> pvote, b1, NA
> > verbal -> obama.therm, b2, NA
> > verbal -> mccain.therm, b3, NA
> > verbal -> jmposaff, b4, NA
> > verbal -> jmnegaff, b5, NA
> > verbal -> boposaff, b6, NA
> > verbal -> bonegaff, b7, NA
> > imp -> oddcand.D, b8, NA
> > imp -> evencand.D, b9, NA
> > imp <-> imp, NA, 1
> > verbal <-> verbal, NA, 1
> > pvote <-> pvote, v1, NA
> > obama.therm <-> obama.therm, v2, NA
> > mccain.therm <-> mccain.therm, v3, NA
> > jmposaff <-> jmposaff, v4, NA
> > jmnegaff <-> jmnegaff, v5, NA
> > boposaff <-> boposaff, v6, NA
> > bonegaff <-> bonegaff, v7, NA
> > oddcand.D <-> oddcand.D, v8, NA
> > evencand.D <-> evencand.D, v9, NA
> > imp <-> verbal, civ, NA
> >
> > sem.vi <- sem(model.vi, correl, 1100)
> > summary(sem.vi)
> >
> > --------------------------- R CODE END
> > ------------------------------------------------
> >
> > Thanks very much.
> >
> > -Solomon
> > --
> > View this message in context: http://www.nabble.com/Confirmatory-factor-
> > analysis-problems-using-sem-package-%28works-in-Amos%29-
> > tp23664618p23664618.html
> > Sent from the R help mailing list archive at Nabble.com.
> >
> > ______________________________________________
> > R-help at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-help
> > PLEASE do read the posting guide
> http://www.R-project.org/posting-guide.html
> > and provide commented, minimal, self-contained, reproducible code.
>
More information about the R-help
mailing list