Dear Allen,
There's enough wrong here that it's easier for me just to fix it.
First, if you're going to fit to a raw moment matrix, then (unless you want
to do regression through the origin), you'll need to specify an intercept:
------- snip ---------
> df <- read.csv('NDVI_lep_data.csv', header=TRUE)
> sem.mod.raw <- rawMoments(~ NDVI + All.S, data = df)
> sem.mod.raw
Raw Moments
(Intercept) NDVI All.S
(Intercept) 1.000000 0.8841920 15.28161
NDVI 0.884192 0.7820657 13.53573
All.S 15.281609 13.5357259 245.71360
N = 29 >
> sem.mod <- specifyModel()
1: NDVI -> All.S, n2S, NA
2: (Intercept) -> All.S, intercept, NA
3: All.S <-> All.S, s2S, NA
4:
Read 3 records
> sem.mod
Path Parameter
1 NDVI -> All.S n2S
2 (Intercept) -> All.S intercept
3 All.S <-> All.S s2S
> sem1 <- sem(sem.mod, sem.mod.raw, N=29, raw=TRUE,
+ fixed.x=c("(Intercept)",
"NDVI"))> summary(sem1)
Model fit to raw moment matrix.
Model Chisquare = 1.0174e-11 Df = 0 Pr(>Chisq) = NA
AIC = 6
AICc = 0.96
BIC = 10.102
CAIC = 1.0174e-11
Normalized Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
-7.72e-13 -3.99e-13 -3.97e-13 -2.63e-13 0.00e+00 0.00e+00
Parameter Estimates
Estimate Std Error z value Pr(>|z|)
n2S 88.262 35.8674 2.4608 0.01386342 All.S <--- NDVI
intercept -62.759 31.7191 -1.9786 0.04786344 All.S <--- (Intercept)
s2S 10.081 2.6474 3.8079 0.00014016 All.S <--> All.S
Iterations = 0
> summary(lm(All.S ~ NDVI, data=df)) # check
Call:
lm(formula = All.S ~ NDVI, data = df)
Residuals:
Min 1Q Median 3Q Max
-5.474 -2.009 -1.337 3.978 5.298
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -62.76 32.87 -1.909 0.0669 .
NDVI 88.26 37.17 2.374 0.0249 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
'.' 0.1 ' ' 1
Residual standard error: 3.291 on 27 degrees of freedom
Multiple R-squared: 0.1727, Adjusted R-squared: 0.1421
F-statistic: 5.638 on 1 and 27 DF, p-value: 0.02494
----------- snip -----------------
Here are other, more convenient, ways of getting the same thing (with output
not shown):
----------- snip -----------------
summary(sem(sem.mod, data=df, raw=TRUE, formula= ~ NDVI + All.S,
fixed.x=c("(Intercept)", "NDVI")))
sem.mod <- specifyEquations()
All.S = intercept*(Intercept) + n2S*NDVI
summary(sem(sem.mod, data=df, raw=TRUE, formula= ~ NDVI + All.S,
fixed.x=c("(Intercept)", "NDVI")))
----------- snip -----------------
To get standardized coefficients, however, you should fit to a covariance
matrix, rather than a raw moment matrix. As it turns out, this won't work
currently for your model specified as I've done it above, because the model
contains only one explicit parameter. This is unusual, which is probably why
the problem didn't surface previously. I've fixed specifyModel() so that
it
will work in this situation, which then produces the following:
----------- snip -----------------
> sem.mod.noint <- specifyModel()
1: NDVI -> All.S, n2S
2:
Read 1 record
NOTE: adding 1 variances to the model
> sem2 <- sem(sem.mod.noint, data=df, formula= ~ NDVI + All.S,
+ fixed.x="NDVI")> summary(sem2)
Model Chisquare = 0 Df = 0 Pr(>Chisq) = NA
Chisquare (null model) = 5.3098 Df = 1
Goodness-of-fit index = 1
AIC = 4
AICc = 0.46154
BIC = 6.7346
CAIC = 0
Normalized Residuals
Min. 1st Qu. Median Mean 3rd Qu. Max.
0 0 0 0 0 0
Parameter Estimates
Estimate Std Error z value Pr(>|z|)
n2S 88.262 36.5022 2.4180 0.01560679 All.S <--- NDVI
V[All.S] 10.441 2.7905 3.7417 0.00018281 All.S <--> All.S
Iterations = 0
> standardizedCoefficients(sem2)
Std. Estimate
n2S n2S 0.4156191 All.S <--- NDVI
V[All.S] V[All.S] 0.8272607 All.S <--> All.S
1.0000000 NDVI <--> NDVI
> with(df, cor(NDVI, All.S)) # check
[1] 0.4156191
----------- snip -----------------
I'll fix specifyModel() so that it works for a model with one explicit
parameter in the next version of the sem package. In the meantime, you can
duplicate what I've done as follows (output not shown):
----------- snip -----------------
sem.mod.noint.2 <- specifyModel()
NDVI -> All.S, n2S
All.S <-> All.S, err.All.S
sem3 <- sem(sem.mod.noint.2, data=df, formula= ~ NDVI + All.S,
fixed.x="NDVI")
summary(sem3)
----------- snip -----------------
Best,
John
-------------------------------------------------------
John Fox
Senator William McMaster Professor of Social Statistics
Department of Sociology
McMaster Univeristy
Hamilton, Ontario, Canada
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at
r-project.org]
> On Behalf Of Hurlbert, Allen Hartley
> Sent: August-03-12 5:37 PM
> To: r-help at r-project.org
> Subject: [R] SEM standardized path coefficients
>
> Hello,
>
> I have conducted an SEM in which the resultant standardized path
> coefficients are much higher than would be expected from the raw
> correlation matrix. To explore further, I stripped the model down to a
simple> bivariate relationship between two variables (NDVI, and species richness),
> where it's my understanding that the SEM's standardized path
coefficient
> should equal the correlation coefficient.
>
> The datafile can be downloaded at
> http://labs.bio.unc.edu/Hurlbert/pubs/NDVI_lep_data.csv
> The model specification file can be downloaded at
> http://labs.bio.unc.edu/Hurlbert/pubs/SEM_NDVI_AllS_model.txt
> and is simply:
>
> NDVI -> All.S, n2S, NA
> All.S -> All.S, S2S, NA
> NDVI <-> NDVI, n2n, 1
>
> df = read.csv('NDVI_lep_data.csv',header=T)
> cor(df[,c('NDVI','All.S')])
>
> yields:
>
>
> NDVI All.S
>
> NDVI 1.0000000 0.4156191
>
> All.S 0.4156191 1.0000000
>
> But, conducting an SEM yields:
>
> sem.mod = specifyModel('SEM_NDVI_AllS_model.txt')
> sem.mod.cov = rawMoments(~ NDVI + All.S, data = df) sem.mod.cov >
sem.mod.cov[-1,-1] sem.mod.cov
>
>
> NDVI All.S
>
> NDVI 0.7820657 13.53573
>
> All.S 13.5357259 245.71360
>
> sem1 = sem(sem.mod, sem.mod.cov, N=29)
> stdCoef(sem1)
>
>
> n2S n2S 0.97643950 All.S <--- NDVI
>
> n2n n2n 1.00000000 NDVI <--> NDVI
>
> s2S s2S 0.04656591 All.S <--> All.S
>
> I am using version 3.0 of sem on R-14.0.
>
> Thanks for any insight anyone can provide.
>
> Allen Hurlbert, Ph.D.
> Department of Biology
> University of North Carolina
> Chapel Hill, NC 27599-3280
> Tel: 919.843.9930
> Wilson Hall 331
>
>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> 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.