[R] Total (un)standardized effects in SEM?
John Fox
jfox at mcmaster.ca
Mon Aug 28 01:16:12 CEST 2006
Dear Rense,
I already wrote a function, appended below, to compute effects for SEMs, as
a method for the effects() generic. I've been thinking about adding this to
the sem package, possibly with asymptotic standard errors computed by the
delta method, as suggested by Michael Sobel. If I do so, I can also
incorporate standardized effects.
As I mentioned I have serious doubts about the utility of these ideas, and
worry about encouraging them, but I guess that I could say the same thing
for SEMs in general. My original purpose in writing the sem package was to
have something in R that I use in teaching.
Regards,
John
---------- snip ---------------
effects.sem <- function(object, ...){
A <- object$A
exog <- apply(A, 1, function(x) all(x == 0))
endog <- !exog
B <- A[endog, endog, drop=FALSE]
C <- A[endog, exog, drop=FALSE]
I <- diag(nrow(B))
IBinv <- solve(I - B)
TotEndog <- IBinv - I
TotExog <- IBinv %*% C
result <- list(B=B, C=C, TotEndog=TotEndog, TotExog=TotExog)
class(result) <- "semEffect"
result
}
print.semEffect <- function(x, digits=5, ...){
B <- x$B
C <- x$C
cat("\nDirect Effects of Exogenous on Endogenous Variables:\n")
print(C, digits=digits)
cat("\n\nDirect Effects of Endogenous Variables on Each Other:\n")
print(B, digits=digits)
cat("\n\nIndirect Effects of Exogenous on Endogenous Variables:\n")
print(x$TotExog - C, digits=digits)
cat("\n\nIndirect Effects of Endogenous Variables on Each Other:\n")
print(x$TotEndog - B, digits=digits)
cat("\n\nTotal Effects of Exogenous on Endogenous Variables:\n")
print(x$TotExog, digits=digits)
cat("\n\nTotal Effects of Endogenous Variables on Each Other:\n")
print(x$TotEndog, digits=digits)
invisible(x)
}
--------------------------------
John Fox
Department of Sociology
McMaster University
Hamilton, Ontario
Canada L8S 4M4
905-525-9140x23604
http://socserv.mcmaster.ca/jfox
--------------------------------
> -----Original Message-----
> From: r-help-bounces at stat.math.ethz.ch
> [mailto:r-help-bounces at stat.math.ethz.ch] On Behalf Of Rense
> Nieuwenhuis
> Sent: Sunday, August 27, 2006 2:13 PM
> To: r-help at stat.math.ethz.ch
> Cc: John Fox
> Subject: Re: [R] Total (un)standardized effects in SEM?
>
> 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]]
>
> ______________________________________________
> R-help at stat.math.ethz.ch 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