R中的線性回歸函數(shù)比較簡單，就是lm()蚜枢，比較復(fù)雜的是對(duì)線性模型的診斷和調(diào)整皆刺。這里結(jié)合Statistical Learning和杜克大學(xué)的Data Analysis and Statistical Inference的章節(jié)以及《R語言實(shí)戰(zhàn)》的OLS(Ordinary Least Square)回歸模型章節(jié)來總結(jié)一下，診斷多元線性回歸模型的操作分析步驟耕突。
1笤成、選擇預(yù)測變量
因變量比較容易確定，多元回歸模型中難在自變量的選擇眷茁。自變量選擇主要可分為向前選擇（逐次加使RSS最小的自變量）炕泳，向后選擇（逐次扔掉p值最大的變量）。個(gè)人傾向于向后選擇法上祈，一來p值比較直觀培遵，模型返回結(jié)果直接給出了各變量的p值，卻沒有直接給出RSS登刺；二來當(dāng)自變量比較多時(shí)籽腕，一個(gè)個(gè)加比較麻煩。
Call:
lm(formula = Sales ~ . + Income:Advertising + Age:Price, data = Carseats)
Residuals:
Min 1Q Median 3Q Max
-2.9208 -0.7503 0.0177 0.6754 3.3413
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.5755654 1.0087470 6.519 2.22e-10 ***
CompPrice 0.0929371 0.0041183 22.567 < 2e-16 ***
Income 0.0108940 0.0026044 4.183 3.57e-05 ***
Advertising 0.0702462 0.0226091 3.107 0.002030 **
Population 0.0001592 0.0003679 0.433 0.665330
Price -0.1008064 0.0074399 -13.549 < 2e-16 ***
ShelveLocGood 4.8486762 0.1528378 31.724 < 2e-16 ***
ShelveLocMedium 1.9532620 0.1257682 15.531 < 2e-16 ***
Age -0.0579466 0.0159506 -3.633 0.000318 ***
Education -0.0208525 0.0196131 -1.063 0.288361
UrbanYes 0.1401597 0.1124019 1.247 0.213171
USYes -0.1575571 0.1489234 -1.058 0.290729
Income:Advertising 0.0007510 0.0002784 2.698 0.007290 **
Price:Age 0.0001068 0.0001333 0.801 0.423812
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.011 on 386 degrees of freedom
Multiple R-squared: 0.8761, Adjusted R-squared: 0.8719
F-statistic: 210 on 13 and 386 DF, p-value: < 2.2e-16
構(gòu)建一個(gè)回歸模型后纸俭，先看F統(tǒng)計(jì)量的p值皇耗，這是對(duì)整個(gè)模型的假設(shè)檢驗(yàn)，原假設(shè)是各系數(shù)都為0揍很，如果連這個(gè)p值都不顯著郎楼，無法證明至少有一個(gè)自變量對(duì)因變量有顯著性影響矾瘾，這個(gè)模型便不成立。然后看Adjusted R2箭启，每調(diào)整一次模型壕翩，應(yīng)該力使它變大；Adjusted R2越大說明模型中相關(guān)的自變量對(duì)因變量可解釋的變異比例越大傅寡，模型的預(yù)測性就越好放妈。
構(gòu)建了線性模型后，如果是一元線性回歸荐操，可以畫模型圖初步判斷一下線性關(guān)系（多元回歸模型不好可視化）：
par(mfrow=c(1,1))
plot(medv~lstat,Boston)
fit1=lm(medv~lstat,data=Boston)
abline(fit1,col="red")
2芜抒、模型診斷
確定了回歸模型的自變量并初步得到一個(gè)線性回歸模型，并不是直接可以拿來用的托启，還要進(jìn)行驗(yàn)證和診斷宅倒。診斷之前，先回顧多元線性回歸模型的假設(shè)前提（by Data Analysis and Statistical Inference）：
（數(shù)值型）自變量要與因變量有線性關(guān)系屯耸；
殘差基本呈正態(tài)分布拐迁；
殘差方差基本不變（同方差性）；
殘差（樣本）間相關(guān)獨(dú)立疗绣。
一個(gè)好的多元線性回歸模型應(yīng)當(dāng)盡量滿足這4點(diǎn)假設(shè)前提线召。
用lm()構(gòu)造一個(gè)線性模型fit后，plot(fit)即可返回4張圖（可以par(mfrow=c(2,2))一次展現(xiàn)）多矮，這4張圖可作初步檢驗(yàn)：
左上圖用來檢驗(yàn)假設(shè)1缓淹，如果散點(diǎn)看不出什么規(guī)律，則表示線性關(guān)系良好塔逃，若有明顯關(guān)系讯壶，則說明非線性關(guān)系明顯。右上圖用來檢驗(yàn)假設(shè)2湾盗，若散點(diǎn)大致都集中在QQ圖中的直線上伏蚊，則說明殘差正態(tài)性良好。左下圖用來檢驗(yàn)假設(shè)3淹仑，若點(diǎn)在曲線周圍隨機(jī)分布丙挽，則可認(rèn)為假設(shè)3成立；若散點(diǎn)呈明顯規(guī)律匀借，比如方差隨均值而增大，則越往右的散點(diǎn)上下間距會(huì)越大平窘，方差差異就越明顯吓肋。假設(shè)4的獨(dú)立性無法通過這幾張圖來檢驗(yàn)，只能通過數(shù)據(jù)本身的來源的意義去判斷瑰艘。
右下圖是用來檢驗(yàn)異常值是鬼。異常值與三個(gè)概念有關(guān)：
離群點(diǎn)：y遠(yuǎn)離散點(diǎn)主體區(qū)域的點(diǎn)
杠桿點(diǎn)：x遠(yuǎn)離散點(diǎn)主體區(qū)域的點(diǎn)肤舞，一般不影響回歸直線的斜率
強(qiáng)影響點(diǎn)：影響回歸直線的斜率，一般是高杠桿點(diǎn)均蜜。
對(duì)于多元線性回歸李剖，高杠桿點(diǎn)不一定就是極端點(diǎn)，有可能是各個(gè)變量的取值都正常囤耳，但仍然偏離散點(diǎn)主體篙顺。
對(duì)于異常值，可以謹(jǐn)慎地刪除充择，看新的模型是否效果更好德玫。
《R語言實(shí)戰(zhàn)》里推薦了更好的診斷方法，總結(jié)如下椎麦。
1宰僧、多元線性回歸假設(shè)驗(yàn)證：
gvlma包的gvlma()函數(shù)可對(duì)擬合模型的假設(shè)作綜合驗(yàn)證，并對(duì)峰度观挎、偏度進(jìn)行驗(yàn)證琴儿。
states <- as.data.frame(state.x77[, c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income +
Frost, data = states)
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)
Call:
gvlma(x = fit)
Value p-value Decision
Global Stat 2.7728 0.5965 Assumptions acceptable.
Skewness 1.5374 0.2150 Assumptions acceptable.
Kurtosis 0.6376 0.4246 Assumptions acceptable.
Link Function 0.1154 0.7341 Assumptions acceptable.
Heteroscedasticity 0.4824 0.4873 Assumptions acceptable.
最后的Global Stat是對(duì)4個(gè)假設(shè)條件進(jìn)行綜合驗(yàn)證，通過了即表示4個(gè)假設(shè)驗(yàn)證都通過了嘁捷。最后的Heterosceasticity是進(jìn)行異方差檢測凤类。注意這里假設(shè)檢驗(yàn)的原假設(shè)都是假設(shè)成立，所以當(dāng)p>0.05時(shí)普气，假設(shè)才能能過驗(yàn)證谜疤。
如果綜合驗(yàn)證不通過，也有其他方法對(duì)4個(gè)假設(shè)條件分別驗(yàn)證：
線性假設(shè)
library(car)
crPlots(fit)
返回的圖是各個(gè)自變量與殘差（因變量）的線性關(guān)系圖现诀，若存著明顯的非線性關(guān)系夷磕，則需要對(duì)自變量作非線性轉(zhuǎn)化。書中說這張圖表明線性關(guān)系良好仔沿。
正態(tài)性
library(car)
qqPlot(fit,labels = row.names(states),id.method = "identify",simulate = TRUE,main = "Q-Q Plot")
qqPlot()可以生成交互式的qq圖坐桩，選中異常點(diǎn)，就返回該點(diǎn)的名稱封锉。該圖中除了Nevad點(diǎn)绵跷，其他點(diǎn)都在直線附近，可見正態(tài)性良好成福。
同方差性
library(car)
ncvTest(fit)
Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.746514 Df = 1 p = 0.1863156
p值大于0.05碾局，可認(rèn)為滿足方差相同的假設(shè)。
獨(dú)立性
library(car) durbinWatsonTest(fit) lag Autocorrelation D-W Statistic p-value 1? ? ? -0.2006929? ? ? 2.317691? 0.258 Alternative hypothesis: rho != 0
p值大于0.05奴艾，可認(rèn)為誤差之間相互獨(dú)立净当。
除了以上4點(diǎn)基本假設(shè)，還有其他方面需要進(jìn)行診斷——
2、多重共線性
理想中的線性模型各個(gè)自變量應(yīng)該是線性無關(guān)的像啼，若自變量間存在共線性俘闯，則會(huì)降低回歸系數(shù)的準(zhǔn)確性。一般用方差膨脹因子VIF(Variance Inflation Factor)來衡量共線性忽冻，《統(tǒng)計(jì)學(xué)習(xí)》中認(rèn)為VIF超過5或10就存在共線性真朗，《R語言實(shí)戰(zhàn)》中認(rèn)為VIF大于4則存在共線性。理想中的線性模型VIF=1僧诚，表完全不存在共線性遮婶。
library(car)
vif(fit)
Population Illiteracy Income Frost
1.245282 2.165848 1.345822 2.082547
可見這4個(gè)自變量VIF都比較小，可認(rèn)為不存在多重共線性的問題振诬。
3蹭睡、異常值檢驗(yàn)
離群點(diǎn)
離群點(diǎn)有三種判斷方法：一是用qqPlot()畫QQ圖，落在置信區(qū)間（上圖中兩條虛線）外的即可認(rèn)為是離群點(diǎn)赶么，如上圖中的Nevad點(diǎn)肩豁；一種是判斷學(xué)生標(biāo)準(zhǔn)化殘差值，絕對(duì)值大于2（《R語言實(shí)戰(zhàn)》中認(rèn)為2辫呻，《統(tǒng)計(jì)學(xué)習(xí)》中認(rèn)為3）的可認(rèn)為是離群點(diǎn)清钥。
plot(x=fitted(fit),y=rstudent(fit))
abline(h=3,col="red",lty=2)
abline(h=-3,col="red",lty=2)
which(abs(rstudent(fit))>3)
Nevada
28
還有一種方法是利用car包里的outlierTest()函數(shù)進(jìn)行假設(shè)檢驗(yàn)：
library(car)
outlierTest(fit)
rstudent unadjusted p-value Bonferonni p
Nevada 3.542929 0.00095088 0.047544
這個(gè)函數(shù)用來檢驗(yàn)最大的標(biāo)準(zhǔn)化殘差值，如果p>0.05放闺，可以認(rèn)為沒有離群點(diǎn)祟昭；若p<0.05，則該點(diǎn)是離群點(diǎn)怖侦，但不能說明只有一個(gè)離群點(diǎn)篡悟，可以把這個(gè)點(diǎn)刪除之后再作檢驗(yàn)。第三種方法可以與第二種方法結(jié)合起來使用匾寝。
高杠桿點(diǎn)
高杠桿值觀測點(diǎn)搬葬，即是與其他預(yù)測變量有關(guān)的離群點(diǎn)。換句話說艳悔，它們是由許多異常的預(yù)測變量值組合起來的急凰，與響應(yīng)變量值沒有關(guān)系〔履辏《統(tǒng)計(jì)學(xué)習(xí)》中給出了一個(gè)杠桿統(tǒng)計(jì)量抡锈，《R語言實(shí)戰(zhàn)》中給出了一種具體的操作方法。（兩本書也稍有出入乔外，《統(tǒng)計(jì)學(xué)習(xí)》中平均杠桿值為(p+1)/n床三，而在《R語言實(shí)戰(zhàn)》中平均杠桿值為p/n；事實(shí)上在樣本量n比較大時(shí)袁稽，幾乎沒有差別勿璃。）
hat.plot <- function(fit){
p <- length(coefficients(fit))
n <- length(fitted(fit))
plot(hatvalues(fit),main = "Index Plot of Hat Values")
abline(h=c(2,3)*p/n,col="red",lty=2)
identify(1:n, hatvalues(fit), names(hatvalues(fit))) #這句產(chǎn)生交互效果，選中某個(gè)點(diǎn)后推汽，關(guān)閉后返回點(diǎn)的名稱
}
hat.plot(fit)
超過2倍或3倍的平均杠桿值即可認(rèn)為是高杠桿點(diǎn)补疑，這里把Alaska和California作為高杠桿點(diǎn)。
強(qiáng)影響點(diǎn)
強(qiáng)影響點(diǎn)是那種若刪除則模型的系數(shù)會(huì)產(chǎn)生明顯的變化的點(diǎn)歹撒。一種方法是計(jì)算Cook距離莲组，一般來說， Cook’s D值大于4/(n-k -1)暖夭，則表明它是強(qiáng)影響點(diǎn)锹杈，其中n 為樣本量大小， k 是預(yù)測變量數(shù)目迈着。
cutoff <- 4/(nrow(states-length(fit$coefficients)-2)) #coefficients加上了截距項(xiàng)竭望，因此要多減1
plot(fit,which=4,cook.levels = cutoff)
abline(h=cutoff,lty=2,col="red")
實(shí)際上這就是前面診斷的4張圖之一，語句還是plot(fit)裕菠，which=4表示指定第4張圖咬清，cook.levels可設(shè)定標(biāo)準(zhǔn)值。紅色虛線以上就返回了強(qiáng)影響點(diǎn)奴潘。
car包里的influencePlot()函數(shù)能一次性同時(shí)檢查離群點(diǎn)旧烧、高杠桿點(diǎn)、強(qiáng)影響點(diǎn)画髓。
library(car)?
influencePlot(fit,id.method = "identity", main="Influence Plot",
sub="Circle size is proportional to Cook's distance")
StudRes Hat CookD
Alaska 1.753692 0.43247319 0.4480510
Nevada 3.542929 0.09508977 0.2099157
縱坐標(biāo)超過+2或小于-2的點(diǎn)可被認(rèn)為是離群點(diǎn)掘剪，水平軸超過0.2或0.3的州有高杠桿值（通常為預(yù)測值的組合）。圓圈大小與影響成比例奈虾，圓圈很大的點(diǎn)可能是對(duì)模型參數(shù)的估計(jì)造成的不成比例影響的強(qiáng)影響點(diǎn)夺谁。
3、模型調(diào)整
到目前為止肉微，《統(tǒng)計(jì)學(xué)習(xí)》中提到的多元線性回歸模型潛在的問題匾鸥，包括4個(gè)假設(shè)不成立、異常值浪册、共線性的診斷方法在上面已經(jīng)全部得到解決扫腺。這里總結(jié)、延伸《R語言實(shí)戰(zhàn)》里提到的調(diào)整方法——
刪除觀測點(diǎn)
對(duì)于異常值一個(gè)最簡單粗暴又高效的方法就是直接刪除村象，不過有兩點(diǎn)要注意笆环。一是當(dāng)數(shù)據(jù)量大的時(shí)候可以這么做，若數(shù)據(jù)量較小則應(yīng)慎重厚者；二是根據(jù)數(shù)據(jù)的意義判斷躁劣，若明顯就是錯(cuò)誤就可以直接刪除，否則需判斷是否會(huì)隱藏著深層的現(xiàn)象库菲。
另外刪除觀測點(diǎn)后要與刪除之前的模型作比較账忘，看模型是否變得更好。
變量變換
?
在進(jìn)行非線性變換之前，先看看4個(gè)假設(shè)是否成立鳖擒，如果成立可以不用變換溉浙；沒必要追求更好的擬合效果而把模型搞得太復(fù)雜，這有可能出現(xiàn)過擬合現(xiàn)象蒋荚。如果連假設(shè)檢驗(yàn)都不通過戳稽，可以通過變量變換來調(diào)整模型。這里只討論線性關(guān)系不佳的情況期升，其他情況遇到了再說惊奇。
（1）多項(xiàng)式回歸
如果殘差圖中呈現(xiàn)明顯的非線性關(guān)系，可以考慮對(duì)自變量進(jìn)行多項(xiàng)式回歸播赁。舉一個(gè)例子：
library(MASS)
library(ISLR)
fit1=lm(medv~lstat,data=Boston)
plot(fit1,which = 1)
可以看到這個(gè)一元線性回歸模型的殘差圖中颂郎，散點(diǎn)的規(guī)律還是比較明顯，說明線性關(guān)系較弱容为。
fit1_1 <- lm(medv~poly(lstat,2),data = Boston)
plot(fit1_1,which = 1)
將自變量進(jìn)行2次多項(xiàng)式回歸后乓序，發(fā)現(xiàn)現(xiàn)在的殘差圖好多了，散點(diǎn)基本無規(guī)律舟奠，線性關(guān)系較明顯竭缝。
再看看兩個(gè)模型的整體效果——
summary(fit1)
summary(fit1_1)
可見多項(xiàng)式回歸的模型Adjusted R2也增大了，模型的解釋性也變強(qiáng)了沼瘫。
多項(xiàng)式回歸在《統(tǒng)計(jì)學(xué)習(xí)》后面的非線性模型中還會(huì)提到抬纸，到時(shí)候再討論。
（2）Box-Tidwell變換
car包中的boxTidwell() 函數(shù)通過獲得預(yù)測變量冪數(shù)的最大似然估計(jì)來改善線性關(guān)系耿戚。
library(car)
boxTidwell(Murder~Population+Illiteracy,data=states)
Score Statistic p-value MLE of lambda
Population -0.3228003 0.7468465 0.8693882
Illiteracy 0.6193814 0.5356651 1.3581188
iterations = 19
#這里看lambda湿故，表示各個(gè)變量的冪次數(shù)
lmfit <- lm(Murder~Population+Illiteracy,data=states)
lmfit2 <- lm(Murder~I(Population^0.87)+I(Illiteracy^1.36),data=states)
plot(lmfit,which = 1)
plot(lmfit2,which = 1)
summary(lmfit)
summary(lmfit2)
可以發(fā)現(xiàn)殘差圖和Adjusted R2的提升都甚微，因此沒有必要作非線性轉(zhuǎn)換膜蛔。
4坛猪、模型分析
（1）模型比較
前面只是簡單得用Adjusted R2來比較模型，《R語言實(shí)戰(zhàn)》里介紹了可以用方差分析來比較嵌套模型（即它的一些項(xiàng)完全包含在另一個(gè)模型中）有沒有顯著性差異皂股。方差分析的思想是：如果線性模型y~x1+x2+x3與y~x1+x2沒有顯著性差異墅茉，若同時(shí)x3變量對(duì)模型也不顯著，那就沒必要加上變量x3呜呐。下面進(jìn)行試驗(yàn)：
aovfit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
aovfit2 <- lm(Murder~Population+Illiteracy,data=states)
anova(aovfit1,aovfit2)
Analysis of Variance Table
Model 1: Murder ~ Population + Illiteracy + Income + Frost
Model 2: Murder ~ Population + Illiteracy
Res.Df RSS Df Sum of Sq F Pr(>F)
1 45 289.17
2 47 289.25 -2 -0.078505 0.0061 0.9939
summary(aovfit1)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.235e+00 3.866e+00 0.319 0.7510
Population 2.237e-04 9.052e-05 2.471 0.0173 *
Illiteracy 4.143e+00 8.744e-01 4.738 2.19e-05 ***
Income 6.442e-05 6.837e-04 0.094 0.9253
Frost 5.813e-04 1.005e-02 0.058 0.9541
Residual standard error: 2.535 on 45 degrees of freedom
Multiple R-squared: 0.567, Adjusted R-squared: 0.5285
F-statistic: 14.73 on 4 and 45 DF, p-value: 9.133e-08
summary(aovfit2)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.652e+00 8.101e-01 2.039 0.04713 *
Population 2.242e-04 7.984e-05 2.808 0.00724 **
Illiteracy 4.081e+00 5.848e-01 6.978 8.83e-09 ***
Residual standard error: 2.481 on 47 degrees of freedom
Multiple R-squared: 0.5668, Adjusted R-squared: 0.5484
F-statistic: 30.75 on 2 and 47 DF, p-value: 2.893e-09
Income和Frost兩個(gè)變量不顯著就斤，兩個(gè)模型之間沒有顯著性差異，就可以不加這兩個(gè)變量蘑辑。刪去這兩個(gè)不顯著的變量后洋机，R2略微減少，Adjusted R2增大洋魂，這也符合二者的定義绷旗。
《R語言實(shí)戰(zhàn)》里還介紹到了用AIC(Akaike Information Criterion喜鼓，赤池信息準(zhǔn)則)值來比較模型，AIC值越小的模型優(yōu)先選擇衔肢，原理不明庄岖。
aovfit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
aovfit2 <- lm(Murder~Population+Illiteracy,data=states)
AIC(aovfit1,aovfit2)
df AIC
aovfit1 6 241.6429
aovfit2 4 237.6565
第二個(gè)模型AIC值更小，因此選第二個(gè)模型（真是簡單粗暴）膀懈。注：ANOVA需限定嵌套模型顿锰，AIC則不需要谨垃∑袈В可見AIC是更簡單也更實(shí)用的模型比較方法。
（2）變量選擇
這里的變量選擇與最開始的變量選擇同也不同刘陶，雖然是一回事胳赌，但一開始是一個(gè)粗略的變量的選擇，主要是為了構(gòu)建模型匙隔；這里則要進(jìn)行細(xì)致的變量選擇來調(diào)整模型疑苫。
逐步回歸
前面提到的向前或向后選擇或者是同時(shí)向前向后選擇變量都是逐步回歸法。MASS包中的stepAIC() 函數(shù)可以實(shí)現(xiàn)逐步回歸模型（向前纷责、向后和向前向后）捍掺，依據(jù)的是精確AIC準(zhǔn)則。以下實(shí)例是向后回歸法：
library(MASS)
aovfit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data=states)
stepAIC(aovfit1,direction = "backward") # “forward”為向前選擇再膳，"backward"為向后選擇挺勿，"both"為混合選擇
Start: AIC=97.75
Murder ~ Population + Illiteracy + Income + Frost
Df Sum of Sq RSS AIC
- Frost 1 0.021 289.19 95.753
- Income 1 0.057 289.22 95.759
289.17 97.749
- Population 1 39.238 328.41 102.111
- Illiteracy 1 144.264 433.43 115.986
Step: AIC=95.75
Murder ~ Population + Illiteracy + Income
Df Sum of Sq RSS AIC
- Income 1 0.057 289.25 93.763
289.19 95.753
- Population 1 43.658 332.85 100.783
- Illiteracy 1 236.196 525.38 123.605
Step: AIC=93.76
Murder ~ Population + Illiteracy
Df Sum of Sq RSS AIC
289.25 93.763
- Population 1 48.517 337.76 99.516
- Illiteracy 1 299.646 588.89 127.311
Call:
lm(formula = Murder ~ Population + Illiteracy, data = states)
Coefficients:
(Intercept) Population Illiteracy
1.6515497 0.0002242 4.0807366
可見原本的4元回歸模型向后退了兩次，最終穩(wěn)定成了2元回歸模型喂柒，與前面模型比較的結(jié)果一致不瓶。
全子集回歸
《R語言實(shí)戰(zhàn)》里提到了逐步回歸法的局限：不是每個(gè)模型都評(píng)價(jià)了，不能保證選擇的是“最佳”模型灾杰。比如上例中蚊丐，從Murder ~ Population + Illiteracy + Income + Frost到Murder ~ Population + Illiteracy + Income再到Murder~Population+Illiteracy雖然AIC值確實(shí)在減少，但Murder ~ Population + Illiteracy + Frost沒評(píng)價(jià)艳吠，如果遇到變量很多的情況下麦备，逐步回歸只沿一個(gè)方向回歸，就有可能錯(cuò)過最優(yōu)的回歸方向昭娩。
library(leaps)
leaps <- regsubsets(Murder~Population+Illiteracy+Income+Frost,data=states,nbest=4)
plot(leaps,scale = "adjr2")
橫坐標(biāo)是變量凛篙，縱坐標(biāo)是Adjusted R2，可見除截距項(xiàng)以外题禀，只選定Population和Illiteracy這兩個(gè)變量鞋诗，可以使線性模型有最大的Adjusted R2。
全子集回歸比逐步回歸范圍更廣迈嘹，模型優(yōu)化效果更好黍翎，但是一旦變量數(shù)多了之后，全子集回歸迭代的次數(shù)就很多材原，就會(huì)很慢显歧。
事實(shí)上，變量的選擇不是機(jī)械式地只看那幾個(gè)統(tǒng)計(jì)指標(biāo)嗡午，更主要的是根據(jù)數(shù)據(jù)的實(shí)際意義，從業(yè)務(wù)角度上來選擇合適的變量。
線性模型變量的選擇在《統(tǒng)計(jì)學(xué)習(xí)》后面的第6章還會(huì)繼續(xù)講到覆劈，到時(shí)繼續(xù)綜合討論。
（3）交互項(xiàng)
交互項(xiàng)《統(tǒng)計(jì)學(xué)習(xí)》中花了一定篇幅來描寫沛励，但在《R語言實(shí)戰(zhàn)》是在方差分析章節(jié)中討論责语。添加變量間的交互項(xiàng)有時(shí)可以改善線性關(guān)系，提高Adjusted R2目派。針對(duì)數(shù)據(jù)的實(shí)際意義坤候，如果兩個(gè)基本上是獨(dú)立的，也很難產(chǎn)生交互企蹭、產(chǎn)生協(xié)同效應(yīng)的變量白筹，那就不必考慮交互項(xiàng)；只有從業(yè)務(wù)角度分析谅摄，有可能產(chǎn)生協(xié)同效應(yīng)的變量間才考慮交互項(xiàng)徒河。
涉及到交互項(xiàng)有一個(gè)原則：如果交互項(xiàng)是顯著的，那么即使變量不顯著送漠，也要放在回歸模型中顽照；若變量和交互項(xiàng)都不顯著，則可以都不放螺男。
（4）交叉驗(yàn)證
Andrew Ng的Machine Learning中就提到了棒厘，模型對(duì)舊數(shù)據(jù)擬合得好不一定就對(duì)新數(shù)據(jù)預(yù)測得好。因此一個(gè)數(shù)據(jù)集應(yīng)當(dāng)被分兩訓(xùn)練集和測試集兩部分（或者訓(xùn)練集下隧、交叉驗(yàn)證集奢人、測試集三部分），訓(xùn)練好的模型還要在新數(shù)據(jù)中測試性能淆院。
所謂交叉驗(yàn)證何乎，即將一定比例的數(shù)據(jù)挑選出來作為訓(xùn)練樣本，另外的樣本作保留樣本土辩，先在訓(xùn)練樣本上獲取回歸方程支救，然后在保留樣本上做預(yù)測。由于保留樣本不涉及模型參數(shù)的選擇拷淘，該樣本可獲得比新數(shù)據(jù)更為精確的估計(jì)各墨。
在k 重交叉驗(yàn)證中，樣本被分為k個(gè)子樣本启涯，輪流將k-1個(gè)子樣本組合作為訓(xùn)練集贬堵，另外1個(gè)子樣本作為保留集恃轩。這樣會(huì)獲得k 個(gè)預(yù)測方程，記錄k 個(gè)保留樣本的預(yù)測表現(xiàn)結(jié)果黎做，然后求其平均值叉跛。
bootstrap包中的crossval()函數(shù)可以實(shí)現(xiàn)k重交叉驗(yàn)證。
shrinkage <- function(fit, k = 10) {
require(bootstrap)
# define functions
theta.fit <- function(x, y) {
lsfit(x, y)
}
theta.predict <- function(fit, x) {
cbind(1, x) %*% fit$coef
}
# matrix of predictors
x <- fit$model[, 2:ncol(fit$model)]
# vector of predicted values
y <- fit$model[, 1]
results <- crossval(x, y, theta.fit, theta.predict, ngroup = k)
r2 <- cor(y, fit$fitted.values)^2
r2cv <- cor(y, results$cv.fit)^2
cat("Original R-square =", r2, "\n")
cat(k, "Fold Cross-Validated R-square =", r2cv, "\n")
cat("Change =", r2 - r2cv, "\n")
}
這個(gè)自定義的shrinkage()函數(shù)用來做k重交叉驗(yàn)證蒸殿，比計(jì)算訓(xùn)練集和交叉驗(yàn)證集的R方差異筷厘。這個(gè)函數(shù)里涉及到一個(gè)概念：復(fù)相關(guān)系數(shù)。復(fù)相關(guān)系數(shù)實(shí)際上就是y和fitted(y)的簡單相關(guān)系數(shù)宏所。對(duì)于一元線性回歸酥艳，R2就是簡單相關(guān)系數(shù)的平方；對(duì)于多元線性回歸楣铁，R2是復(fù)相關(guān)系數(shù)的平方玖雁。這個(gè)我沒有成功地從公式上推導(dǎo)證明成立，就記下吧盖腕。這個(gè)方法用到了自助法的思想，這個(gè)在統(tǒng)計(jì)學(xué)習(xí)后面會(huì)細(xì)致講到浓镜。
fit <- lm(Murder ~ Population + Income + Illiteracy +
Frost, data = states)
shrinkage(fit)
Original R-square = 0.5669502
10 Fold Cross-Validated R-square = 0.441954
Change = 0.1249963
可見這個(gè)4元回歸模型在交叉驗(yàn)證集中的R2下降了0.12之多溃列。若換成前面分析的2元回歸模型——
fit2 <- lm(Murder ~ Population + Illiteracy , data = states)
shrinkage(fit2)
Original R-square = 0.5668327
10 Fold Cross-Validated R-square = 0.517304
Change = 0.04952868
這次R2下降只有約0.05。R2減少得越少膛薛，則預(yù)測得越準(zhǔn)確听隐。
5、模型應(yīng)用
（1）預(yù)測
最重要的應(yīng)用毫無疑問就是用建立的模型進(jìn)行預(yù)測了哄啄。構(gòu)建好模型后雅任，可用predict()函數(shù)進(jìn)行預(yù)測——
fit2 <- lm(Murder ~ Population + Illiteracy , data = states)
predict(fit2,
newdata = data.frame(Population=c(2000,3000),Illiteracy=c(1.7,2.2)),
interval = "confidence")
fit lwr upr
1 9.037174 8.004911 10.06944
2 11.301729 9.866851 12.73661
這里newdata提供了兩個(gè)全新的點(diǎn)供模型來預(yù)測。還可以用interval指定返回置信區(qū)間（confidence）或者預(yù)測區(qū)間（prediction）咨跌，這也反映了統(tǒng)計(jì)與機(jī)器學(xué)習(xí)的一個(gè)差異——可解釋性沪么。注意置信區(qū)間考慮的是平均值，而預(yù)測區(qū)間考慮的是單個(gè)觀測值锌半，所以預(yù)測區(qū)間永遠(yuǎn)比置信區(qū)間廣禽车，因此預(yù)測區(qū)間考慮了單個(gè)觀測值的不可約誤差；而平均值同時(shí)也把不可約誤差給抵消掉了刊殉。
（2）相對(duì)重要性
有的時(shí)候需要解釋模型中各個(gè)自變量對(duì)因變量的重要程度殉摔，簡單處理可以直接看系數(shù)即可，《R語言實(shí)戰(zhàn)》里自定義了一個(gè)relweights()函數(shù)可以計(jì)算各個(gè)變量的權(quán)重：
relweights <- function(fit, ...) {
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar, 2:nvar]
rxy <- R[2:nvar, 1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
# correlations between original predictors and new orthogonal variables
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda^2
# regression coefficients of Y on orthogonal variables
beta <- solve(lambda) %*% rxy
rsquare <- colSums(beta^2)
rawwgt <- lambdasq %*% beta^2
import <- (rawwgt/rsquare) * 100
lbls <- names(fit$model[2:nvar])
rownames(import) <- lbls
colnames(import) <- "Weights"
# plot results
barplot(t(import), names.arg = lbls, ylab = "% of R-Square",
xlab = "Predictor Variables", main = "Relative Importance of Predictor Variables",
sub = paste("R-Square = ", round(rsquare, digits = 3)),
...)
return(import)
}
不要在意算法原理和代碼邏輯這種細(xì)節(jié)记焊，直接看結(jié)果：
fit <- lm(Murder ~ Population + Illiteracy + Income +
Frost, data = states)
relweights(fit, col = "lightgrey")
Weights
Population 14.723401
Illiteracy 59.000195
Income 5.488962
Frost 20.787442
fit$coefficients
(Intercept) Population Illiteracy Income Frost
1.2345634112 0.0002236754 4.1428365903 0.0000644247 0.0005813055
barplot(fit$coefficients[2:5])
?
在本例中逸月，相對(duì)權(quán)重與系數(shù)的排序結(jié)果一致。推薦用相對(duì)權(quán)重遍膜。
?
?
?
?
?
?
