路徑分析(The path analysis),也是結構方程模型的一種形式,具體信息,請參見見前文:這可能是結構方程模型最詳細的文章了。
直接上例子:
上圖中,作者構建了:不同土壤P庫與P肥料應用之間關係的路徑圖,Errors代表每個變量的方差。圖中的Fertiliser P為磷肥,注意其箭頭指向,它是能直接影響土壤中不同形態P的直接原因。作者根據已有研究和知識體系認為:有機磷(OP)、Fe-Al bound P,Ca-bound P和殘餘P是影響不穩定磷含量的原因(即它們指向Labile P),其餘的箭頭指向也是根據知識體系來做解釋。記住,這個是觀念模型,作者也對這個模型進行了相應的調整。對於建立路徑分析模型的理由,作者是這樣描述的:Path analysis has been shown to be a valuable tool to evaluate the cause-and-effect interrelations between P pools in different soil types and under a variety of management regimes。(參考譯文:路徑分析已被證明是一種有價值的工具,可以評估不同土壤類型和各種管理制度下的P庫之間的因果關係)。
對於結論部分,作者對路徑模型再一次進行解釋:Path analysis uses the data generated by the P fractionation technique in a theoretical model to differentiate between direct and indirect effects from one variable on others.These
interrelations are evaluated as partial correlations between the variables
but, unlike routine multiple regression analysis in which a single dependent variable is considered, path analysis conducts a multivariate multiple regression analysis where several dependent variables are subjected to regression analysis simultaneously on one or more independent variable.(參考譯文:路徑分析使用理論模型中由P區分技術生成的數據來區分一個變量對其他變量的直接和間接影響,這些相互關係被評估為變量之間的偏相關性,但是,與考慮單個因變量的常規多元回歸分析不同,路徑分析進行了多元多元回歸分析,其中多個因變量同時對一個或多個因變量進行回歸分析)。這段話其實很難懂,不懂也不影響,其實最主要的是看作者如何建模和做出解釋的即可。結論部分該文(見參考文獻)對路徑分析的結果並結合實際P形態轉化做了大量闡述,感興趣的可以仔細閱讀該文。其實在我看來,在披露了路徑係數、檢驗p值以及模型擬合度後,結合具體問題進行描述一下即可。方式也和R語言|PLS-PM結構方程模型構建及解釋 的解釋類似,這裡不再重複。
接下來再來看一例:
對於這個路徑分析模型,作者的目的是:挑選的5個變量來預測Langmuir P sorption maxima,即S(max)。單箭頭上的係數P(ij)代表路徑係數(變量與預測變量間),雙向箭頭上的係數r(ij)代表簡單相關係數(變量間)。這裡與前一個例子不同的是:作者假定5個變量間還存在相互作用,互為因果。此外,作者還闡述了變量間相互關係的r值求法,不過這個對於現在的我們來說,借鑑意義不太大了。
本文中涉及的包為plspl包,數據集以vegan內置數據varechem為例。
本文中模型建模篩選步奏如下:首先計算各變量間的相關係數,根據相關性正負來決定箭頭所指方向。
library(vegan)library(PerformanceAnalytics)#加載包,沒有的話需要安裝data("varechem")chart.Correlation(varechem, histogram=TRUE, pch=19)#計算相關係數並出圖如果覺得上圖不好看變量關係,那麼換一個R包出個圖:
install.packages("ggcorrplot")library(ggcorrplot)#加載包# 計算相關性M <- cor(varechem,method = "spearman")
ggcorrplot(M, method = "square", #設置相關性圖展示類型 show.legend = T, #設置是否展示圖例 legend.title = "Corr", #設置圖例的標題 type = "lower", #設置只展示定部panel colors = c("#6D9EC1", "white", "#E46726"), #設置相關性圖的顏色 lab = T, #設置是否顯示顯關係數 hc.order = T)圖如下:
接下來看著相關係數值構建初始模型,取一個方向看,設定絕對值大於0.5的變量進行建模,低於這個值的則捨去。構建的初始模型如下(本文僅做單向箭頭示例,雙向箭頭大家自行嘗試):
注意:上圖只是我人為定義的箭頭指向關係,規定反方向作用無效,並不代表實際情況如此,並沒有實際意義,這裡只是作為例子進行說明而已,且為了省事,比如變量S與K的關係我並沒有進行標註,並且數據是否滿足多元正態性也並沒檢驗,這裡只是為了說明一個理念而已,不必當真。
下面開始在R中建模:
#假定模型建模model <-'Mn ~ Mg+S+P+K+Ca+Humdepth+ Zn Ca ~ P+K+Mg+S+Zn'library(vegan)library(lavaan)mypath <- sem(model = model, data = varechem, se = 'bootstrap', bootstrap = 100)#上面這步如果是筆記本的話會很慢,臺式機無壓力,等一會就好,論電腦好壞的重要性mypath #查看簡單結果summary(mypath, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)結果很多,從上到下,取通常關注的幾個數來做解釋,Estimator為
ML代表估計方法使用最大似然估計法;P-value (Chi-square)的值與0.5作比較,如果低於0.5,則表明模型有待提升;Comparative Fit Index (CFI)為擬合度的指標,值在0-1之間,大於0.9的值是「良好」模型擬合的標準,越接近1越好;Akaike (AIC)值通常越低越好,不過需要在有改進模型的情況下進行比較才有用;Regressions下代表了各項回歸係數與p值。fitmeasures(mypath, c('chisq', 'rmsea', 'cfi', 'aic'))#查看模型擬合度library(semPlot) #沒有這個包則需要安裝
semPaths(mypath, what = 'std', layout = 'tree', residuals = FALSE, edge.label.cex = 1)CFI指數為1,rmsea指數為0,chisq為0.731,說明這個模型達到「非常好」的程度。但需要指出的是,由於變量已經進行挑選,存在人為確定的因素,一些被捨去的因子是否更重要且不能丟失,選取的指標是否具有現實意義等問題,需要數據分析的人去做取捨,因為數學方法是死的,它只會按照按照你給定的步奏進行輸出,是否需要考慮重新建模,增添指標,還是要根據實際研究情況來進行分析取捨。如果想要改進模型,那麼可以使用 modificationindices()命令來計算模型的MI值,排序後進行挑選改進模型,變量的MI值越大,說明越需要將該變量加入模型中,但並非添加得越多越好,完全擬合的模型還不如不要,畢竟還是要回到解釋模型背後隱藏的實際意義,因此還是要結合專業背景知識進行取捨。
mf <- modificationindices(mypath)mf <- mf[order(mf$mi, decreasing = TRUE), ]head(mf)說明:上圖表示了「可以添加Ca和Humdepth和Mn之間的關係」,用以改進模型。大家可以嘗試一下。
#加載包,導入數據,計算相關性,篩選建模變量library(vegan)library(PerformanceAnalytics)#加載包,沒有的話需要安裝data("varechem")chart.Correlation(varechem, histogram=TRUE, pch=19)#計算相關係數並出圖 install.packages("ggcorrplot")library(ggcorrplot)#加載包# 計算相關性M <- cor(varechem,method = "spearman")
ggcorrplot(M, method = "square", #設置相關性圖展示類型 show.legend = T, #設置是否展示圖例 legend.title = "Corr", #設置圖例的標題 type = "lower", #設置只展示定部panel colors = c("#6D9EC1", "white", "#E46726"), #設置相關性圖的顏色 lab = T, #設置是否顯示顯關係數 hc.order = T)#假定模型建模model <-'Mn ~ Mg+S+P+K+Ca+Humdepth+ Zn Ca ~ P+K+Mg+S+Zn'library(vegan)library(lavaan)mypath <- sem(model = model, data = varechem, se = 'bootstrap', bootstrap = 100)#上面這步如果是筆記本的話會很慢,臺式機無壓力,等一會就好,論電腦好壞的重要性mypath #查看簡單結果summary(mypath, standardized = TRUE, fit.measures = TRUE, rsquare = TRUE)#直接提取模型擬合度信息並出圖fitmeasures(mypath, c('chisq', 'rmsea', 'cfi', 'aic'))#查看模型擬合度library(semPlot) #沒有這個包則需要安裝
semPaths(mypath, what = 'std', layout = 'tree', residuals = FALSE, edge.label.cex = 1)mf <- modificationindices(mypath)mf <- mf[order(mf$mi, decreasing = TRUE), ]head(mf)#查看MI值排序參考文章:
1.González Jiménez, J.L, Healy M G , Daly K . Effects of fertiliser on phosphorus pools in soils with contrasting organic matter content: A fractionation and path analysis study[J]. Geoderma, 2018, 338:128-135.
2.Kang J , D Hesterberg, Osmond D L . Soil Organic Matter Effects on Phosphorus Sorption: A Path Analysis[J]. Soil Science Society of America Journal, 2009, 73(2):360-366.
3.Assessing linkage between soil phosphorus forms in contrasting tillage systems by path analysis[J]. Soil & Tillage Research, 2018.
R計算:
均值和方差
R數據分析:
一行代碼做相關分析
單因素方差分析
雙因素方差分析
NMDS分析
方差分解分析
結構方程SEM詳解