Re: [R] Total (un)standardized effects in SEM?

From: Rense Nieuwenhuis <r.nieuwenhuis_at_student.ru.nl>
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.