1当犯、OLS線性回歸的基本原則
最優(yōu)擬合曲線應(yīng)該使各點到直線的距離的平方和（即殘差平方和米绕，簡稱RSS）最小瑟捣。
2馋艺、OLS線性回歸模型假設(shè)：
- 正態(tài)性：對于固定的自變量值，因變量值成正態(tài)分布迈套；
- 獨立性：個體之間相互獨立;
- 線性相關(guān)：因變量和自變量之間為線性相關(guān);
- 同方差性：因變量的方差不隨自變量的水平不同而變化捐祠，即因變量的方差是不變的。
3桑李、lm()線性回歸函數(shù)
lm(formula, data)
formula中的符號注釋：
~分割符號踱蛀，左邊為因變量，右邊為自變量芙扎，例如星岗， z~x+y，表示通過x和y來預(yù)測z
+分割預(yù)測變量
:表示預(yù)測變量的交互項戒洼，例如俏橘，z~x+y+x:y
* 表示所有可能的交互項，例如圈浇，z~x*y 展開為 z~x+y+x:y
^表示交互項的次數(shù)寥掐，例如，z ~ (x+y)^2磷蜀，展開為z~x+y+x:y
.表示包含除因變量之外的所有變量召耘，例如，如果只有三個變量x褐隆，y和z污它，那么代碼 z~. 展開為z~x+y+x:y
-1表示刪除截距項，強制回歸的直線通過原點
I()從算術(shù)的角度來解釋括號中的表達式庶弃，例如衫贬，z~y+I(x^2) 表示的擬合公式是 z=a+by+cx2
function 可以在表達式中應(yīng)用數(shù)學(xué)函數(shù)，例如歇攻，log(z) ~x+y
對于擬合后的模型（lm函數(shù)返回的對象）固惯，可以應(yīng)用下面的函數(shù)，得到模型的更多額外的信息：
summary() 展示擬合模型的詳細結(jié)果
coefficients()列出模型的參數(shù)（截距項intercept和斜率）
confint()提供模型參數(shù)的置信區(qū)間
residuals() 列出擬合模型的殘差值
fitted() 列出擬合模型的預(yù)測值
anova() 生成一個擬合模型的方差分析表
predict() 用擬合模型對新的數(shù)據(jù)預(yù)測響應(yīng)變量
4缴守、舉個例子
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, # 模型對象葬毫、數(shù)據(jù)框、向量屡穗、矩陣贴捡、列表
type = "text", # 可設(shè)置為text、html村砂、latex
title = "y ~ x", # 表名
style = "default",
summary = F, # 是否做數(shù)據(jù)統(tǒng)計
align = T, # 是否居中
no.space = T, # 刪除空行
single.row = T,# 使估計量和置信區(qū)間并排展示
keep.stat = "adj.rsq", # keep.stat="n"只展示樣本量的大小烂斋，并移除其他統(tǒng)計量
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 多項式回歸
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源祈，比直線回歸改善很多。
畫個擬合圖看看效果：
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 對數(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色迂，比多項式回歸又有了改善香缺。
畫個擬合圖看看效果，可以看到擬合效果更好：
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歇僧，比多項式回歸還差图张。
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，比對數(shù)回歸又有了改善诈悍。
畫個擬合圖看看效果祸轮，可以看到冪函數(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)生較大誤差适袜？
- 線性模型是否合理？
- 殘差是否滿足獨立性舷夺、等方差苦酱、正態(tài)分布等假設(shè)條件？
- 是否存在多重共線性给猾？
回歸診斷的總結(jié)的函數(shù)為：influence.measures()
5.1 正態(tài)性檢驗
- 函數(shù)shapiro.test()
- 直方圖或散點圖目測檢驗
- 殘差計算函數(shù)residuals()疫萤，對殘差作正態(tài)性檢驗（殘差散點圖）
# 殘差正態(tài)性檢驗，P值>0.05敢伸，W接近1說明符合正態(tài)分布
shapiro.test(lm5$residuals)
##
## Shapiro-Wilk normality test
##
## data: lm5$residuals
## W = 0.93033, p-value = 0.3837
# 回歸診斷的標(biāo)準方法
par(mfrow=c(2,2))
plot(lm5)
正態(tài)性：當(dāng)預(yù)測變量值固定時扯饶，因變量成正態(tài)分布，則殘差值也應(yīng)該是一個均值為0的正態(tài)分布池颈。正態(tài)QQ圖（Normal Q-Q）是在正態(tài)分布對應(yīng)的值下尾序，標(biāo)準化殘差的概率圖。圖中的數(shù)據(jù)點按對角直線排列饶辙，趨于一條直線蹲诀，并被對角直接穿過，直觀上符合正態(tài)分布弃揽。對于近似服從正態(tài)分布的標(biāo)準化殘差脯爪，應(yīng)該有 95% 的樣本點落在 [-2,2] 區(qū)間內(nèi)；若不是如此矿微，就違反了正態(tài)假設(shè)痕慢。
獨立性：無法判斷因變量和自變量的值是否相互獨立，只能從收集的數(shù)據(jù)中來驗證涌矢。
線性相關(guān)性：若因變量與自變量線性相關(guān)掖举，那么殘差值與預(yù)測（擬合）值除了系統(tǒng)誤差之外，就沒有任何關(guān)聯(lián)娜庇。在殘差和擬合圖（Residuals VS Fitted）中塔次，可以看到一條曲線方篮，這暗示著可能需要對擬合模型加上一個二次項。（所以我們最后選擇冪函數(shù)回歸）
同方差性：若滿足不變方差假設(shè)励负，那么在位置尺度圖（Scale-Location）中藕溅，水平線周圍的點應(yīng)該隨機分布。
殘差和杠桿圖（Residuals VS Leverage）提供了特殊的單個觀測點（離群點继榆、高杠桿點巾表、強影響點）的信息，本圖中出現(xiàn)了紅色的等高線略吨，說明數(shù)據(jù)中存在特別影響回歸結(jié)果的異常點集币。殘差和杠桿圖的可讀性差，不夠?qū)嵱谩?/strong>
- 離群點：表明擬合回歸模型對其預(yù)測效果不佳（產(chǎn)生巨大的殘差）
- 高杠桿點：指自變量因子空間中的離群點（異常值）翠忠，是由許多異常的自變量值組合起來的鞠苟，與因變量沒有關(guān)系。
- 強影響點：表明它對模型參數(shù)的估計產(chǎn)生的影響過大负间，非常不成比例偶妖。強影響點可以通過Cook距離（Cook distance）統(tǒng)計量來鑒別≌＃看到上面4幅中趾访，每幅圖上都有一些點被特別的標(biāo)記出來了，這些點是可能存在的異常值點董虱，如果要對模型進行優(yōu)化扼鞋，我們可以從這些來入手，從圖中發(fā)現(xiàn)愤诱，索引編號為1云头、3和12的3個點在多幅圖中出現(xiàn)。這3個點可能為異常點淫半，可以從數(shù)據(jù)中去掉這3個點溃槐，再進行顯著性檢驗和殘差分析。但本次殘差分析的結(jié)果已經(jīng)很好了科吭，所以對于異常點的優(yōu)化昏滴，可能并不能明顯的提升模型的效果。
5.2 殘差正態(tài)性檢驗（qqplot）
car包qqplot()函數(shù)提供了精確的正態(tài)假設(shè)檢驗方法对人，置信區(qū)間通過虛線劃定谣殊，當(dāng)絕大多數(shù)點都落在置信區(qū)間時，說明正態(tài)性假設(shè)符合得很好牺弄。
5.3 離群值（強影響點）
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)準過濾最異常值，將該過濾標(biāo)準應(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 誤差的獨立性
car包提供了durbinWatsonTest()函數(shù)蛇捌，用于做Durbin-Watson檢驗抚恒，檢測誤差的序列相關(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络拌，不顯著柑爸，誤差項之間獨立，不存在自相關(guān)性盒音。
5.5 線性相關(guān)性
通過成分殘差圖（component + residual plots）檢查因變量和自變量之間是否呈線性關(guān)系。若圖形存在非線性馅而，則說明可能對預(yù)測變量的函數(shù)形式建模不夠充分祥诽，那么需要添加一些曲線成分，比如多項式瓮恭，對一個或多個變量進行變換（log(x)代替x）雄坪，或用其他回歸變體形式而不是線性回歸。
crPlots(lm5)
5.6 同方差性
判斷方差是否恒定屯蹦，nvcTest()函數(shù)生成一個記分檢驗维哈，原假設(shè)為誤方差不變，備擇假設(shè)為誤差方差隨著擬合值水平的變化而變化登澜。若檢驗顯著阔挠，說明存在誤方差不恒定。
spreadLevelPlot()函數(shù)創(chuàng)建了一個添加了最佳擬合曲線的散點圖脑蠕，展示了標(biāo)準化殘差絕對值與擬合值的關(guān)系购撼。
ncvTest(lm5)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.4044358, Df = 1, p = 0.52481
記分檢驗不顯著：p=0.52481，說明滿足方差不變假設(shè)谴仙，也可以通過分布水平看到這一點迂求，點在水平的最佳擬合曲線周圍呈隨機分布。
spreadLevelPlot(lm5)
5.7 多重共線性（一個自變量不存在）
在回歸分析的時候晃跺，古典模型中有一個基本的假定就是自變量之間是不相關(guān)的揩局，但是如果我們在擬合出來的回歸模型出現(xiàn)了自變量之間高度相關(guān)的話，可能對結(jié)果又產(chǎn)生影響掀虎，我們稱這個問題為多重共線性凌盯。
根據(jù)經(jīng)驗表明，多重共線性存在的一個標(biāo)志就是模型存在較大的標(biāo)準差和較小的T統(tǒng)計量涩盾，如果一個模型的可決系數(shù)R2很大十气，F(xiàn)檢驗高度限制，但偏回歸系數(shù)的T檢驗幾乎都不顯著春霍，那么模型很可能是存在多重共線性砸西。
- 利用自變量之間的簡單相關(guān)系數(shù)檢驗，一般而言，如果兩個解釋變量的簡單相關(guān)系數(shù)較高芹枷，則可以認為是存在著嚴重的多重共線性衅疙。
- 隨著多重共線性的嚴重程度增強，方差膨脹因子會逐漸的變大鸳慈，在回歸中我們用VIF表示方差膨脹因子饱溢，VIF=1/(1-R2)。一般當(dāng)VIF>=10的時候走芋，我們就可以認為存在嚴重多重共線性绩郎。
car包中的vif()函數(shù)可以幫我們算出這個方差膨脹：
library(car)
vif(lm5)
存在多重共線性的解決方法之一：利用MASS包中的函數(shù)lm.ridge()來建立嶺回歸模型。
6翁逞、回歸診斷（gvlma包）
gvlma()函數(shù)肋杖，用于對線性模型假設(shè)進行綜合驗證，同時還能驗證偏斜度挖函，峰度和異方差的評價状植。從輸出項中可以看出，假設(shè)檢驗的顯著性水平是5%怨喘。如果p<0.05津畸，說明違反了假設(shè)條件。從Global Stat 的Decision 欄中必怜，可以看到數(shù)據(jù)滿足OLS回歸模型所有統(tǒng)計假設(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.
