Dear Jarrett et al., I apologize for missing Ralf Finne's original posting -- I was out of town at the time.
A few comments on this code: (1) You should be able to get standard errors from the GLS solution (though there's a missing factor 0.5 in the definition of objective.2 that should deflate the standard errors by a factor of sqrt(2), and unfortunately this is what one would get by default), but other output from summary.sem() won't be meaningful. (2) The difference between objective.1 and objective.2 in the original code for sem.default() is that the latter uses an analytic gradient. The implementation here of objective.2 for GLS and ULS fits doesn't supply a gradient and consequently it would be clearer just to use objective.2 <- objective.1 for these cases. Even better would be to set analytic.gradient to FALSE in these cases. Of course, one could figure out the gradient for the GLS and ULS estimators, but that's probably not worth the effort. (3) Most of the output from summary(), etc., including standard errors, for the ULS fit won't be sensible. The ULS "fitting function" never made sense to me since it's not invariant with respect to rescaling of the observed variables, but perhaps I'm missing something. (4) I think that it would make most sense to save the fitting function in the object returned by sem(), and then to tailor the output of other functions such as summary.sem() and vcov.sem() to each case. Regards, John > -----Original Message----- > From: r-help-boun...@r-project.org [mailto:r-help-boun...@r-project.org] On > Behalf Of Jarrett Byrnes > Sent: December-15-09 5:07 PM > To: r-help@r-project.org > Subject: Re: [R] Structural Equation Models(SEM) > > Joerg Everman has a great solution to this. He changed the middle of > the sem.mod code to include a variable, fit, and then used the > following approach around where you define the objectives: > > if (fit=="ml") { > objective.1 <- function(par){ > A <- P <- matrix(0, m, m) > val <- ifelse (fixed, ram[,5], par[sel.free]) > A[arrows.1] <- val[one.head] > P[arrows.2t] <- P[arrows.2] <- val[!one.head] > I.Ainv <- solve(diag(m) - A) > C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J) > Cinv <- solve(C) > f <- sum(diag(S %*% Cinv)) + log(det(C)) > attributes(f) <- list(C=C, A=A, P=P) > f > } > objective.2 <- function(par){ > A <- P <- matrix(0, m, m) > val <- ifelse (fixed, ram[,5], par[sel.free]) > A[arrows.1] <- val[one.head] > P[arrows.2t] <- P[arrows.2] <- val[!one.head] > I.Ainv <- solve(diag(m) - A) > C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J) > Cinv <- solve(C) > f <- sum(diag(S %*% Cinv)) + log(det(C)) > grad.P <- correct * t(I.Ainv) %*% t(J) %*% Cinv %*% (C - S) %* > % Cinv %*% J %*% I.Ainv > grad.A <- grad.P %*% P %*% t(I.Ainv) > gradient <- rep(0, t) > gradient[unique.free.1] <- tapply(grad.A[arrows. > 1.free],ram[ram[,1]==1 & ram[,4]!=0, 4], sum) > gradient[unique.free.2] <- tapply(grad.P[arrows. > 2.free],ram[ram[,1]==2 & ram[,4]!=0, 4], sum) > attributes(f) <- list(C=C, A=A, P=P, gradient=gradient) > f > } > } > if (fit=="gls") { > objective.1 <- function(par){ > A <- P <- matrix(0, m, m) > val <- ifelse (fixed, ram[,5], par[sel.free]) > A[arrows.1] <- val[one.head] > P[arrows.2t] <- P[arrows.2] <- val[!one.head] > I.Ainv <- solve(diag(m) - A) > C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J) > Sinv <- solve(S) > f <- 0.5 * sum(diag( ((S - C) %*% Sinv) %*% ((S - C) %*% > Sinv) )) > attributes(f) <- list(C=C, A=A, P=P) > f > } > objective.2 <- function(par){ > A <- P <- matrix(0, m, m) > val <- ifelse (fixed, ram[,5], par[sel.free]) > A[arrows.1] <- val[one.head] > P[arrows.2t] <- P[arrows.2] <- val[!one.head] > I.Ainv <- solve(diag(m) - A) > C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J) > Sinv <- solve(S) > f <- sum(diag( ((S - C) %*% Sinv) %*% ((S - C) %*% Sinv) )) > # grad.P <- correct * t(I.Ainv) %*% t(J) %*% Cinv %*% (C - S) %* > % Cinv %*% J %*% I.Ainv > # grad.A <- grad.P %*% P %*% t(I.Ainv) > # gradient <- rep(0, t) > # gradient[unique.free.1] <- tapply(grad.A[arrows. > 1.free],ram[ram[,1]==1 & ram[,4]!=0, 4], sum) > # gradient[unique.free.2] <- tapply(grad.P[arrows. > 2.free],ram[ram[,1]==2 & ram[,4]!=0, 4], sum) > # attributes(f) <- list(C=C, A=A, P=P, gradient=gradient) > attributes(f) <- list(C=C, A=A, P=P) > f > } > } > if (fit=="uls") { > objective.1 <- function(par){ > A <- P <- matrix(0, m, m) > val <- ifelse (fixed, ram[,5], par[sel.free]) > A[arrows.1] <- val[one.head] > P[arrows.2t] <- P[arrows.2] <- val[!one.head] > I.Ainv <- solve(diag(m) - A) > C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J) > Sinv <- solve(S) > f <- 0.5 * sum(diag( (S - C) %*% (S - C) )) > attributes(f) <- list(C=C, A=A, P=P) > f > } > objective.2 <- function(par){ > A <- P <- matrix(0, m, m) > val <- ifelse (fixed, ram[,5], par[sel.free]) > A[arrows.1] <- val[one.head] > P[arrows.2t] <- P[arrows.2] <- val[!one.head] > I.Ainv <- solve(diag(m) - A) > C <- J %*% I.Ainv %*% P %*% t(I.Ainv) %*% t(J) > Sinv <- solve(S) > f <- 0.5 * sum(diag( (S - C) %*% (S - C) )) > # grad.P <- correct * t(I.Ainv) %*% t(J) %*% Cinv %*% (C - S) %* > % Cinv %*% J %*% I.Ainv > # grad.A <- grad.P %*% P %*% t(I.Ainv) > # gradient <- rep(0, t) > # gradient[unique.free.1] <- tapply(grad.A[arrows. > 1.free],ram[ram[,1]==1 & ram[,4]!=0, 4], sum) > # gradient[unique.free.2] <- tapply(grad.P[arrows. > 2.free],ram[ram[,1]==2 & ram[,4]!=0, 4], sum) > # attributes(f) <- list(C=C, A=A, P=P, gradient=gradient) > attributes(f) <- list(C=C, A=A, P=P) > f > } > } > > > On Dec 2, 2009, at 10:22 AM, Jeremy Miles wrote: > > > In the world of SEM, GLS has pretty much fallen by the wayside - I > > can't recall anything I've seen arguing for it's use in the past 10 > > years, and I also can't recall anyone using it over ML. The > > recommendations for non-normal distributions tend to be robust-ML, or > > robust weighted least squares. These are more computationally > > intensive, and I *think* that John Fox (author of sem) has written > > somewhere that it wouldn't be possible to implement them within R, > > without using a lower level language - or rather that it might be > > possible, but it would be really, really slow. > > > > However, ML and GLS are pretty similar, if you dug around in the > > source code, you could probably make the change (see, > > http://www2.gsu.edu/~mkteer/discrep.html for example, for the > > equations; in fact GLS is somewhat computationally simpler, as you > > don't need to invert the implied covariance matrix at each iteration). > > However, the fact that it's not hard to make the change, and that no > > one has made the change, is another argument that it's not a change > > that needs to be made. > > > > Jeremy > > > > > > > > 2009/12/2 Ralf Finne <ralf.fi...@novia.fi>: > >> Hi R-colleagues. > >> > >> I have been using the sem(sem) function. It uses > >> maximum likelyhood as optimizing. method. > >> According to simulation study in Umeå Sweden > >> (http://www.stat.umu.se/kursweb/vt07/stad04mom3/?download=UlfHolmberg.pdf > >> Sorry it is in swedish, except the abstract) > >> maximum likelihood is OK for large samples and normal distribution > >> the SEM-problem should be optimized by GLS (Generalized Least > >> Squares). > >> > >> > >> So to the question: > >> > >> Is there any R-function that solves SEM with GLS? > >> > >> > >> Ralf Finne > >> Novia University of Applied Science > >> Vasa Finland > >> > >> ______________________________________________ > >> R-help@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. > >> > > > > > > > > -- > > Jeremy Miles > > Psychology Research Methods Wiki: www.researchmethodsinpsychology.com > > > > ______________________________________________ > > R-help@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. > > ______________________________________________ > R-help@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. ______________________________________________ R-help@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.
