Dear John,
thank you very much for your reply. The suggestions you make for
calculating the direct and indirect effects are exactly what I was
looking for. Although I'm very new to SEM and not at all very
experienced in R, I tried to put them together in a function (called
decomp) and expanded it to be able to calculate standardized effects
as well. For that, I made a few changes to your
standardized.coefficients() and added the little function std.matrix
() which converts the output of standardized.coefficients. It's not a
very coherent set of functions, but it works (for now). I have
performed a preliminary test using LISREL. The results where identical.
Again, I would like to thank you for the reply and the work on the
SEM package.
With highest regards,
Rense Nieuwenhuis
(Note: The syntax below consists of three functions. The first one
decomposes SEM-objects. It needs the other two functions to work. I'm
sure this can be programmed a lot cleaner, but for now, it serves my
needs.)
decomp <- function(x, std=FALSE)
{
if(std==FALSE)
{
A <- x$A
# unstandardized structural coefficients
}
if(std==TRUE)
{
A <- std.matrix(std.coef(x)) # standardized
structural coefficients
}
exog <- apply(A, 1, function(x) all(x == 0))
endog <- !exog
B <- A[endog, endog, drop=FALSE] # direct
effects, endogenous ->
endogenous
C <- A[endog, exog, drop=FALSE] # direct
effects, exogenous ->
endogenous
I <- diag(nrow(B))
IBinv <- solve(I - B)
total.endo.endo <- IBinv - I # total
effects, endogenous ->
endogenous
total.exo.endo <- IBinv %*% C # total
effects, exogenous ->
endogenous
ind.endo.endo <- total.endo.endo - B # indirect effects,
endogenous -> endogenous
ind.exo.endo <- total.exo.endo - C # indirect
effects, exogenous -
> endogenous
temp <- list( total.endo.endo = total.endo.endo,
total.exo.endo = total.exo.endo,
ind.endo.endo = ind.endo.endo,
ind.exo.endo = ind.exo.endo)
return (temp)
}
std.matrix <- function(x)
{
i <- length(x$D)
ii <- length(x$A)
zero <- rep(0,i^2)
temp.matrix <- matrix(zero,nrow=i,ncol=i)
colnames(temp.matrix) <- x$D
rownames(temp.matrix) <- x$D
for (t in c(1:ii))
{
temp.matrix[x$A[t],x$B[t]] <- x$C[t,1]
}
return(temp.matrix)
}
std.coef <- function(object, digits=5)
{
old.digits <- options(digits = digits)
on.exit(options(old.digits))
P <- object$P
A <- object$A
t <- object$t
par <- object$coeff
par.posn <- object$par.posn
IAinv <- solve(diag(nrow(A)) - A)
C <- IAinv %*% P %*% t(IAinv)
ram <- object$ram
par.names <- rep(" ", nrow(ram))
for (i in 1:t) {
which.par <- ram[, 4] == i
ram[which.par, 5] <- par[i]
par.names[which.par] <- names(par)[i]
}
one.head <- ram[, 1] == 1
coeff <- ram[one.head, 5]
coeff <- coeff * sqrt(diag(C[ram[one.head, 3], ram[one.head,
3]])/diag(C[ram[one.head, 2], ram[one.head, 2]]))
var.names <- rownames(A)
par.code <- paste(var.names[ram[one.head, 2]], c("<---",
"<-->")[ram[one.head, 1]], var.names[ram[one.head, 3]])
coeff <- data.frame(par.names[one.head], coeff, par.code)
names(coeff) <- c(" ", "Std. Estimate", " ")
## From here on added / changed by me
## print(coeff, rowlab = rep(" ", nrow(coeff)), right = FALSE)
temp <- list( A = var.names[ram[one.head,2]],
B = var.names[ram[one.head,3]],
C = coeff[2],
D = colnames(object$A) )
return(temp)
}
On Aug 25, 2006, at 16:30 , John Fox wrote:
> Dear Rense,
>
> (This question was posted a few days ago when I wasn't reading my
> email.)
>
> So-called effect decompositions are simple functions of the structural
> coefficients of the model, which in a model fit by sem() are
> contained in
> the $A component of the returned object. (See ?sem.) One approach,
> therefore, would be to put the coefficients in the appropriate
> locations of
> the estimated Beta, Gamma, Lamda-x, and Lambda-y matrices of the
> LISREL
> model, and then to compute the "effects" in the usual manner.
>
> It should be possible to do this for the RAM formulation of the
> model as
> well, simply by distinguishing exogenous from endogenous variables.
> Here's
> an illustration using model C in the LISREL 7 Manual, pp. 169-177,
> for the
> Wheaton et al. "stability of alienation" data (a common example--I
> happen to
> have an old LISREL manual handy):
>
>> S.wh <- matrix(c(
> + 11.834, 0, 0, 0, 0, 0,
> + 6.947, 9.364, 0, 0, 0, 0,
> + 6.819, 5.091, 12.532, 0, 0, 0,
> + 4.783, 5.028, 7.495, 9.986, 0, 0,
> + -3.839, -3.889, -3.841, -3.625, 9.610, 0,
> + -2.190, -1.883, -2.175, -1.878, 3.552, 4.503),
> + 6, 6)
>>
>> rownames(S.wh) <- colnames(S.wh) <-
> + c
> ('Anomia67','Powerless67','Anomia71','Powerless71','Education','SEI')
>
>>
>> model.wh <- specify.model()
> 1: Alienation67 -> Anomia67, NA, 1
> 2: Alienation67 -> Powerless67, lam1, NA
> 3: Alienation71 -> Anomia71, NA, 1
> 4: Alienation71 -> Powerless71, lam2, NA
> 5: SES -> Education, NA, 1
> 6: SES -> SEI, lam3, NA
> 7: SES -> Alienation67, gam1, NA
> 8: Alienation67 -> Alienation71, beta, NA
> 9: SES -> Alienation71, gam2, NA
> 10: Anomia67 <-> Anomia67, the1, NA
> 11: Anomia71 <-> Anomia71, the3, NA
> 12: Powerless67 <-> Powerless67, the2, NA
> 13: Powerless71 <-> Powerless71, the4, NA
> 14: Education <-> Education, thd1, NA
> 15: SEI <-> SEI, thd2, NA
> 16: Anomia67 <-> Anomia71, the13, NA
> 17: Alienation67 <-> Alienation67, psi1, NA
> 18: Alienation71 <-> Alienation71, psi2, NA
> 19: SES <-> SES, phi, NA
> 20:
> Read 19 records
>>
>> sem.wh <- sem(model.wh, S.wh, 932)
>>
>> summary(sem.wh)
>
> Model Chisquare = 6.3349 Df = 5 Pr(>Chisq) = 0.27498
> Chisquare (null model) = 17973 Df = 15
> Goodness-of-fit index = 0.99773
> Adjusted goodness-of-fit index = 0.99046
> RMSEA index = 0.016934 90 % CI: (NA, 0.05092)
> Bentler-Bonnett NFI = 0.99965
> Tucker-Lewis NNFI = 0.99978
> Bentler CFI = 0.99993
> BIC = -27.852
>
> Normalized Residuals
> Min. 1st Qu. Median Mean 3rd Qu. Max.
> -9.57e-01 -1.34e-01 -4.24e-02 -9.17e-02 6.43e-05 5.47e-01
>
> Parameter Estimates
> Estimate Std Error z value Pr(>|z|)
> lam1 1.02656 0.053424 19.2152 0.0000e+00 Powerless67 <---
> Alienation67
> lam2 0.97089 0.049608 19.5712 0.0000e+00 Powerless71 <---
> Alienation71
> lam3 0.51632 0.042247 12.2214 0.0000e+00 SEI <--- SES
> gam1 -0.54981 0.054298 -10.1258 0.0000e+00 Alienation67 <--- SES
> beta 0.61732 0.049486 12.4746 0.0000e+00 Alienation71 <---
> Alienation67
> gam2 -0.21151 0.049862 -4.2419 2.2164e-05 Alienation71 <--- SES
> the1 5.06546 0.373464 13.5635 0.0000e+00 Anomia67 <--> Anomia67
> the3 4.81176 0.397345 12.1098 0.0000e+00 Anomia71 <--> Anomia71
> the2 2.21438 0.319740 6.9256 4.3423e-12 Powerless67 <-->
> Powerless67
> the4 2.68322 0.331274 8.0997 4.4409e-16 Powerless71 <-->
> Powerless71
> thd1 2.73051 0.517737 5.2739 1.3353e-07 Education <--> Education
> thd2 2.66905 0.182260 14.6442 0.0000e+00 SEI <--> SEI
> the13 1.88739 0.241627 7.8112 5.7732e-15 Anomia71 <--> Anomia67
> psi1 4.70477 0.427511 11.0050 0.0000e+00 Alienation67 <-->
> Alienation67
> psi2 3.86642 0.343971 11.2406 0.0000e+00 Alienation71 <-->
> Alienation71
> phi 6.87948 0.659208 10.4360 0.0000e+00 SES <--> SES
>
> Iterations = 58
>>
>> A <- sem.wh$A # structural coefficients
>> exog <- apply(A, 1, function(x) all(x == 0))
>> endog <- !exog
>
>> (B <- A[endog, endog, drop=FALSE]) # direct effects, endogenous ->
> endogenous
> Anomia67 Powerless67 Anomia71 Powerless71 Education SEI
> Anomia67 0 0 0 0 0 0
> Powerless67 0 0 0 0 0 0
> Anomia71 0 0 0 0 0 0
> Powerless71 0 0 0 0 0 0
> Education 0 0 0 0 0 0
> SEI 0 0 0 0 0 0
> Alienation67 0 0 0 0 0 0
> Alienation71 0 0 0 0 0 0
> Alienation67 Alienation71
> Anomia67 1.0000000 0.000000
> Powerless67 1.0265597 0.000000
> Anomia71 0.0000000 1.000000
> Powerless71 0.0000000 0.970892
> Education 0.0000000 0.000000
> SEI 0.0000000 0.000000
> Alienation67 0.0000000 0.000000
> Alienation71 0.6173153 0.000000
>
>> (C <- A[endog, exog, drop=FALSE]) # direct effects, exogenous ->
> endogenous
> SES
> Anomia67 0.0000000
> Powerless67 0.0000000
> Anomia71 0.0000000
> Powerless71 0.0000000
> Education 1.0000000
> SEI 0.5163168
> Alienation67 -0.5498096
> Alienation71 -0.2115088
>
>> I <- diag(nrow(B))
>> IBinv <- solve(I - B)
>> (Ty <- IBinv - I) # total effects, endogenous -> endogenous
> Anomia67 Powerless67 Anomia71 Powerless71 Education SEI
> Anomia67 0 0 0 0 0 0
> Powerless67 0 0 0 0 0 0
> Anomia71 0 0 0 0 0 0
> Powerless71 0 0 0 0 0 0
> Education 0 0 0 0 0 0
> SEI 0 0 0 0 0 0
> Alienation67 0 0 0 0 0 0
> Alienation71 0 0 0 0 0 0
> Alienation67 Alienation71
> Anomia67 1.0000000 0.000000
> Powerless67 1.0265597 0.000000
> Anomia71 0.6173153 1.000000
> Powerless71 0.5993465 0.970892
> Education 0.0000000 0.000000
> SEI 0.0000000 0.000000
> Alienation67 0.0000000 0.000000
> Alienation71 0.6173153 0.000000
>
>> (Tx <- IBinv %*% C) # total effects, exogenous -> endogenous
> SES
> Anomia67 -0.5498096
> Powerless67 -0.5644124
> Anomia71 -0.5509147
> Powerless71 -0.5348786
> Education 1.0000000
> SEI 0.5163168
> Alienation67 -0.5498096
> Alienation71 -0.5509147
>
>> Ty - B # indirect effects, endogenous -> endogenous
> Anomia67 Powerless67 Anomia71 Powerless71 Education SEI
> Anomia67 0 0 0 0 0 0
> Powerless67 0 0 0 0 0 0
> Anomia71 0 0 0 0 0 0
> Powerless71 0 0 0 0 0 0
> Education 0 0 0 0 0 0
> SEI 0 0 0 0 0 0
> Alienation67 0 0 0 0 0 0
> Alienation71 0 0 0 0 0 0
> Alienation67 Alienation71
> Anomia67 0.0000000 0
> Powerless67 0.0000000 0
> Anomia71 0.6173153 0
> Powerless71 0.5993465 0
> Education 0.0000000 0
> SEI 0.0000000 0
> Alienation67 0.0000000 0
> Alienation71 0.0000000 0
>
>> Tx - C # indirect effects, exogenous -> endogenous
> SES
> Anomia67 -0.5498096
> Powerless67 -0.5644124
> Anomia71 -0.5509147
> Powerless71 -0.5348786
> Education 0.0000000
> SEI 0.0000000
> Alienation67 0.0000000
> Alienation71 -0.3394059
>
> These results agree with those in the LISREL manual (and for
> another example
> there as well), but I haven't checked the method carefully.
>
> It would, of course, be simple to encapsulate the steps above in a
> function,
> but here's a caveat: The idea of indirect and total effects makes
> sense to
> me for a recursive model, and for the exogenous variables in a
> nonrecursive
> model, where they are the reduced-form coefficients (supposing, of
> course,
> that the model makes sense in the first place, which is often
> problematic),
> but not for the endogenous variables in a nonrecursive model. That
> is why I
> haven't put such a function in the sem package; perhaps I should
> reconsider.
>
> Having said that, I'm ashamed to add that I believe that I was the
> person
> who suggested the definition of total and indirect effects
> currently used
> for these models.
>
> Finally, you can get standardized effects similarly by using
> standardized
> structural coefficients. In the sem package, these are computed and
> printed
> by standardized.eoefficients(). This function doesn't return the
> standardized A matrix in a usable form, but could be made to do so.
>
> Regards,
> John
>
> --------------------------------
> John Fox
> Department of Sociology
> McMaster University
> Hamilton, Ontario
> Canada L8S 4M4
> 905-525-9140x23604
> http://socserv.mcmaster.ca/jfox
> --------------------------------
>
>
[[alternative HTML version deleted]]
______________________________________________
[email protected] 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.
