R語言與簡單的回歸分析

回歸模型是計量里最基礎也最常見的模型之一。究其原因，我想是因為在實際問題中我們并不知道總體分布如何，而且只有一組數據，那么試著對數據作回歸分析將會是一個不錯的選擇。

一、簡單線性回歸

簡單的線性回歸涉及到兩個變量：一個是解釋變量，通常稱為x；另一個是被解釋變量，通常稱為y?；貧w會用常見的最小二乘算法擬合線性模型：

yi = β0 + β1xi +εi

其中β0和β1是回歸系數，εi表示誤差。

在R中，你可以通過函數lm()去計算他。Lm()用法如下：

lm(formula, data, subset, weights, na.action,

method = "qr", model = TRUE, x = FALSE, y = FALSE, qr = TRUE,

singular.ok = TRUE, contrasts = NULL, offset, ...)

參數是formula模型公式，例如y ~ x。公式中波浪號（~）左側的是響應變量，右側是預測變量。函數會估計回歸系數β0和β1，分別以截距（intercept）和x的系數表示。

有三種方式可以實現最小二乘法的簡單線性回歸，假設數據wage1（可以通過names函數查看數據框各項名稱）

（1）lm(wage1$wage ~ wage1$educ + wage1$exper)

（2）lm (wage ~ educ + exper, data= wage1)

（3）attach(wage1)

lm(wage~educ+exper)#不要忘記處理完后用detach（）解出關聯

我們以數據wage1為例，可以看到工資與教育水平的線性關系：

運行下列代碼：

library(foreign)

A<-read.dta("D:/R/data/WAGE1.dta")#導入數據

lm(wage~educ,data=A)

>lm(wage~educ,data=A)

Call:

lm(formula = wage~ educ, data = A)

Coefficients:

(Intercept)         educ 

-0.9049      0.5414

當然得到這些數據是不夠的，我們必須要有足夠的證據去證明我們所做的回歸的合理性。那么如何獲取回歸的信息呢？

嘗試運行以下代碼：

result<-lm(wage~educ,data=A)

summary(result)

我們可以得到以下結果：

Call:

lm(formula = wage~ educ, data = A)

Residuals:

Min     1Q     Median     3Q     Max

-5.3396   -2.1501    -0.9674     1.1921    16.6085

Coefficients:

Estimate   Std.Error    t value  Pr(>|t|)   

(Intercept)   -0.90485   0.68497   -1.321   0.187   

educ         0.54136    0.05325 10.167   <2e-16 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

Residual standarderror: 3.378 on 524 degrees of freedom

MultipleR-squared: 0.1648,     AdjustedR-squared: 0.1632

F-statistic: 103.4on 1 and 524 DF,   p-value: < 2.2e-16

解讀上述結果，我們不難看出，單從判決系數R-squared上看，回歸結果是不理想的，但是，從p值來看，我們還是可以得到回歸系數是很顯著地（注意，這里的P<0.05就可以認為拒絕回歸系數為0，即回歸變量與被解釋變量無關的原擇假設，選擇備擇假設）所以說我們的回歸的效果不好但還是可以接受的。當然，這一點也可以通過做散點圖給我們直觀的印象：

但是影響薪酬的因素不只是education，可能還有其他的，比如工作經驗，工作任期。為了更好地解釋影響薪酬的因素，我們就必須用到多元線性回歸。

如果不需要那么多的信息，比如我們只需要F統計量，那么運行anova（result）即可得到方差分析表，更適合閱讀。

另外，再介紹一下函數predict，用法如下：

predict(object, newdata, se.fit = FALSE, scale = NULL, df = Inf, interval = c("none", "confidence", "prediction"), level = 0.95, type = c("response", "terms"), terms = NULL, na.action = na.pass, pred.var = res.var/weights, weights = 1, ...)說明一下，newdata的數據結構是一個數據框。           這里還值得一提的參數時interval，他有三個選項：none代表不作區間預測，僅給出響應變量的估計；confidence是給出E（Y|X=x）的置信區間；prediction是給出真實Y的置信區間，運行代碼你就會發現兩者的差別。出于穩健性考慮，給出自變量取值時，預測Y值通常會采用prediction參數，但是對于X=x時Y的均值的預測就應該用confidence參數（因為回歸模型是Y=βX+e，所以自變量相同，響應變量也未必一樣）

二、多元線性回歸

還是使用lm函數。在公式的右側指定多個預測變量，用加號（+）連接：

> lm(y ~ u + v+ w)

顯然，多元線性回歸是簡單的線性回歸的擴展?？梢杂卸鄠€預測變量，還是用OLS計算多項式的系數。三變量的回歸等同于這個線性模型：

yi = β0 + β1ui +β2vi + β3wi + εi

在R中，簡單線性回歸和多元線性回歸都是用lm函數。只要在模型公式的右側增加變量即可。輸出中會有擬合的模型的系數：

>result1<-lm(wage~educ+exper+tenure,data=A)

>summary(result1)

Call:

lm(formula = wage~ educ + exper + tenure, data = A)

Residuals:

Min     1Q    Median      3Q    Max

-866.29    -249.23   -51.07   189.62   2190.01

Coefficients:

Estimate    Std.Error    t value   Pr(>|t|)   

(Intercept)    -276.240   106.702   -2.589   0.009778 **

educ          74.415      6.287  11.836   <2e-16 ***

exper         14.892      3.253  4.578   5.33e-06 ***

tenure         8.257      2.498  3.306   0.000983 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

Residual standarderror: 374.3 on 931 degrees of freedom

MultipleR-squared: 0.1459,     AdjustedR-squared: 0.1431

F-statistic:    53 on 3 and 931 DF,   p-value: < 2.2e-16

我們將數據稍作平穩化處理，將wage換成log（wage），再來看看。

>plot(wage~educ,data=A)

>A$logwage<-log(A$wage)

>result1<-lm(logwage~educ+exper+tenure,data=A)

>summary(result1)

Call:

lm(formula =logwage ~ educ + exper + tenure, data = A)

Residuals:

Min      1Q    Median      3Q      Max

-2.05802     -0.29645   -0.03265  0.28788   1.42809

Coefficients:

Estimate   Std. Error   t value   Pr(>|t|)   

(Intercept)   0.284360  0.104190   2.729   0.00656**

educ        0.092029   0.007330 12.555  < 2e-16 ***

exper       0.004121   0.001723  2.391  0.01714 * 

tenure      0.022067  0.003094   7.133 3.29e-12 ***

---

Signif.codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05‘.’ 0.1 ‘ ’ 1

Residual standarderror: 0.4409 on 522 degrees of freedom

MultipleR-squared: 0.316,      AdjustedR-squared: 0.3121

F-statistic: 80.39on 3 and 522 DF,   p-value: < 2.2e-16

看得出，平穩化后的數據線性性是更加好的。

下面我們來提取回歸分析的各項統計數據：

一些統計量和參數都被存儲在lm或者summary中

output <-summary (result1)

SSR<- deviance（result1）＃殘差平方和；（另一種方法：RSquared <- output$r.squared)

LL<-logLik(result1) #對數似然統計量

DegreesOfFreedom<-result1$df #自由度

Yhat<- result1$fitted.values#擬合值向量

Resid<- result1$residuals

s<-output$sigma #誤差標準差的估計值（假設同方差）

CovMatrix <-s^2*output$cov #系數的方差－協方差矩陣（與vcov(result1)同）

ANOVA<-anova(result1)#F統計量

confidentinterval<-confint(result1)#回歸系數的置信區間,level=0.95

effects(result1)#計算正交效應向量（Vector of orthogonal effects ）

predict函數用法與一元完全相同

三、檢查結果

檢查回歸結果是一件復雜而痛苦地事情，需要檢驗的東西也很多，當然有不少事是應該在數據進行回歸分析之前就該處理的，比如檢查復共線性；也有處理中需要考慮的，比如模型的選擇，數據的變換；也有事后需要做的，比如殘差正態性檢驗；還有需要關注與特別處理的數據，比如離群點，杠桿點。

這里我們只提最簡單與最常見的事后處理的基本分析。

通過圖形我們可以以一種十分直觀的辦法檢測我們的擬合效果：

plot(result1)

通過擴展包car中的函數來檢測異常值與主要影響因子。代碼如下：
library(car) 
outlierTest(result1)
influence.measures(result1)

在R中，線性回歸計算變得無比簡單，一個lm函數（或glm函數）基本上就擺平了OLS的一切。但擬合數據還僅僅是萬里長征第一步。最終決定成敗的是擬合的模型是否能真正地派上用場。這樣對結果的檢測與分析就顯得尤為重要。