From: Rense Nieuwenhuis <r.nieuwenhuis_at_student.ru.nl>

Date: Mon 28 Aug 2006 - 05:13:26 EST

> endogenous

}

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. Received on Mon Aug 28 05:31:29 2006

Date: Mon 28 Aug 2006 - 05:13:26 EST

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@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. Received on Mon Aug 28 05:31:29 2006

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.1.8, at Mon 28 Aug 2006 - 10:23:41 EST.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*