piecewiseSEM的內建資料集keeley，記錄了加利福尼亞某些地區遭受火災之後，火災地區的位置以及植物（灌木林）多樣性等資訊。
library(piecewiseSEM)
#資料集，詳情 ?keeley
data(keeley)
head(keeley)
distance，地區到海岸的距離；
elev，海拔高度；
abiotic，非生物條件；
age，災前的林分年齡；
hetero，地區內異質性
firesev，火災嚴重性；
cover，灌木覆蓋度；
rich，灌木物種豐富度。
研究人員希望瞭解火災程度和林分年齡的關係，以及災後灌木物種多樣性的恢復模式：林分年齡是否導致更嚴重的火災程度，火災程度嚴重的地區災後灌木叢恢復能力是否更低。
也就是構建如下模型：
如果使用分段SEM，則需分別構建兩組方程：
firesev ~ age，
cover ~ firesev。
所謂“分段”，即每組關係都是獨立估計的，最後合併以生成有關全域性SEM的推論。
接下來，我們使用piecewiseSEM包中的分段SEM方法建模。
線性迴歸是大多數情形下直接考慮的方法，此時擬合分段SEM過程中即包含了各組獨立的線性迴歸模型。
#分段 SEM，詳情 ?psem
#這裡，對於每個獨立的響應方程，直接使用簡單線性迴歸，確定響應關係
#其它情況，如果已知變數間的某種非線性關係，即可以使用非線性模型
keeley_psem <- psem(
lm(firesev ~ age, data = keeley),
lm(cover ~ firesev, data = keeley),
data = keeley)
summary(keeley_psem)
plot(keeley_psem)
概要中包括了分段SEM的擬合度評估指標，以及表徵變數因果關係的引數估計值及其顯著性等資訊。
從概要中我們可以獲知，兩組獨立的線性迴歸模型顯著，表明更老的灌木林以更嚴重的大火為特徵，並與災後植物覆蓋率的低恢復率和豐富度有關。
路徑圖展示了模型中變數間的因果關係，圖中的數值為標準化後的引數估計值（標準化的迴歸係數）。
那麼另一問題，林分年齡和災後植物恢復率之間是否存在直接的顯著因果關係（二者之間是否也可以得到一個直接指向）呢？也就是如下模型：
另一種方法，透過Akaike資訊準則（AIC）評估。由於AIC通常在比較多種模型時使用，因此我們不妨建立含age和cover關係的新模型，並比較它和原模型的AIC值水平。或者，透過全域性的卡方檢驗比較新舊模型的差異是否顯著。
#將 age 和 cover 的因果關係考慮在內，建立新模型
keeley_psem2 <- psem(
lm(cover ~ firesev + age, data = keeley),
lm(firesev ~ age, data = keeley),
data = keeley
)
#比較新（含 age 和 cover 的關係）舊（不含 age 和 cover 的關係）模型
#全域性模型的卡方檢驗
anova(keeley_psem, keeley_psem2)
儘管AIC、BIC值等顯示新模型優於舊模型，但似乎二者數值差異並不很大，暗示新模型相對於舊模型，整體擬合度沒有較大的提升。並且，新舊模型間全域性的卡方檢驗p值不顯著，表明兩個模型沒有顯著差異。綜合考慮，傾向於選擇簡約性的原模型。
此外，對於主要統計量，儘管它們均在summary()概要中有所展示，但也可使用特定函式單獨檢視。
#引數估計值（迴歸係數）
coefs(keeley_psem, standardize = 'none')
coefs(keeley_psem, standardize = 'scale') #可顯示標準化後的
#定向分離測試
dSep(keeley_psem, .progressBar = FALSE)
#Fisher’s C 統計量，p 值等同於上述定向分離測試
fisherC(keeley_psem)
#AIC 值，用於比較多個模型
AIC(keeley_psem, keeley_psem2)
上述過程展示了在分段SEM中使用最簡單的線性迴歸分步擬合模型。
實際應用中，變數間的關係往往比較複雜，有時可能還需引入較複雜的方法。
最後再順便舉一個小示例，在分段SEM中使用廣義線性模型。
#一個帶廣義線性模型的 SEM 示例
#直接用隨機數展示在分段 SEM 中使用到 glm 時的書寫方式了
set.seed(1)
dat <- data.frame(
x = runif(100),
y1 = runif(100),
y2 = rpois(100, 1),
y3 = runif(100)
)
modelList <- psem(
lm(y1 ~ x, data = dat),
glm(y2 ~ x, family = 'poisson', data = dat),
lm(y3 ~ y1 + y2, data = dat),
data = dat
)
summary(modelList, conserve = TRUE)
https://jslefche.github.io/sem_book/local-estimation.html
