[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