R實現多元線性回歸，主要利用的就是lm()函數
熟悉其他統計回歸量的函數，對做回歸分析也是很有幫助的。
- anova(m)： ANOVA表
- coefficients(m)： 模型的系數
- coef(m)： 跟coefficients(m)一樣
- confint(m)： 回歸系數的置信區間
- deviance(m)： 殘差平方和
- effects(m)： 正交效應向量（Vector of orthogonal effects ）
- fitted(m)： 擬合的Y值向量Vector of fitted y values
- residuals(m)： 模型殘差Model residuals
- resid(m)： 跟residuals(m)一樣
- summary(m)：關鍵統計量，例如R2、F統計量和殘差標准差(σ)
- vcov(m)：主參數的協防差矩陣
以下是R做多元線性回歸的幾個基本步驟：
1.讀入數據，R-STUDIO直接有按鈕，否則就
> zsj <- read.csv("D:/Paper/data/zsj.csv")
數據一般從excel的CSV或者txt里讀取，實現整理好以符合R的數據框的結構
ps1：這塊有很多包提供從不同來源讀取數據的方法，筆者還得慢慢學。。
2.畫相關圖選擇回歸方程的形式
> plot(Y~X1);abline(lm(Y~X1))
> plot(Y~X2);abline(lm(Y~X2))
可見X1與Y的關系是明顯的線性的，X2也類似此處省略
3.做回歸,並檢視回歸結果
> lm.test<-lm(Y~X1+X2,data=zsj)
> summary(lm.test)
Call:
lm(formula = Y ~ X1 + X2, data = zsj)
Residuals:
Min 1Q Median 3Q Max
-0.21286 -0.05896 -0.01450 0.05556 0.30795
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0931750 0.0109333 8.522 5.85e-16 ***
X1 0.0109935 0.0003711 29.625 < 2e-16 ***
X2 0.0099941 0.0010459 9.555 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.08109 on 327 degrees of freedom
Multiple R-squared: 0.7953, Adjusted R-squared: 0.7941
F-statistic: 635.3 on 2 and 327 DF, p-value: < 2.2e-16
可見各項顯著性檢驗都是得到通過的
4.用殘差分析剔除異常點
> plot(lm.test,which=1:4)
得到的四個圖依次為：
4.1普通殘差與擬合值的殘差圖
4.2正態QQ的殘差圖（若殘差是來自正態總體分布的樣本，則QQ圖中的點應該在一條直線上）
4.3標准化殘差開方與擬合值的殘差圖（對於近似服從正態分布的標准化殘差，應該有95%的樣本點落在[-2,2]的區間內。這也是判斷異常點的直觀方法）
4.4cook統計量的殘差圖（cook統計量值越大的點越可能是異常值，但具體閥值是多少較難判別）
從圖中可見，54,65,295三個樣本存在異常，需要剔除。
5.檢驗異方差
5.1GQtest，H0（誤差平方與自變量，自變量的平方和其交叉相都不相關），p值很小時拒絕H0，認為上訴公式有相關性，存在異方差
> res.test<-residuals(lm.test)
> library(lmtest)
> gqtest(lm.test)
Goldfeld-Quandt test
data: lm.test
GQ = 0.9353, df1 = 162, df2 = 162, p-value = 0.6647
5.2BPtest,H0(同方差),p值很小時認為存在異方差
> bptest(lm.test)
studentized Breusch-Pagan test
data: lm.test
BP = 3.0757, df = 2, p-value = 0.2148
兩個檢驗都可以看出異方差不存在，不過為了總結所有情況這里還是做了一下修正。。
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則相反。
轉自：http://blog.sina.com.cn/s/blog_6ee39c3901017fpd.html
其他參考資料：
R入門25招 （其中第20~24招）