This page is under the contruction! More to come!
piecewiseSEM is the right package if you are thinking to fit a nested SEM models. piecewiseSEM can handles other classes of linear models like Generalized Linear Model, Mixed Effects Model, or Generalized Least Squares, which this function make piecewiseSEM differ from the lavaan package.
Here we will be using a dataset from Grace & Keeley (2006) to look get the effects of biotic and abiotic factors on plant species richness in shrublands.
Grace, J. B., & Keeley, J. E. (2006). A structural equation model analysis of postfire plant diversity in California shrublands. Ecological Applications, 16(2), 503–514. https://doi.org/10.1890/1051-0761(2006)016[0503:ASEMAO]2.0.CO;2
#here we will be using two packages piecewiseSEM and nlme
library(piecewiseSEM)##
## This is piecewiseSEM version 2.3.0.
##
##
## Questions or bugs can be addressed to <LefcheckJ@si.edu>.#getting keeley dataset from piecewiseSEM package
data("keeley")
head(keeley)## distance elev abiotic age hetero firesev cover rich
## 1 53.40900 1225 60.67103 40 0.757065 3.50 1.0387974 51
## 2 37.03745 60 40.94291 25 0.491340 4.05 0.4775924 31
## 3 53.69565 200 50.98805 15 0.844485 2.60 0.9489357 71
## 4 53.69565 200 61.15633 15 0.690847 2.90 1.1949002 64
## 5 51.95985 970 46.66807 23 0.545628 4.30 1.2981890 68
## 6 51.95985 970 39.82357 24 0.652895 4.00 1.1734866 34summary(keeley)## distance elev abiotic age
## Min. :37.04 Min. : 60.0 Min. :32.59 Min. : 3.00
## 1st Qu.:39.46 1st Qu.: 202.5 1st Qu.:43.81 1st Qu.:15.00
## Median :51.77 Median : 400.0 Median :48.04 Median :25.00
## Mean :49.23 Mean : 424.7 Mean :49.24 Mean :25.57
## 3rd Qu.:58.40 3rd Qu.: 630.0 3rd Qu.:54.90 3rd Qu.:35.00
## Max. :60.72 Max. :1225.0 Max. :70.46 Max. :60.00
## hetero firesev cover rich
## Min. :0.3842 Min. :1.200 Min. :0.05558 Min. :15.00
## 1st Qu.:0.6246 1st Qu.:3.700 1st Qu.:0.48769 1st Qu.:37.00
## Median :0.6843 Median :4.300 Median :0.63712 Median :50.00
## Mean :0.6833 Mean :4.565 Mean :0.69123 Mean :49.23
## 3rd Qu.:0.7684 3rd Qu.:5.550 3rd Qu.:0.91468 3rd Qu.:62.00
## Max. :0.8779 Max. :9.200 Max. :1.53541 Max. :85.00# fiting structural equation modeling using piecewiseSEM
library(nlme)
model1 = psem (
gls (rich ~ elev + abiotic + age + firesev + cover, data = keeley),
gls (abiotic ~ elev + age + firesev + cover, data = keeley),
gls (cover ~ age + firesev + hetero, data = keeley),
data = keeley
)
summary (model1)##
|
| | 0%
|
|======================= | 33%
|
|=============================================== | 67%
|
|======================================================================| 100%##
## Structural Equation Model of model1
##
## Call:
## rich ~ elev + abiotic + age + firesev + cover
## abiotic ~ elev + age + firesev + cover
## cover ~ age + firesev + hetero
##
## AIC
## 1397.714
##
## ---
## Tests of directed separation:
##
## Independ.Claim Test.Type DF Crit.Value P.Value
## cover ~ elev + ... coef 90 2.8319 0.0058 **
## rich ~ hetero + ... coef 90 4.6231 0.0000 ***
## abiotic ~ hetero + ... coef 90 1.9412 0.0556
##
## --
## Global goodness-of-fit:
##
## Chi-Squared = 26.839 with P-value = 0 and on 3 degrees of freedom
## Fisher's C = 38.479 with P-value = 0 and on 6 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate
## rich elev 0.0071 0.0056 90 1.2790 0.2044 0.1217
## rich abiotic 0.7906 0.1830 90 4.3209 0.0000 0.4019 ***
## rich age -0.1762 0.1213 90 -1.4523 0.1501 -0.1465
## rich firesev -1.2860 0.9461 90 -1.3593 0.1777 -0.1407
## rich cover 6.6921 4.7544 90 1.4075 0.1630 0.1405
## abiotic elev 0.0099 0.0031 90 3.1861 0.0020 0.3341 **
## abiotic age -0.0732 0.0715 90 -1.0240 0.3087 -0.1198
## abiotic firesev -0.6564 0.5563 90 -1.1799 0.2413 -0.1412
## abiotic cover -1.3001 2.8151 90 -0.4618 0.6454 -0.0537
## cover age -0.0053 0.0026 90 -2.0164 0.0469 -0.2102 *
## cover firesev -0.0677 0.0200 90 -3.3925 0.0010 -0.3526 **
## cover hetero -0.5719 0.2570 90 -2.2249 0.0287 -0.2070 *
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method R.squared
## rich none 0.38
## abiotic none 0.15
## cover none 0.26# fiting structural equation modeling using piecewiseSEM
model2 = psem (
gls (rich ~ elev + abiotic + age + firesev + cover, data = keeley),
gls (abiotic ~ elev + age + firesev + cover, data = keeley),
gls (cover ~ age + firesev + elev, data = keeley),
data = keeley
)
summary(model2)##
## Structural Equation Model of model2
##
## Call:
## rich ~ elev + abiotic + age + firesev + cover
## abiotic ~ elev + age + firesev + cover
## cover ~ age + firesev + elev
##
## AIC
## 1413.477
##
## ---
## Tests of directed separation:
##
## No independence claims present. Tests of directed separation not possible.
##
## --
## Global goodness-of-fit:
##
## Chi-Squared = 0 with P-value = 1 and on 0 degrees of freedom
## Fisher's C = NA with P-value = NA and on 0 degrees of freedom
##
## ---
## Coefficients:
##
## Response Predictor Estimate Std.Error DF Crit.Value P.Value Std.Estimate
## rich elev 0.0071 0.0056 90 1.2790 0.2044 0.1217
## rich abiotic 0.7906 0.1830 90 4.3209 0.0000 0.4019 ***
## rich age -0.1762 0.1213 90 -1.4523 0.1501 -0.1465
## rich firesev -1.2860 0.9461 90 -1.3593 0.1777 -0.1407
## rich cover 6.6921 4.7544 90 1.4075 0.1630 0.1405
## abiotic elev 0.0099 0.0031 90 3.1861 0.0020 0.3341 **
## abiotic age -0.0732 0.0715 90 -1.0240 0.3087 -0.1198
## abiotic firesev -0.6564 0.5563 90 -1.1799 0.2413 -0.1412
## abiotic cover -1.3001 2.8151 90 -0.4618 0.6454 -0.0537
## cover age -0.0058 0.0027 90 -2.1674 0.0330 -0.2289 *
## cover firesev -0.0594 0.0203 90 -2.9236 0.0044 -0.3095 **
## cover elev 0.0002 0.0001 90 2.1381 0.0353 0.2026 *
##
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05
##
## ---
## Individual R-squared:
##
## Response method R.squared
## rich none 0.38
## abiotic none 0.15
## cover none 0.26library(lavaan)## This is lavaan 0.6-17
## lavaan is FREE software! Please report any bugs.library(lavaanPlot)
library(piecewiseSEM)
data("keeley")
fit = '
rich ~ elev + abiotic + age + firesev + cover
abiotic ~ elev + age + firesev + cover
cover ~ age + firesev + elev'
lavan1 = sem (fit, data = keeley ) #estimator: ML## Warning in lav_data_full(data = data, group = group, cluster = cluster, :
## lavaan WARNING: some observed variances are (at least) a factor 1000 times
## larger than others; use varTable(fit) to investigatesummary(lavan1, fit.measures = TRUE, standardized = TRUE) # Std.all: standardized regression coefficients## lavaan 0.6.17 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 15
##
## Number of observations 90
##
## Model Test User Model:
##
## Test statistic 0.000
## Degrees of freedom 0
##
## Model Test Baseline Model:
##
## Test statistic 84.779
## Degrees of freedom 12
## P-value 0.000
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000
## Tucker-Lewis Index (TLI) 1.000
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -663.707
## Loglikelihood unrestricted model (H1) -663.707
##
## Akaike (AIC) 1357.414
## Bayesian (BIC) 1394.911
## Sample-size adjusted Bayesian (SABIC) 1347.570
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.000
## P-value H_0: RMSEA <= 0.050 NA
## P-value H_0: RMSEA >= 0.080 NA
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.000
##
## Parameter Estimates:
##
## Standard errors Standard
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## rich ~
## elev 0.007 0.005 1.324 0.186 0.007 0.122
## abiotic 0.791 0.177 4.473 0.000 0.791 0.402
## age -0.176 0.117 -1.503 0.133 -0.176 -0.147
## firesev -1.286 0.914 -1.407 0.159 -1.286 -0.141
## cover 6.692 4.593 1.457 0.145 6.692 0.141
## abiotic ~
## elev 0.010 0.003 3.279 0.001 0.010 0.334
## age -0.073 0.069 -1.054 0.292 -0.073 -0.120
## firesev -0.656 0.541 -1.214 0.225 -0.656 -0.141
## cover -1.300 2.736 -0.475 0.635 -1.300 -0.054
## cover ~
## age -0.006 0.003 -2.217 0.027 -0.006 -0.229
## firesev -0.059 0.020 -2.991 0.003 -0.059 -0.309
## elev 0.000 0.000 2.187 0.029 0.000 0.203
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .rich 139.572 20.806 6.708 0.000 139.572 0.619
## .abiotic 49.636 7.399 6.708 0.000 49.636 0.851
## .cover 0.074 0.011 6.708 0.000 0.074 0.740fitMeasures(lavan1, c("cfi", "rmsea", "srmr", "GFI")) #cfi and GFI > 0.95; rmsea and srmr < 0.08## cfi rmsea srmr gfi
## 1 0 0 1modificationindices(lavan1, minimum.value = 5) #this is to check if your model need to be modifies## [1] lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## <0 rows> (or 0-length row.names)lavaanPlot(model = lavan1, coefs = TRUE, stand = TRUE,
sig = 0.05) #standardized regression paths, showing only paths with p<= .05