目錄

某地區(qū)土壤所含可給態(tài)磷回歸分析

一、研究目的

二、數(shù)據(jù)來源和相關(guān)說明

三、描述性分析

3.1 樣本描述

3.2 數(shù)據(jù)可視化

四、數(shù)據(jù)建模

4.1 回歸模型A

4.2 回歸模型B

4.3 回歸模型B模型診斷

4.4 回歸模型C

五、結(jié)論及建議

5.1 結(jié)論

5.2 建議

六、代碼

---

某地區(qū)土壤所含可給態(tài)磷回歸分析

摘要： 本文建立了多個回歸模型，分析土壤內(nèi)可給態(tài)磷影響因素。經(jīng)過分析得出土壤內(nèi)可給態(tài)林濃度主要與土壤內(nèi)無機(jī)磷濃度、土壤內(nèi)溶于K2CO3溶液但不溶于溴化物水解的有機(jī)磷有關(guān)，且前者對可給態(tài)磷濃度影響較大。逐步回歸模型滿足正態(tài)性與方差齊性，且通過了顯著性檢驗(yàn)。根據(jù)土壤內(nèi)無機(jī)磷濃度增加的0-1變量對土壤內(nèi)可給態(tài)磷濃度不存在顯著性影響?；诖?#xff0c;可以根據(jù)回歸方程模型，對玉米的生產(chǎn)與產(chǎn)量起到指導(dǎo)意義。

一、研究目的

結(jié)合土壤中可給態(tài)磷的濃度，可以預(yù)測玉米體內(nèi)可供態(tài)磷濃度的狀態(tài)，預(yù)測玉米的生長態(tài)勢，指導(dǎo)玉米生產(chǎn)；另一方面，可以求出土壤中無機(jī)磷濃度被期望的范圍，從而對改善玉米的土壤條件，對提高玉米產(chǎn)量起到重要作用。結(jié)合此背景，本文欲分析某地區(qū)土壤所含可給態(tài)磷的影響因素，進(jìn)而對提出合理化建議。

二、數(shù)據(jù)來源和相關(guān)說明

本文依據(jù)表1某地區(qū)土壤所含可給態(tài)磷的情況，建立多元線性回歸方程。數(shù)據(jù)共計(jì)18條，涉及4個變量，分別是：

表 2-1? 變量說明

變量	解釋
X1	土壤內(nèi)所含無機(jī)磷濃度
X2	土壤內(nèi)溶于K2CO3溶液并受溴化物水解的有機(jī)磷
X3 X4	土壤內(nèi)溶于K2CO3溶液但不溶于溴化物水解的有機(jī)磷 以變量X1中位數(shù)為分界點(diǎn)，將X1化為0-1變量
Y	一種在20℃土壤內(nèi)的玉米中的可給態(tài)磷

三、描述性分析

為了獲取對數(shù)據(jù)的直觀了解，本文先對數(shù)據(jù)進(jìn)行了描述性統(tǒng)計(jì)分析。

3.1 樣本描述

表 3-1? 樣本描述

	MEAN	SD	MIN	MED	MAX
X1	11.94	10.15	0.40	10.50	29.90
X2	42.06	13.58	19.00	44.00	65.00
X3	123.00	45.74	37.00	123.50	202.00
Y	81.28	27.00	51.00	77.00	168.00

????? 從表3-1可以得出：土壤內(nèi)所含無機(jī)磷濃度（X1）介于0.40-29.90之間，其平均水平約為11.94（平均值）和10.50（中位數(shù)），變異水平約為10.15（標(biāo)準(zhǔn)差）；土壤內(nèi)溶于K2CO3溶液并受溴化物水解的有機(jī)磷（X2）介于19.00-65.00之間，其平均水平約為42.06（平均值）和44.00（中位數(shù)），變異水平約為13.58（標(biāo)準(zhǔn)差）；土壤內(nèi)溶于K2CO3溶液但不溶于溴化物水解的有機(jī)磷（X3）介于37.00-202.00之間，其平均水平約為123.00（平均值）和123.50（中位數(shù)），變異水平約為45.74（標(biāo)準(zhǔn)差）；一種在20℃土壤內(nèi)的玉米中的可給態(tài)磷（Y）介于51.00-168.00之間，其平均水平約為81.28（平均值）和77.00（中位數(shù)），變異水平約為27.00（標(biāo)準(zhǔn)差）。

3.2 數(shù)據(jù)可視化

????? 為了更直觀分析數(shù)據(jù)分布情況與數(shù)據(jù)間的相關(guān)關(guān)系，本文繪制了各個變量的直方圖與各變量的相關(guān)性熱力圖，結(jié)果如圖3-1和3-2所示。

圖 3-1? 直方圖

????? 從圖3-1可以得出：四個均不服從正態(tài)分布，其中無機(jī)磷（X1）濃度主要集中在0~5；土壤內(nèi)溶于K2CO3溶液并受溴化物水解的有機(jī)磷（X2）濃度主要集中在30-60；土壤內(nèi)溶于K2CO3溶液但不溶于溴化物水解的有機(jī)磷（X3）濃度主要集中在100-175；可給態(tài)磷濃度（Y）主要集中在60-100。

圖 3-2? 相關(guān)性熱力圖

????? 由圖3-2可以得出：變量X1與Y之間存在較強(qiáng)的線性相關(guān)關(guān)系，變量X2、X3與Y之間的線性相關(guān)關(guān)系較弱。

四、數(shù)據(jù)建模

4.1 回歸模型A

????? 為了分析可給態(tài)磷濃度（Y）的影響因素，本文首先建立了全部變量的多元線性回歸方程A，即以Y作為因變量，X1-X3作為自變量，建立線性回歸方程：

Y=β0+β1X1+β2X2+β3X3+?

模型結(jié)果如表4-1所示。

表 4-1? 回歸模型A結(jié)果

	Estimate	Std.Error	t value	Pr(>|t|)
(Intercept)	43.650	18.054	2.418	0.030*
X1	1.785	0.540	3.308	0.052**
X2	-0.083	0.420	-0.198	0.846
X3	0.161	0.112	1.443	0.171
R-squared	0.5493
Adjusted R-squared	0.4527
p-value	0.009227

????? 從表4-1可以得出：在顯著性水平α=0.05 下，三個變量中，只有X1的P值<0.05，即對Y有顯著影響，而變量X2和X3的P值>0.05，即對Y不存在顯著性影響。模型的檢驗(yàn)P值為0.009<0.05，即通過檢驗(yàn)，說明方程具有統(tǒng)計(jì)學(xué)意義。

4.2 回歸模型B

????? 由于模型A中X2和X3對Y不存在顯著性影響，因而本文又建立了逐步回歸模型B，對自變量進(jìn)行選擇，模型結(jié)果如表4-2所示。

表4-2 回歸模型B結(jié)果

	Estimate	Std.Error	t value	Pr(>|t|)
(Intercept)	41.479	13.883	2.988	0.009**
X1	1.737	0.467	3.721	0.002**
X3	0.155	0.104	1.494	0.156
R-squared	0.5481
Adjusted R-squared	0.4878
p-value	0.002589

????? 由表4-2可以得出：經(jīng)過逐步回歸分析，剔除了變量X2，保留了變量X1和X3。模型P值為0.003<0.05，通過檢驗(yàn)?；貧w方程為：

Y=41.479+1.737X1+0.155X3 ??? (1)

即X1每增加1，因變量Y增加1.737；X3每增加1，因變量Y增加0.155，自變量X1對Y影響更大。

4.3 回歸模型B模型診斷

????? 首先對自變量進(jìn)行共線性診斷，結(jié)果如表4-3所示。X1和X3的VIF值均小于4，即變量之間不存在多重共線性。

表 4-3? 共線性檢驗(yàn)

	VIF
X1	1.023639
X3	1.023639

????? 其次，對模型的正態(tài)性、方差齊性與異常值進(jìn)行診斷，結(jié)果如表4-1所示。由QQ圖可以得出，殘差項(xiàng)基本服從正態(tài)分布；由左下角圖形可以得出，方差基本保持水平，即滿足等方差性；由右下角圖形可以得出，大部分點(diǎn)的cook距離均小于0.5，只有第17個點(diǎn)的cook距離>0.5，位于1.0-1.5之間，可以考慮剔除此點(diǎn)，以優(yōu)化模型。

圖 4-1? 回歸模型B診斷

4.4 回歸模型C

????? 考慮到自變量X1對Y影響較大，本文又根據(jù)變量X1增設(shè)變量X4，具體操作是以X1中位數(shù)為分界點(diǎn)將X1數(shù)據(jù)分為兩類，若>X1中位數(shù)，則X4=1，否則X4=0，即變量X4為0-1變量。將變量X4那個模型，建立回歸模型C。模型結(jié)果如表4-4所示。

表4-4? 回歸模型C結(jié)果

	Estimate	Std.Error	t value	Pr(>|t|)
(Intercept)	39.069	17.687	2.209	0.04574*
X1	2.850	0.904	3.153	0.00763**
X2	-0.088	0.405	-0.218	0.83048
X3	0.200	0.111	1.802	0.09475
X4	-25.363	17.609	-1.44	0.17342
R-squared	0.6113
Adjusted R-squared	0.4918
p-value	0.01069

由表4-4可以得出：模型P值為0.01<0.05，即在顯著性水平α=0.05 下，模型C具有統(tǒng)計(jì)學(xué)意義；自變量X1-X4中，只有變量X1通過了檢驗(yàn)，即對因變量Y存在顯著性影響，X2-X3未通過檢驗(yàn)。下面考慮對變量X1-X4進(jìn)行選擇，即建立逐步回歸模型，模型結(jié)果如表4-5所示。

表4-5? 逐步回歸模型結(jié)果

	Estimate	Std.Error	t value	Pr(>|t|)
(Intercept)	41.479	13.883	2.988	0.009**
X1	1.737	0.467	3.721	0.002**
X3	0.155	0.104	1.494	0.156
R-squared	0.5481
Adjusted R-squared	0.4878
p-value	0.002589

??由表4-5可以得出，經(jīng)過逐步回歸分析，剔除了變量X2和X4，保留了X1和X3，模型結(jié)果與回歸模型B結(jié)果相同，模型解讀與診斷與模型B相同。

五、結(jié)論及建議

5.1 結(jié)論

????? 通過建立多個回歸模型，得出如下結(jié)論：土壤內(nèi)可給態(tài)林濃度主要與土壤內(nèi)無機(jī)磷濃度、土壤內(nèi)溶于K2CO3溶液但不溶于溴化物水解的有機(jī)磷有關(guān)，且前者對可給態(tài)磷濃度影響較大。逐步回歸模型滿足正態(tài)性與方差齊性，且通過了顯著性檢驗(yàn)。根據(jù)變量X1增加的0-1變量X4對土壤內(nèi)可給態(tài)磷濃度不存在顯著性影響。

5.2 建議

????? 土壤內(nèi)所含可給態(tài)磷濃度與土壤內(nèi)無機(jī)磷濃度、土壤內(nèi)溶于K2CO3溶液但不溶于溴化物水解的有機(jī)磷濃度存在線性回歸關(guān)系?？梢愿鶕?jù)土壤內(nèi)兩者濃度，預(yù)測土壤內(nèi)可給態(tài)磷濃度，進(jìn)而預(yù)測玉米體內(nèi)可供態(tài)磷濃度的狀態(tài)與生長態(tài)勢，指導(dǎo)玉米生產(chǎn)；另一方面，可以通過改善土壤內(nèi)所含可給態(tài)磷濃度，改善玉米的土壤條件，進(jìn)而提高玉米產(chǎn)量。

六、代碼

a=read.csv("D:/個人成長/學(xué)業(yè)/課程/大三下課程/統(tǒng)計(jì)模型/作業(yè)/第二次作業(yè)/k2co3.csv",header=1)
a[c(1:5),]N = sapply(a,length)
MU = sapply(a,mean)
SD = sapply(a,sd)
MIN = sapply(a,min)
MED = sapply(a,median)
MAX = sapply(a,max)
result = cbind(N,MU,SD,MIN,MED,MAX)
resultpar(mfrow = c(2,2))
hist(a$X1, xlab = "無機(jī)磷", ylab = "頻數(shù)",main=NULL)
hist(a$X2, xlab = "雙溶有機(jī)磷", ylab = "頻數(shù)",main=NULL)
hist(a$X3, xlab = "單溶有機(jī)磷", ylab = "頻數(shù)",main=NULL)
hist(a$Y, xlab = "可給太磷", ylab = "頻數(shù)",main=NULL)
vif(Model.AIC)
library(corrplot)
k=cor(a,use='everything',method='pearson')
par(mfrow=c(1,1))
corrplot(k,addCoef.col = "black")fit = lm(Y~X1+X2+X3, data = a)
summary(fit)Model.AIC=step(fit,trace=F)
summary(Model.AIC)#共線性檢驗(yàn)
library(car)
vif(Model.AIC)#模型診斷
par(mfrow = c(2,2))
plot(Model.AIC,which=1)
plot(Model.AIC,which = 2)
plot(Model.AIC,which = 3)
plot(Model.AIC,which = 4)b=read.csv("D:/個人成長/學(xué)業(yè)/課程/大三下課程/統(tǒng)計(jì)模型/作業(yè)/第二次作業(yè)/k2c032.csv",header=1)fit2 = lm(Y~X1+X2+X3+X4, data = b)
summary(fit2)Model.AIC2=step(fit,trace=F)
summary(Model.AIC2)

---

個人見解，歡迎大家批評指正