R語言回歸分析總結

04-28

一元線性回歸

X <- matrix(c(194.5, 20.79, 1.3179, 131.79,194.3, 20.79, 1.3179, 131.79,197.9, 22.40, 1.3502, 135.02,198.4, 22.67, 1.3555, 135.55,199.4, 23.15, 1.3646, 136.46,199.9, 23.35, 1.3683, 136.83,200.9, 23.89, 1.3782, 137.82,201.1, 23.99, 1.3800, 138.00,201.4, 24.02, 1.3806, 138.06,201.3, 24.01, 1.3805, 138.05,203.6, 25.14, 1.4004, 140.04,204.6, 26.57, 1.4244, 142.44,209.5, 28.49, 1.4547, 145.47,208.6, 27.76, 1.4434, 144.34,210.7, 29.04, 1.4630, 146.30,211.9, 29.88, 1.4754, 147.54,212.2, 30.06, 1.4780, 147.80),ncol=4, byrow=T,dimnames = list(1:17, c("F", "h", "log", "log100")))forbes<-as.data.frame(X)plot(forbes$F, forbes$log100)

畫自變數和因變數之間的散點圖，數據探索分析。

回歸方程建模。

lm.sol <- lm(log100 ~ F, data=forbes)summary(lm.sol)

參數T檢驗：回歸方程F檢驗

abline(lm.sol)

殘差分析：

y.res<-residuals(lm.sol);plot(y.res)

text(12,y.res[12], labels=12,adj=1.2)

做簡單處理去掉12號樣本點。

i<-1:17; forbes12<-as.data.frame(X[i!=12, ])lm12<-lm(log100~F, data=forbes12)summary(lm12)

去掉12號樣本回歸方程係數沒有太大的變化，但是係數的標準差和殘差標準差有很大的變化。減少約3倍左右，相關係數 $R^{2}$ 也有提高。

預測：

lm.pred<-predict(lm.sol, new, interval="prediction", level=0.95)

多元線性回歸分析回歸

修正擬合：

new.model <- update(old.model, new.formula)

fm5 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = production)

fm6 <- update(fm5, . ~ . + x6)

smf6 <- update(fm6, sqrt(.) ~ .)

逐步回歸：

函數step(),

step(object, scope, scale = 0,direction = c("both", "backward", "forward"),trace = 1, keep = NULL, steps = 1000, k = 2, ...)

它是以AIC信息統計量為準則，通過選擇最小的AIC信息統計量，來達到刪除或增加變數的目的。

add1()

drop1(lm.step)

（回歸係數檢驗，粗略說就是檢驗某個變數Xi的係數是否為0,參數通過t檢驗）

（回歸方程檢驗，簡單說就是檢驗該組數據是否適用線性方程做回歸,回歸方程通過F檢驗）

回歸診斷

X<-scan()679 292 1012 493 582 1156 997 2189 1097 20781818 1700 747 2030 1643 414 354 1276 745 435540 874 1543 1029 710 1434 837 1748 1381 14281255 1777 370 2316 1130 463 770 724 808 790783 406 1242 658 1746 468 1114 413 1787 35601495 2221 1526Y<-scan()0.79 0.44 0.56 0.79 2.70 3.64 4.73 9.50 5.34 6.855.84 5.21 3.25 4.43 3.16 0.50 0.17 1.88 0.77 1.390.56 1.56 5.28 0.64 4.00 0.31 4.20 4.88 3.48 7.582.63 4.99 0.59 8.19 4.79 0.51 1.74 4.10 3.94 0.963.29 0.44 3.24 2.14 5.71 0.64 1.90 0.51 8.33 14.945.11 3.85 3.93lm.sol<-lm(Y~X); summary(lm.sol)

y.rst<-rstandard(lm.sol); y.fit<-predict(lm.sol)plot(y.rst~y.fit)abline(0.1,0.5);abline(-0.1,-0.5)

殘差的Q-Q圖

plot(model, 2) --model是由lm生成的對象

1.讀入數據（R-STUDIO按鈕選擇數據文件）

2.畫相關圖選擇回歸方程的形式

> plot(Y~X1);abline(lm(Y~X1))

> plot(Y~X2);abline(lm(Y~X2))

3.做回歸,建模型，檢查回歸效果。

參數通過t檢驗（回歸係數檢驗，粗略說就是檢驗某個變數Xi的係數是否為0）

回歸方程通過F檢驗（回歸方程檢驗，簡單說就是檢驗該組數據是否適用線性方程做回歸）

修正擬合：

new.model <- update(old.model, new.formula)

fm5 <- lm(y ~ x1 + x2 + x3 + x4 + x5, data = production)

fm6 <- update(fm5, . ~ . + x6)

smf6 <- update(fm6, sqrt(.) ~ .)

逐步回歸：

函數step(),

step(object, scope, scale = 0,direction = c("both", "backward", "forward"),trace = 1, keep = NULL, steps = 1000, k = 2, ...)

它是以AIC信息統計量為準則，通過選擇最小的AIC信息統計量，來達到刪除或增加變數的目的。

add1()

drop1(lm.step)

4.用殘差分析，剔除異常點

殘差是否滿足：獨立性，等方差性，正態性

回歸模型建立後，計算其殘差，並用shapiro.test()函數做殘差的正態性檢驗

y.res<-residuals(lm.sol)

shapiro.test(y.res)

p-value P值小於0.05,我們有足夠理由拒絕其符合正態分布。大於0.05我們沒辦法拒絕其符合正態分布.

得到的四個圖依次為：

4.1普通殘差與擬合值的殘差圖

4.2正態QQ的殘差圖（若殘差是來自正態總體分布的樣本，則QQ圖中的點應該在一條直線上）

4.3標準化殘差開方與擬合值的殘差圖（對於近似服從正態分布的標準化殘差，應該有95%的樣本點落在[-2,2]的區間內。這也是判斷異常點的直觀方法）

4.4cook統計量的殘差圖（cook統計量值越大的點越可能是異常值，但具體閥值是多少較難判別）

5.檢驗異方差

5.1GQtest，H0（誤差平方與自變數，自變數的平方和其交叉相都不相關），p值很小時拒絕H0，認為上訴公式有相關性，存在異方差

> res.test<-residuals(lm.test)

> library(lmtest)

> gqtest(lm.test)

5.2BPtest,H0(同方差),p值很小時認為存在異方差

> bptest(lm.test)

如果p值很大，兩個檢驗都可以看出異方差不存在，不過為了總結所有情況這裡還是做了一下修正。。

6.修正異方差

修正的方法選擇FGLS即可行廣義最小二乘

6.1修正步驟

6.1.1將y對xi求回歸，算出res--u

6.1.2計算log(u^2)

6.1.3做log(u^2)對xi的輔助回歸 log(u^2),得到擬合函數g=b0+b1x1+..+b2x2

6.1.4計算擬合權數1/h=1/exp(g)，並以此做wls估計

> lm.test2<-lm(log(resid(lm.test)^2)~X1+X2,data=zsj)

> lm.test3<-lm(Y~X1+X2,weights=1/exp(fitted(lm.test2)),data=zsj)

> summary(lm.test3)

7.檢驗多重共線性

7.1計算解釋變數相關稀疏矩陣的條件數k，k<100多重共線性程度很小，100<k<1000較強，>1000嚴重

> XX<-cor(zsj[5:6])

> kappa(XX)

[1] 2.223986

7.2尋找共線性強的解釋變數組合

> eigen(XX)#用於發現共線性強的解釋變數組合#

$values

[1] 1.3129577 0.6870423

$vectors

[,1] [,2]

[1,] 0.7071068 -0.7071068

[2,] 0.7071068 0.7071068

8.修正多重共線性---逐步回歸法

> step(lm.test)

Start: AIC=-1655.03

Y ~ X1 + X2

Df Sum of Sq RSS AIC

<none> 2.1504 -1655.0

- X2 1 0.6005 2.7509 -1575.8

- X1 1 5.7714 7.9218 -1226.7

Call:

lm(formula = Y ~ X1 + X2, data = zsj)

Coefficients:

(Intercept) X1 X2

0.093175 0.010994 0.009994

可見X2，X1都不去掉的時候AIC值最小，模型最佳。

ps2：step中可進行參數設置：direction=c("both","forward","backward")來選擇逐步回歸的方向，默認both，forward時逐漸增加解釋變兩個數，backward則相反。