R語言與簡單的回歸分析

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

一、簡單線性回歸

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

yi = β0 + β1xi +εi

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

在R中，你可以通過函數(shù)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, ...)

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

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

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

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

（3）attach(wage1)

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

我們以數(shù)據(jù)wage1為例，可以看到工資與教育水平的線性關(guān)系：

運行下列代碼：

library(foreign)

A<-read.dta("D:/R/data/WAGE1.dta")#導(dǎo)入數(shù)據(jù)

lm(wage~educ,data=A)

>lm(wage~educ,data=A)

Call:

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

Coefficients:

(Intercept)         educ 

-0.9049      0.5414

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

嘗試運行以下代碼：

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

summary(result)

我們可以得到以下結(jié)果：

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

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

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

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

另外，再介紹一下函數(shù)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的數(shù)據(jù)結(jié)構(gòu)是一個數(shù)據(jù)框。           這里還值得一提的參數(shù)時interval，他有三個選項：none代表不作區(qū)間預(yù)測，僅給出響應(yīng)變量的估計；confidence是給出E（Y|X=x）的置信區(qū)間；prediction是給出真實Y的置信區(qū)間，運行代碼你就會發(fā)現(xiàn)兩者的差別。出于穩(wěn)健性考慮，給出自變量取值時，預(yù)測Y值通常會采用prediction參數(shù)，但是對于X=x時Y的均值的預(yù)測就應(yīng)該用confidence參數(shù)（因為回歸模型是Y=βX+e，所以自變量相同，響應(yīng)變量也未必一樣）

二、多元線性回歸

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

> lm(y ~ u + v+ w)

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

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

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

>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

我們將數(shù)據(jù)稍作平穩(wěn)化處理，將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

看得出，平穩(wěn)化后的數(shù)據(jù)線性性是更加好的。

下面我們來提取回歸分析的各項統(tǒng)計數(shù)據(jù)：

一些統(tǒng)計量和參數(shù)都被存儲在lm或者summary中

output <-summary (result1)

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

LL<-logLik(result1) #對數(shù)似然統(tǒng)計量

DegreesOfFreedom<-result1$df #自由度

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

Resid<- result1$residuals

s<-output$sigma #誤差標(biāo)準(zhǔn)差的估計值（假設(shè)同方差）

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

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

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

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

predict函數(shù)用法與一元完全相同

三、檢查結(jié)果

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

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

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

plot(result1)

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

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