1、OLS線性回歸的基本原則
最優(yōu)擬合曲線應(yīng)該使各點(diǎn)到直線的距離的平方和（即殘差平方和，簡(jiǎn)稱RSS）最小。
2、OLS線性回歸模型假設(shè)：
- 正態(tài)性：對(duì)于固定的自變量值，因變量值成正態(tài)分布；
- 獨(dú)立性：個(gè)體之間相互獨(dú)立;
- 線性相關(guān)：因變量和自變量之間為線性相關(guān);
- 同方差性：因變量的方差不隨自變量的水平不同而變化，即因變量的方差是不變的。
3、lm()線性回歸函數(shù)
lm(formula, data)
formula中的符號(hào)注釋：
~分割符號(hào)，左邊為因變量，右邊為自變量，例如， z~x+y，表示通過(guò)x和y來(lái)預(yù)測(cè)z
+分割預(yù)測(cè)變量
:表示預(yù)測(cè)變量的交互項(xiàng)，例如，z~x+y+x:y
* 表示所有可能的交互項(xiàng)，例如，z~x*y 展開(kāi)為 z~x+y+x:y
^表示交互項(xiàng)的次數(shù)，例如，z ~ (x+y)^2，展開(kāi)為z~x+y+x:y
.表示包含除因變量之外的所有變量，例如，如果只有三個(gè)變量x，y和z，那么代碼 z~. 展開(kāi)為z~x+y+x:y
-1表示刪除截距項(xiàng)，強(qiáng)制回歸的直線通過(guò)原點(diǎn)
I()從算術(shù)的角度來(lái)解釋括號(hào)中的表達(dá)式，例如，z~y+I(x^2) 表示的擬合公式是 z=a+by+cx2
function 可以在表達(dá)式中應(yīng)用數(shù)學(xué)函數(shù)，例如，log(z) ~x+y
對(duì)于擬合后的模型（lm函數(shù)返回的對(duì)象），可以應(yīng)用下面的函數(shù)，得到模型的更多額外的信息：
summary() 展示擬合模型的詳細(xì)結(jié)果
coefficients()列出模型的參數(shù)（截距項(xiàng)intercept和斜率）
confint()提供模型參數(shù)的置信區(qū)間
residuals() 列出擬合模型的殘差值
fitted() 列出擬合模型的預(yù)測(cè)值
anova() 生成一個(gè)擬合模型的方差分析表
predict() 用擬合模型對(duì)新的數(shù)據(jù)預(yù)測(cè)響應(yīng)變量
4、舉個(gè)例子
4.1 數(shù)據(jù)
library(pacman)
p_load(stargazer)
x <- c(1.5,2.8,4.5,7.5,10.5,13.5,15.1,16.5,19.5,22.5,24.5,26.5)
y <- c(7.0,5.5,4.6,3.6,2.9,2.7,2.5,2.4,2.2,2.1,1.9,1.8)
4.2 模型選擇
4.2.1 直線回歸
lm1 <- lm(y~x)
stargazer(lm1, # 模型對(duì)象、數(shù)據(jù)框、向量、矩陣、列表
type = "text", # 可設(shè)置為text、html、latex
title = "y ~ x", # 表名
style = "default",
summary = F, # 是否做數(shù)據(jù)統(tǒng)計(jì)
align = T, # 是否居中
no.space = T, # 刪除空行
single.row = T,# 使估計(jì)量和置信區(qū)間并排展示
keep.stat = "adj.rsq", # keep.stat="n"只展示樣本量的大小，并移除其他統(tǒng)計(jì)量
font.size = "large",
model.names = T,
header = F # 避免輸出關(guān)于包作者的一些信息
)
## y ~ x
## =======================================
## Dependent variable:
## ---------------------------
## y
## OLS
## ---------------------------------------
## x -0.170*** (0.027)
## Constant 5.603*** (0.435)
## ---------------------------------------
## Adjusted R2 0.776
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.776，不是很理想。
4.2.2 多項(xiàng)式回歸
x1=x
x2=x2
x1 <- x
x2 <- x^2
lm2 <- lm(y~x1+x2)
stargazer(lm2,
type = "text",
title = "y ~ x1 + x2",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## y ~ x1 + x2
## =======================================
## Dependent variable:
## ---------------------------
## y
## OLS
## ---------------------------------------
## x1 -0.466*** (0.057)
## x2 0.011*** (0.002)
## Constant 6.915*** (0.332)
## ---------------------------------------
## Adjusted R2 0.941
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.941，比直線回歸改善很多。
畫(huà)個(gè)擬合圖看看效果：
library(tidyverse)
df1 <- cbind(y=y,x=x) %>% as.data.frame()
df2 <- cbind(y=fitted(lm2),x=x) %>% as.data.frame()
ggplot() +
geom_point(data=df1,aes(x,y)) +
geom_line(data=df2,aes(x,y))
4.2.3 對(duì)數(shù)回歸
x=log(x)
lm3 <- lm(y ~ log(x))
stargazer(lm3,
type = "text",
title = "y ~ log(x)",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## y ~ log(x)
## =======================================
## Dependent variable:
## ---------------------------
## y
## OLS
## ---------------------------------------
## log(x) -1.757*** (0.068)
## Constant 7.364*** (0.169)
## ---------------------------------------
## Adjusted R2 0.984
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.984，比多項(xiàng)式回歸又有了改善。
畫(huà)個(gè)擬合圖看看效果，可以看到擬合效果更好：
df3 <- cbind(y=fitted(lm3),x=x) %>% as.data.frame()
ggplot() +
geom_point(data=df1,aes(x,y)) +
geom_line(data=df3,aes(x,y))
4.2.4 指數(shù)回歸
y=log(y)
lm4 <- lm(log(y) ~ x)
stargazer(lm4,
type = "text",
title = "log(y) ~ x",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## log(y) ~ x
## =======================================
## Dependent variable:
## ---------------------------
## log(y)
## OLS
## ---------------------------------------
## x -0.049*** (0.005)
## Constant 1.760*** (0.075)
## ---------------------------------------
## Adjusted R2 0.907
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.907，比多項(xiàng)式回歸還差。
4.2.5 冪函數(shù)回歸
y=log(y)
x=log(x)
lm5 <- lm(log(y) ~ log(x))
stargazer(lm5,
type = "text",
title = "log(y) ~ log(x)",
style = "default",
summary = F,
align = T,
no.space = T,
single.row = T,
keep.stat = "adj.rsq",
font.size = "large",
model.names = T,
header = F
)
##
## log(y) ~ log(x)
## =======================================
## Dependent variable:
## ---------------------------
## log(y)
## OLS
## ---------------------------------------
## log(x) -0.472*** (0.012)
## Constant 2.191*** (0.030)
## ---------------------------------------
## Adjusted R2 0.993
## =======================================
## Note: *p<0.1; **p<0.05; ***p<0.01
修正后的R2=0.993，比對(duì)數(shù)回歸又有了改善。
畫(huà)個(gè)擬合圖看看效果，可以看到冪函數(shù)法擬合效果最好：
df4 <- cbind(y=exp(fitted(lm5)),x=x) %>% as.data.frame()
ggplot() +
geom_point(data=df1,aes(x,y)) +
geom_line(data=df4,aes(x,y))
5、回歸診斷內(nèi)容
- 樣本是否符合正態(tài)分布假設(shè)？
- 是否存在離群值導(dǎo)致模型產(chǎn)生較大誤差？
- 線性模型是否合理？
- 殘差是否滿足獨(dú)立性、等方差、正態(tài)分布等假設(shè)條件？
- 是否存在多重共線性？
回歸診斷的總結(jié)的函數(shù)為：influence.measures()
5.1 正態(tài)性檢驗(yàn)
- 函數(shù)shapiro.test()
- 直方圖或散點(diǎn)圖目測(cè)檢驗(yàn)
- 殘差計(jì)算函數(shù)residuals()，對(duì)殘差作正態(tài)性檢驗(yàn)（殘差散點(diǎn)圖）
# 殘差正態(tài)性檢驗(yàn)，P值>0.05，W接近1說(shuō)明符合正態(tài)分布
shapiro.test(lm5$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm5$residuals
## W = 0.93033, p-value = 0.3837
# 回歸診斷的標(biāo)準(zhǔn)方法
par(mfrow=c(2,2))
plot(lm5)
正態(tài)性：當(dāng)預(yù)測(cè)變量值固定時(shí)，因變量成正態(tài)分布，則殘差值也應(yīng)該是一個(gè)均值為0的正態(tài)分布。正態(tài)QQ圖（Normal Q-Q）是在正態(tài)分布對(duì)應(yīng)的值下，標(biāo)準(zhǔn)化殘差的概率圖。圖中的數(shù)據(jù)點(diǎn)按對(duì)角直線排列，趨于一條直線，并被對(duì)角直接穿過(guò)，直觀上符合正態(tài)分布。對(duì)于近似服從正態(tài)分布的標(biāo)準(zhǔn)化殘差，應(yīng)該有 95% 的樣本點(diǎn)落在 [-2,2] 區(qū)間內(nèi)；若不是如此，就違反了正態(tài)假設(shè)。
獨(dú)立性：無(wú)法判斷因變量和自變量的值是否相互獨(dú)立，只能從收集的數(shù)據(jù)中來(lái)驗(yàn)證。
線性相關(guān)性：若因變量與自變量線性相關(guān)，那么殘差值與預(yù)測(cè)（擬合）值除了系統(tǒng)誤差之外，就沒(méi)有任何關(guān)聯(lián)。在殘差和擬合圖（Residuals VS Fitted）中，可以看到一條曲線，這暗示著可能需要對(duì)擬合模型加上一個(gè)二次項(xiàng)。（所以我們最后選擇冪函數(shù)回歸）
同方差性：若滿足不變方差假設(shè)，那么在位置尺度圖（Scale-Location）中，水平線周?chē)狞c(diǎn)應(yīng)該隨機(jī)分布。
殘差和杠桿圖（Residuals VS Leverage）提供了特殊的單個(gè)觀測(cè)點(diǎn)（離群點(diǎn)、高杠桿點(diǎn)、強(qiáng)影響點(diǎn)）的信息，本圖中出現(xiàn)了紅色的等高線，說(shuō)明數(shù)據(jù)中存在特別影響回歸結(jié)果的異常點(diǎn)。殘差和杠桿圖的可讀性差，不夠?qū)嵱谩?/strong>
- 離群點(diǎn)：表明擬合回歸模型對(duì)其預(yù)測(cè)效果不佳（產(chǎn)生巨大的殘差）
- 高杠桿點(diǎn)：指自變量因子空間中的離群點(diǎn)（異常值），是由許多異常的自變量值組合起來(lái)的，與因變量沒(méi)有關(guān)系。
- 強(qiáng)影響點(diǎn)：表明它對(duì)模型參數(shù)的估計(jì)產(chǎn)生的影響過(guò)大，非常不成比例。強(qiáng)影響點(diǎn)可以通過(guò)Cook距離（Cook distance）統(tǒng)計(jì)量來(lái)鑒別。看到上面4幅中，每幅圖上都有一些點(diǎn)被特別的標(biāo)記出來(lái)了，這些點(diǎn)是可能存在的異常值點(diǎn)，如果要對(duì)模型進(jìn)行優(yōu)化，我們可以從這些來(lái)入手，從圖中發(fā)現(xiàn)，索引編號(hào)為1、3和12的3個(gè)點(diǎn)在多幅圖中出現(xiàn)。這3個(gè)點(diǎn)可能為異常點(diǎn)，可以從數(shù)據(jù)中去掉這3個(gè)點(diǎn)，再進(jìn)行顯著性檢驗(yàn)和殘差分析。但本次殘差分析的結(jié)果已經(jīng)很好了，所以對(duì)于異常點(diǎn)的優(yōu)化，可能并不能明顯的提升模型的效果。
5.2 殘差正態(tài)性檢驗(yàn)（qqplot）
car包qqplot()函數(shù)提供了精確的正態(tài)假設(shè)檢驗(yàn)方法，置信區(qū)間通過(guò)虛線劃定，當(dāng)絕大多數(shù)點(diǎn)都落在置信區(qū)間時(shí)，說(shuō)明正態(tài)性假設(shè)符合得很好。
5.3 離群值（強(qiáng)影響點(diǎn)）
hv <- hatvalues(lm5);hv
## 1 2 3 4 5 6 7 8 9
## 0.48254231 0.26578892 0.15707660 0.09415762 0.08337304 0.09120380 0.09907002 0.10721053 0.12714202
## 10 11 12
## 0.14898865 0.16407973 0.17936679
# 如何找到最重要或最有影響的觀察結(jié)果？將hat值切分為四分位數(shù)，應(yīng)用95%標(biāo)準(zhǔn)過(guò)濾最異常值，將該過(guò)濾標(biāo)準(zhǔn)應(yīng)用于觀察結(jié)果
hv.vip <- df4$hv[df4$hv > quantile(hv,0.95)];hv.vip
# 具象化
df5 <- cbind(df1,hv=hv)
plot(df5$x,df5$y,col=ifelse(df5$hv > quantile(hv,0.95),'red','blue'))
5.4 誤差的獨(dú)立性
car包提供了durbinWatsonTest()函數(shù)，用于做Durbin-Watson檢驗(yàn)，檢測(cè)誤差的序列相關(guān)性。
durbinWatsonTest(lm5)
lag Autocorrelation D-W Statistic p-value
1 0.1840413 1.184471 0.048
Alternative hypothesis: rho != 0
p值<0.05，不顯著，誤差項(xiàng)之間獨(dú)立，不存在自相關(guān)性。
5.5 線性相關(guān)性
通過(guò)成分殘差圖（component + residual plots）檢查因變量和自變量之間是否呈線性關(guān)系。若圖形存在非線性，則說(shuō)明可能對(duì)預(yù)測(cè)變量的函數(shù)形式建模不夠充分，那么需要添加一些曲線成分，比如多項(xiàng)式，對(duì)一個(gè)或多個(gè)變量進(jìn)行變換（log(x)代替x），或用其他回歸變體形式而不是線性回歸。
crPlots(lm5)
5.6 同方差性
判斷方差是否恒定，nvcTest()函數(shù)生成一個(gè)記分檢驗(yàn)，原假設(shè)為誤方差不變，備擇假設(shè)為誤差方差隨著擬合值水平的變化而變化。若檢驗(yàn)顯著，說(shuō)明存在誤方差不恒定。
spreadLevelPlot()函數(shù)創(chuàng)建了一個(gè)添加了最佳擬合曲線的散點(diǎn)圖，展示了標(biāo)準(zhǔn)化殘差絕對(duì)值與擬合值的關(guān)系。
ncvTest(lm5)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.4044358, Df = 1, p = 0.52481
記分檢驗(yàn)不顯著：p=0.52481，說(shuō)明滿足方差不變假設(shè)，也可以通過(guò)分布水平看到這一點(diǎn)，點(diǎn)在水平的最佳擬合曲線周?chē)孰S機(jī)分布。
spreadLevelPlot(lm5)
5.7 多重共線性（一個(gè)自變量不存在）
在回歸分析的時(shí)候，古典模型中有一個(gè)基本的假定就是自變量之間是不相關(guān)的，但是如果我們?cè)跀M合出來(lái)的回歸模型出現(xiàn)了自變量之間高度相關(guān)的話，可能對(duì)結(jié)果又產(chǎn)生影響，我們稱這個(gè)問(wèn)題為多重共線性。
根據(jù)經(jīng)驗(yàn)表明，多重共線性存在的一個(gè)標(biāo)志就是模型存在較大的標(biāo)準(zhǔn)差和較小的T統(tǒng)計(jì)量，如果一個(gè)模型的可決系數(shù)R2很大，F(xiàn)檢驗(yàn)高度限制，但偏回歸系數(shù)的T檢驗(yàn)幾乎都不顯著，那么模型很可能是存在多重共線性。
- 利用自變量之間的簡(jiǎn)單相關(guān)系數(shù)檢驗(yàn)，一般而言，如果兩個(gè)解釋變量的簡(jiǎn)單相關(guān)系數(shù)較高，則可以認(rèn)為是存在著嚴(yán)重的多重共線性。
- 隨著多重共線性的嚴(yán)重程度增強(qiáng)，方差膨脹因子會(huì)逐漸的變大，在回歸中我們用VIF表示方差膨脹因子，VIF=1/(1-R2)。一般當(dāng)VIF>=10的時(shí)候，我們就可以認(rèn)為存在嚴(yán)重多重共線性。
car包中的vif()函數(shù)可以幫我們算出這個(gè)方差膨脹：
library(car)
vif(lm5)
存在多重共線性的解決方法之一：利用MASS包中的函數(shù)lm.ridge()來(lái)建立嶺回歸模型。
6、回歸診斷（gvlma包）
gvlma()函數(shù)，用于對(duì)線性模型假設(shè)進(jìn)行綜合驗(yàn)證，同時(shí)還能驗(yàn)證偏斜度，峰度和異方差的評(píng)價(jià)。從輸出項(xiàng)中可以看出，假設(shè)檢驗(yàn)的顯著性水平是5%。如果p<0.05，說(shuō)明違反了假設(shè)條件。從Global Stat 的Decision 欄中，可以看到數(shù)據(jù)滿足OLS回歸模型所有統(tǒng)計(jì)假設(shè)的P值（Assumptions acceptable）。
library(gvlma)
gvmodel <- gvlma(lm5)
summary(gvmodel)
##
## Call:
## lm(formula = log(y) ~ log(x))
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.054727 -0.020805 0.004548 0.024617 0.045896
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.19073 0.02951 74.23 4.81e-15 ***
## log(x) -0.47243 0.01184 -39.90 2.34e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0361 on 10 degrees of freedom
## Multiple R-squared: 0.9938, Adjusted R-squared: 0.9931
## F-statistic: 1592 on 1 and 10 DF, p-value: 2.337e-12
##
##
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance = 0.05
##
## Call:
## gvlma(x = lm5)
##
## Value p-value Decision
## Global Stat 8.46263 0.076028 Assumptions acceptable.
## Skewness 0.26282 0.608186 Assumptions acceptable.
## Kurtosis 0.52317 0.469492 Assumptions acceptable.
## Link Function 7.63997 0.005709 Assumptions NOT satisfied!
## Heteroscedasticity 0.03666 0.848150 Assumptions acceptable.
