1. 程式人生 > 其它 >很棒的R語言迴歸模型和方差模型

很棒的R語言迴歸模型和方差模型

對於初學者,利用R語言自帶的資料進行練習是不錯的選擇,下面這些模型便是最好的例項。

1、迴歸模型

迴歸模型利用自帶的faithful資料來示例,faithful是某位地質學家在黃石公園旅遊景點"Old Faithful"間歇泉所記錄的噴發資料。這個資料包括兩組向量,它們分別是泉水的持續時間按(eruptions)(以分鐘計)和噴發間隔時間 (waiting)(以分鐘計)。下面我們來簡單畫張它的關係圖。

> data(faithful)
> attach(faithful)
>  names(faithful)
[1] "eruptions" "waiting"  
> plot(eruptions,waiting,col="blue")

從這張圖裡可以發現,waiting和eruptions之間基本呈現出正相關,即隨著這次噴發持續時間的增長,下一次的噴發就是相距越遠。我們繼續嘗試用eruptions來解釋waiting。lm函式就是用來建立線性迴歸模型,命令如下:

> lm(waiting~eruptions)
Call:
lm(formula = waiting ~ eruptions)
Coefficients:
(Intercept)    eruptions  
      33.47        10.73  

並建立了一個屬於線性迴歸模型的物件,並傳回各個變數係數和其他不同的資料。當然,這個變數方便的話還是應該儲存起來。下面可以用plot函式對這個迴歸模型作診斷檢驗。

> par(mfrow=c(2,2))

> plot(lm(waiting~eruptions),col="blue")

指令par(mfrow=c(2,2))可以將R的輸出視窗設定成為2行2列,下次輸入par(mfrow=c(1,1))即可恢復預設設定。

這四張圖裡面顯示一些比較有用的診斷資訊:殘餘圖、正態分點陣圖、曲氏距離等等。關於曲氏距離,我自己是第一次涉及,wiki一大概代表的是每一點對迴歸線的影響力的大小,數值越大表示影響力越大。

2、多元迴歸模型

R的內建檔案stackloss,記錄了由氧化氨氣而製造硝酸的資料。資料包括4列:Air.Flow(空氣流量)、Water.Temp(水溫)、Acid.Conc.(硝酸濃度)、stack.loss(氨氣損失之百分比)。

> data(stackloss) 
> attach(stackloss) 
The following object is masked _by_ .GlobalEnv:
    stack.loss
The following object is masked from package:datasets:
    stack.loss
> stackloss.lm=lm(stack.loss~Air.Flow+Water.Temp+Acid.Conc.) 
> summary(stackloss.lm) 
Call:
lm(formula = stack.loss ~ Air.Flow + Water.Temp + Acid.Conc.)
Residuals:
    Min      1Q  Median      3Q     Max 
-7.2377 -1.7117 -0.4551  2.3614  5.6978 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -39.9197    11.8960  -3.356  0.00375 ** 
Air.Flow      0.7156     0.1349   5.307  5.8e-05 ***
Water.Temp    1.2953     0.3680   3.520  0.00263 ** 
Acid.Conc.   -0.1521     0.1563  -0.973  0.34405    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.243 on 17 degrees of freedom
Multiple R-squared:  0.9136,    Adjusted R-squared:  0.8983 
F-statistic:  59.9 on 3 and 17 DF,  p-value: 3.016e-09

從以上結果能夠得到這個多元線性迴歸模型為:

stack.loss=−39.9197+0.7156Air.Flow+1.2953Water.Temp−0.1521Acid.Conc.

最後一個p−的值非常小(3.016e-09),是表示並非所有的自變數都沒用,但也不是每一個自變數都有用。其中Acid.Conc.的p−值非常高(0.344),因此Acid.Conc.應該首先被移除。重新輸入新的迴歸模型:

> stackloss.lm=lm(stack.loss~Air.Flow+Water.Temp) 
> summary(stackloss.lm)
Call:
lm(formula = stack.loss ~ Air.Flow + Water.Temp)
Residuals:
    Min      1Q  Median      3Q     Max 
-7.5290 -1.7505  0.1894  2.1156  5.6588 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -50.3588     5.1383  -9.801 1.22e-08 ***
Air.Flow      0.6712     0.1267   5.298 4.90e-05 ***
Water.Temp    1.2954     0.3675   3.525  0.00242 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.239 on 18 degrees of freedom
Multiple R-squared:  0.9088,    Adjusted R-squared:  0.8986 
F-statistic: 89.64 on 2 and 18 DF,  p-value: 4.382e-10
> stackloss.lm=lm(stack.loss~Air.Flow+Water.Temp) 
> summary(stackloss.lm)
Call:
lm(formula = stack.loss ~ Air.Flow + Water.Temp)
Residuals:
    Min      1Q  Median      3Q     Max 
-7.5290 -1.7505  0.1894  2.1156  5.6588 
Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -50.3588     5.1383  -9.801 1.22e-08 ***
Air.Flow      0.6712     0.1267   5.298 4.90e-05 ***
Water.Temp    1.2954     0.3675   3.525  0.00242 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.239 on 18 degrees of freedom
Multiple R-squared:  0.9088,    Adjusted R-squared:  0.8986 
F-statistic: 89.64 on 2 and 18 DF,  p-value: 4.382e-10

我們可以看到新的擬合的多元迴歸模型為:

stack.loss=−50.3588+0.6712Air.Flow+1.2954Water.Temp

結果也比較理想,最後我們還是對迴歸模型作診斷檢驗:

> par(mfrow=c(2,2))

> plot(stackloss.lm,col="blue")

從上面的圖來看,第21點和第1點的曲式距離非常大。這樣的情況下,我們優先移除這兩點。

> stackloss.lm=lm(stack.loss~Air.Flow+Water.Temp+Acid.Conc.,subset=c(-4,-21))

> plot(stackloss.lm,col="blue")

移除了1和21點之後,基本上就沒什麼問題了。

3、方差分析模型

R內建資料裡面PlantGrowth記錄了用不同肥料種植植物的重量。

> data(PlantGrowth) 
> attach(PlantGrowth) 
> names(PlantGrowth) 
[1] "weight" "group" 
> group=as.factor(group)

這組資料中一共有3個組別,控制組和兩種肥料種植組。我們首先要將group轉換成因子。然後我們用盒形圖來表示,並做簡要的方差分析。

> plot(group,weight,main="植物重量",xlab="肥料") 
> summary(aov(weight~group)) 
           Df Sum Sq Mean Sq F value Pr(>F)  
group        2  3.766  1.8832   4.846 0.0159 *
Residuals   27 10.492  0.3886                 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> 

通過方差分析我們發現,由於p-的值非常小(0.01591),所以這三個組別的植物的重量有著比較顯著的差別。最後照例進行診斷檢驗。