回歸模型是計(jì)量里最基礎(chǔ)也最常見的模型之一。究其原因，我想是因?yàn)樵趯?shí)際問題中我們并不知道總體分布如何，而且只有一組數(shù)據(jù),，那么試著對(duì)數(shù)據(jù)作回歸分析將會(huì)是一個(gè)不錯(cuò)的選擇。

一、簡單線性回歸

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

yi = β0 + β1xi +εi

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

在R中,，你可以通過函數(shù)lm()去計(jì)算他。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,。公式中波浪號(hào)（~）左側(cè)的是響應(yīng)變量，右側(cè)是預(yù)測變量,。函數(shù)會(huì)估計(jì)回歸系數(shù)β0和β1,，分別以截距（intercept）和x的系數(shù)表示。

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

（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)系：

運(yù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ù)去證明我們所做的回歸的合理性,。那么如何獲取回歸的信息呢,？

          嘗試運(yùn)行以下代碼：

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)然,，這一點(diǎn)也可以通過做散點(diǎn)圖給我們直觀的印象：

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

             如果不需要那么多的信息,，比如我們只需要F統(tǒng)計(jì)量，那么運(yùn)行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)是一個(gè)數(shù)據(jù)框,。

二,、多元線性回歸

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

> lm(y ~ u + v+ w)

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

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

              在R中,，簡單線性回歸和多元線性回歸都是用lm函數(shù)。只要在模型公式的右側(cè)增加變量即可,。輸出中會(huì)有擬合的模型的系數(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ù)線性性是更加好的。

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

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

output <-summary (result1)

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

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

DegreesOfFreedom<-result1$df #自由度

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

Resid<- result1$residuals

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

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

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

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

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

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

三,、檢查結(jié)果

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

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

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

plot(result1)

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

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