Implementation of piecewise structural equation modeling (SEM) in R, including estimation of path coefficients and goodness-of-fit statistics.
For more information, see:
Shipley, Bill. "Confirmatory path analysis in a generalized multilevel context." Ecology 90.2 (2009):
363-368.
Shipley, Bill. "The AIC model selection method applied to path analytic models compared using a
d-separation test." Ecology 94.3 (2013): 560-564.
Version: 0.4.4 (2015-03-01)
Author: Jon Lefcheck jslefche@vims.edu
##Examples
###Load package
# library(devtools)
# install_github("jslefche/piecewiseSEM")
library(piecewiseSEM)
###Load data from Shipley 2009
data(shipley2009)
The data is alternately hosted in Ecological Archives E090-028-S1 (DOI: 10.1890/08-1034.1).
###Create model set
The model corresponds to the following hypothesis (Fig. 2, Shipley 2009);
Models are constructed using a mix of the nlme and lmerTest packages, as in the supplements of Shipley 2009.
# Load required libraries
library(lmerTest)
library(nlme)
# Create list of models
shipley2009.modlist = list(
lme(DD~lat, random = ~1|site/tree, na.action = na.omit,
data = shipley2009),
lme(Date~DD, random = ~1|site/tree, na.action = na.omit,
data = shipley2009),
lme(Growth~Date, random = ~1|site/tree, na.action = na.omit,
data = shipley2009),
glmer(Live~Growth+(1|site)+(1|tree),
family=binomial(link = "logit"), data = shipley2009) )
###Run Shipley tests
get.sem.fit returns a list of the following:
(1) the missing paths, whether these paths are conditional on any other variables in the model, and associated p-values;
(2) the Fisher's C statistic and p-value for the model (derived from a Chi-squared distribution);
(3) the AIC, AICc (corrected for small sample size), and associated d.f. for the model.
The argument add.vars allows you to specify a vector of additional variables whose causal independence you also wish to test. This is useful if you are comparing nested models. Default is NULL.
The argument adjust.p allows you to adjust the p-values returned by the function based on the the total degrees of freedom for the model (see supplementary material, Shipley 2013). Default is FALSE (uses the d.f. reported in the summary table).
(See "p-values and all that" for a discussion of p-values from mixed models using the lmer package.)
get.sem.fit(shipley2009.modlist, shipley2009)
The missing paths output differs from Table 2 in Shipley 2009. However, running each d-sep model by hand yields the same answers as this function, leading me to believe that updates to the lme4 and nlme packages are the cause of the discrepancy. Qualitatively, the interpretations are the same.
###Extract path coefficients
Path coefficients can be either unstandardized or standardized (in units of standard deviation of the mean). Default is FALSE. The function returns a data.frame sorted by increasing p-value.
get.sem.coefs(shipley2009.modlist, shipley2009)
###Generate variance-covariance SEM using lavaan
Generate variance-covariance based SEM from the list of linear mixed models. The resulting object can be treated like any other model object constructed using the package lavaan.
lavaan.model = get.lavaan.sem(shipley2009.modlist, shipley2009)
summary(lavaan.model)
The output shows that the variance-covariance SEM is a worse fit, indicating that a hierarchical piecewise approach is justified.
###Plot partial effect between two variables
One might be interested in the partial effects of one variable on another given covariates in the SEM. The function get.partial.resid returns a data.frame of the partial residuals of y ~ x and plots the partial effect.
###Get R2 for individual models
Return R2 and AIC values for component models in the SEM.
get.model.fits(shipley2009.modlist)