[R] Structural Equation Models(SEM)

Jarrett Byrnes byrnes at msi.ucsb.edu
Tue Dec 15 23:07:13 CET 2009


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.Finne at 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 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.
>>
>
>
>
> -- 
> Jeremy Miles
> Psychology Research Methods Wiki: www.researchmethodsinpsychology.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