20180404(從有道遷移)
回歸
-
回歸的多面性
回歸分析的各種變體
回歸類型 用 途 簡(jiǎn)單線性 用一個(gè)量化的解釋變量預(yù)測(cè)一個(gè)量化的響應(yīng)變量 多項(xiàng)式 用一個(gè)量化的解釋變量預(yù)測(cè)一個(gè)量化的響應(yīng)變量叉庐,模型的關(guān)系是n階多項(xiàng)式 多元線性 用兩個(gè)或多個(gè)量化的解釋變量預(yù)測(cè)一個(gè)量化的響應(yīng)變量 多變量 用一個(gè)或多個(gè)解釋變量預(yù)測(cè)多個(gè)響應(yīng)變量 Logistic 用一個(gè)或多個(gè)解釋變量預(yù)測(cè)一個(gè)類別型響應(yīng)變量 泊松 用一個(gè)或多個(gè)解釋變量預(yù)測(cè)一個(gè)代表頻數(shù)的響應(yīng)變量 Cox比例風(fēng)險(xiǎn) 用一個(gè)或多個(gè)解釋變量預(yù)測(cè)一個(gè)事件(死亡、失敗或舊病復(fù)發(fā))發(fā)生的時(shí)間 時(shí)間序列 對(duì)誤差項(xiàng)相關(guān)的時(shí)間序列數(shù)據(jù)建模 非線性 用一個(gè)或多個(gè)量化的解釋變量預(yù)測(cè)一個(gè)量化的響應(yīng)變量瑞凑,不過(guò)模型是非線性的 非參數(shù) 用一個(gè)或多個(gè)量化的解釋變量預(yù)測(cè)一個(gè)量化的響應(yīng)變量末秃,模型的形式源自數(shù)據(jù)形式,不事先設(shè)定 穩(wěn)健 用一個(gè)或多個(gè)量化的解釋變量預(yù)測(cè)一個(gè)量化的響應(yīng)變量籽御,能抵御強(qiáng)影響點(diǎn)的干擾 -
OLS 回歸
-
普通最小二乘(OLS)回歸法的適用情境
OLS回歸是通過(guò)預(yù)測(cè)變量的加權(quán)和來(lái)預(yù)測(cè)量化的因變量练慕,其中權(quán)重是通過(guò)數(shù)據(jù)估計(jì)而得的參數(shù)
-
OLS回歸擬合模型的形式,其中惰匙,n 為觀測(cè)的數(shù)目,k 為預(yù)測(cè)變量的數(shù)目铃将。
\widehat{y}_{i}=\widehat{\beta}_{0}+\widehat{\beta}_{1}X_{1i}+...+\widehat{\beta}_{k}X_{ki} i=1...n
image 為了能夠恰當(dāng)?shù)亟忉孫LS模型的系數(shù)项鬼,數(shù)據(jù)必須滿足以下統(tǒng)計(jì)假設(shè)。
正態(tài)性 對(duì)于固定的自變量值劲阎,因變量值成正態(tài)分布绘盟。
獨(dú)立性 Yi值之間相互獨(dú)立。
線性 因變量與自變量之間為線性相關(guān)悯仙。
-
同方差性 因變量的方差不隨自變量的水平不同而變化龄毡。也可稱作不變方差,但是說(shuō)同方差性感覺(jué)上更犀利。
如果違背了以上假設(shè)茴她,你的統(tǒng)計(jì)顯著性檢驗(yàn)結(jié)果和所得的置信區(qū)間很可能就不精確。注意己沛,OLS回歸還假定自變量是固定的且測(cè)量無(wú)誤差申尼,但在實(shí)踐中通常都放松了這個(gè)假設(shè)
-
用lm()擬合回歸模型
-
在R中,擬合線性模型最基本的函數(shù)就是lm()霹粥,格式為:
myfit <- lm(formula,data)
,其中:- formula指要擬合的模型形式
- data是一個(gè)數(shù)據(jù)框,包含了用于擬合模型的數(shù)據(jù)浩淘。
- 結(jié)果對(duì)象(本例中是myfit)存儲(chǔ)在一個(gè)列表中男旗,包含了所擬合模型的大量信息。
-
表達(dá)式(formula)形式如下:
y ~ x1+x2+x3+...+xk
- ~左邊為響應(yīng)變量
- 右邊為各個(gè)預(yù)測(cè)變量什荣,預(yù)測(cè)變量之間用+符號(hào)分隔
-
R表達(dá)式中常用的符號(hào)
符號(hào) 用 途 ~ 分隔符號(hào),左邊為響應(yīng)變量桅锄,右邊為解釋變量友瘤。</br>例如,要通過(guò)x盟戏、z和w預(yù)測(cè)y,代碼為y ~ x + z + w + 分隔預(yù)測(cè)變量 : 表示預(yù)測(cè)變量的交互項(xiàng)笛求。</br>例如,要通過(guò)x蜂嗽、z及x與z的交互項(xiàng)預(yù)測(cè)y辱揭,代碼為y ~ x + z + x:z * 表示所有可能交互項(xiàng)的簡(jiǎn)潔方式问窃。</br>代碼y~ x * z * w可展開(kāi)為y ~ x + z + w + x:z + x:w + z:w +x:z:w ^ 表示交互項(xiàng)達(dá)到某個(gè)次數(shù)。</br>代碼y ~ (x + z + w)^2可展開(kāi)為y ~ x + z + w + x:z + x:w + z:w . 表示包含除因變量外的所有變量听皿。</br>例如,若一個(gè)數(shù)據(jù)框包含變量x又厉、y、z和w,代碼y ~ .可展開(kāi)為y ~ x +z + w - 減號(hào)婆排,表示從等式中移除某個(gè)變量。</br>例如赞枕,y ~ (x + z + w)^2 – x:w可展開(kāi)為y ~ x + z + w + x:z + z:w -1 刪除截距項(xiàng)炕婶。</br>例如,表達(dá)式y(tǒng) ~ x - 1擬合y在x上的回歸涯贞,并強(qiáng)制直線通過(guò)原點(diǎn) I() 從算術(shù)的角度來(lái)解釋括號(hào)中的元素宋渔。</br>例如严蓖,y ~ x + (z + w)^2將展開(kāi)為y ~ x + z + w + z:w谈飒。</br>相反, 代碼y~ x + I((z + w)^2)將展開(kāi)為y ~ x + h, h是一個(gè)由z和w的平方和創(chuàng)建的新變量 function 可以在表達(dá)式中用的數(shù)學(xué)函數(shù)。</br>例如手素,log(y) ~ x + z + w表示通過(guò)x、z和w來(lái)預(yù)測(cè)log(y) -
對(duì)擬合線性模型非常有用的其他函數(shù)
函 數(shù) 用 途 summary() 展示擬合模型的詳細(xì)結(jié)果 coefficients() 列出擬合模型的模型參數(shù)(截距項(xiàng)和斜率) confint() 提供模型參數(shù)的置信區(qū)間(默認(rèn)95%) fitted() 列出擬合模型的預(yù)測(cè)值 residuals() 列出擬合模型的殘差值 anova() 生成一個(gè)擬合模型的方差分析表崩哩,或者比較兩個(gè)或更多擬合模型的方差分析表 vcov() 列出模型參數(shù)的協(xié)方差矩陣 AIC() 輸出赤池信息統(tǒng)計(jì)量 plot() 生成評(píng)價(jià)擬合模型的診斷圖 predict() 用擬合模型對(duì)新的數(shù)據(jù)集預(yù)測(cè)響應(yīng)變量值
-
-
簡(jiǎn)單線性回歸:當(dāng)回歸模型包含一個(gè)因變量和一個(gè)自變量時(shí),我們稱為簡(jiǎn)單線性回歸
- 示例:基礎(chǔ)安裝中的數(shù)據(jù)集women提供了15個(gè)年齡在30~39歲間女性的身高和體重信息汹押,我們想通過(guò)身高來(lái)預(yù)測(cè)體重,獲得一個(gè)等式可以幫助我們分辨出那些過(guò)重或過(guò)瘦的個(gè)體
- 代碼:
## 擬合回歸模型 fit <- lm(weight ~ height, data=women) ## 查看擬合后的詳細(xì)結(jié)果 summary(fit) Call: lm(formula = weight ~ height, data = women) Residuals: Min 1Q Median 3Q Max -1.7333 -1.1333 -0.3833 0.7417 3.1167 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -87.51667 5.93694 -14.74 1.71e-09 *** height 3.45000 0.09114 37.85 1.09e-14 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.525 on 13 degrees of freedom Multiple R-squared: 0.991, Adjusted R-squared: 0.9903 F-statistic: 1433 on 1 and 13 DF, p-value: 1.091e-14 ## 在Pr(>|t|)欄妙痹,可以看到回歸系數(shù)(3.45)顯著不為0(p<0.001),表明身高每增高1英寸利赋,體重將預(yù)期增加3.45磅 ## R平方項(xiàng)(0.991)表明模型可以解釋體重99.1%的方差媚送,它也是實(shí)際和預(yù)測(cè)值之間的相關(guān)系數(shù)(R2 = r2?Y)。 ## 殘差標(biāo)準(zhǔn)誤(1.53 lbs)則可認(rèn)為是模型用身高預(yù)測(cè)體重的平均誤差吟秩。 ## F統(tǒng)計(jì)量檢驗(yàn)所有的預(yù)測(cè)變量預(yù)測(cè)響應(yīng)變量是否都在某個(gè)幾率水平之上。由于簡(jiǎn)單回歸只有一個(gè)預(yù)測(cè)變量壮池,此處F檢驗(yàn)等同于身高回歸系數(shù)的t檢驗(yàn) women$weight ## 列出擬合模型的預(yù)測(cè)值 fitted(fit) ## 列出擬合模型的殘差值 residuals(fit) ## 提供模型參數(shù)的置信區(qū)間(默認(rèn)95%) confint(fit) plot(women$height,women$weight, main="Women Age 30-39", xlab="Height (in inches)", ylab="Weight (in pounds)") # add the line of best fit abline(fit)
-
多項(xiàng)式回歸,椰憋。當(dāng)只有一個(gè)預(yù)測(cè)變量橙依,但同時(shí)包含變量的冪(比如,X、X2浪读、X3)時(shí),我們稱之為多項(xiàng)式回歸
針對(duì)上訴身高問(wèn)題痘拆,你可以通過(guò)添加一個(gè)二次項(xiàng)(即
X2)來(lái)提高回歸的預(yù)測(cè)精度吐葵,因此可以擬合含二次項(xiàng)的等式:fit2 <- lm(weight ~ height + I(height^2),data=women)
-
代碼
fit2 <- lm(weight ~ height + I(height^2), data=women) summary(fit2) Call: lm(formula = weight ~ height + I(height^2), data = women) Residuals: Min 1Q Median 3Q Max -0.50941 -0.29611 -0.00941 0.28615 0.59706 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 261.87818 25.19677 10.393 2.36e-07 *** height -7.34832 0.77769 -9.449 6.58e-07 *** I(height^2) 0.08306 0.00598 13.891 9.32e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.3841 on 12 degrees of freedom Multiple R-squared: 0.9995, Adjusted R-squared: 0.9994 F-statistic: 1.139e+04 on 2 and 12 DF, p-value: < 2.2e-16 ## 在p<0.001水平下,回歸系數(shù)都非常顯著凤藏。模型的方差解釋率已經(jīng)增加到了99.9%欠雌。 ## 二次項(xiàng)的顯著性(t = 13.89检号,p<0.001)表明包含二次項(xiàng)提高了模型的擬合度齐苛。 plot(women$height,women$weight, main="Women Age 30-39", xlab="Height (in inches)", ylab="Weight (in lbs)") lines(women$height,fitted(fit2))
note
image- 在car包中的scatterplot()函數(shù),它可以很容易、方便地繪制二元關(guān)系圖,==這個(gè)功能加強(qiáng)的圖形擂煞,既提供了身高與體重的散點(diǎn)圖、線性擬合曲線和平滑擬合(loess)曲線蒿涎,還在相應(yīng)邊界展示了每個(gè)變量的箱線圖==
- spread=FALSE選項(xiàng)刪除了殘差正負(fù)均方根在平滑曲線上的展開(kāi)和非對(duì)稱信息。
- lty.smooth=2選項(xiàng)設(shè)置loess擬合曲線為虛線嗽冒。
- pch=19選項(xiàng)設(shè)置點(diǎn)為實(shí)心圓(默認(rèn)為空心圓)
library(car) scatterplot(weight ~ height, data=women, spread=FALSE, smoother.args=list(lty=2), pch=19, main="Women Age 30-39", xlab="Height (inches)", ylab="Weight (lbs.)")
-
多元線性回歸,當(dāng)有不止一個(gè)預(yù)測(cè)變量時(shí)干像,則稱為多元線性回歸
多元回歸分析中帅腌,第一步最好檢查一下變量間的相關(guān)性。cor()函數(shù)提供了二變量之間的相關(guān)系數(shù)麻汰,car包中scatterplotMatrix()函數(shù)則會(huì)生成散點(diǎn)圖矩陣
-
使用lm()函數(shù)擬合多元線性回歸模型
- 當(dāng)預(yù)測(cè)變量不止一個(gè)時(shí)速客,回歸系數(shù)的含義為,一個(gè)預(yù)測(cè)變量增加一個(gè)單位五鲫,其他預(yù)測(cè)變量保持不變時(shí)溺职,因變量將要增加的數(shù)量。
cor(states) library(car) scatterplotMatrix(states,spread = FLASE,smooth = 2,main="Scatter Plot Matrix") fit <- lm(Murder ~ Population+Illiteracy+Income+Frost,data=states) summary(fit)
-
有交互項(xiàng)的多元線性回歸
在初步判斷多個(gè)變量之間存在互相作用時(shí),可以將變量之間的交互項(xiàng)加入擬合模型中规婆,通過(guò)對(duì)交互項(xiàng)在模型中的pvalue驗(yàn)證是否相關(guān)性顯著
-
可以通過(guò)effects包中的effect()函數(shù),你可以用圖形展示交互項(xiàng)的結(jié)果
- 格式:plot(effect(term,mod,xlevels),multiple=TRUE)
- term即模型要畫(huà)的項(xiàng),mod為通過(guò)lm()擬合的模型唆姐,xlevels是一個(gè)列表,指定變量要設(shè)定的常量值伸蚯,multiline=TRUE選項(xiàng)表示添加相應(yīng)直線
library(effects) plot(effect("hp:wt", fitcar, xlevels=list(wt=c(2.2, 3.2, 4.2))), multiline=TRUE)
-
-
回歸診斷
-
標(biāo)準(zhǔn)方法:最常見(jiàn)的方法就是對(duì)lm()函數(shù)返回的對(duì)象使用plot()函數(shù)引瀑,可以生成評(píng)價(jià)模型擬合情況的四幅圖形
image正態(tài)性:當(dāng)預(yù)測(cè)變量值固定時(shí)灭衷,因變量成正態(tài)分布瞳遍,則殘差值也應(yīng)該是一個(gè)均值為0的正態(tài)分布舔箭。正態(tài)Q-Q圖(Normal Q-Q奢讨,右上)是在正態(tài)分布對(duì)應(yīng)的值下描沟,標(biāo)準(zhǔn)化殘差的概率圖佩伤。若滿足正態(tài)假設(shè)邀层,那么圖上的點(diǎn)應(yīng)該落在呈45度角的直線上;若不是如此,那么就違反了正態(tài)性的假設(shè)。
獨(dú)立性:從這些圖中分辨出因變量值是否相互獨(dú)立放刨,只能從收集的數(shù)據(jù)中來(lái)驗(yàn)證答毫。
線性:若因變量與自變量線性相關(guān),那么殘差值與預(yù)測(cè)(擬合)值就沒(méi)有任何系統(tǒng)關(guān)聯(lián)法严。換句話說(shuō),除了白噪聲迁客,模型應(yīng)該包含數(shù)據(jù)中所有的系統(tǒng)方差。在“殘差圖與擬合圖”(Residuals vs Fitted棉浸,左上)中可以清楚的看到一個(gè)曲線關(guān)系知押,這暗示著你可能需要對(duì)回歸模型加上一個(gè)二次項(xiàng)泌射。
同方差性:若滿足不變方差假設(shè),那么在位置尺度圖(Scale-Location Graph鬓照,左下)中熔酷,水平線周圍的點(diǎn)應(yīng)該隨機(jī)分布。該圖似乎滿足此假設(shè)豺裆。
-
最后一幅“殘差與杠桿圖”(Residuals vs Leverage拒秘,右下)提供了你可能關(guān)注的單個(gè)觀測(cè)點(diǎn)的信息。從圖形可以鑒別出離群點(diǎn)臭猜、高杠桿值點(diǎn)和強(qiáng)影響點(diǎn)躺酒。
- 一個(gè)觀測(cè)點(diǎn)是離群點(diǎn),表明擬合回歸模型對(duì)其預(yù)測(cè)效果不佳(產(chǎn)生了巨大的或正或負(fù)的殘差)蔑歌。
- 一個(gè)觀測(cè)點(diǎn)有很高的杠桿值羹应,表明它是一個(gè)異常的預(yù)測(cè)變量值的組合。也就是說(shuō)次屠,在預(yù)測(cè)變量空間中园匹,它是一個(gè)離群點(diǎn)。因變量值不參與計(jì)算一個(gè)觀測(cè)點(diǎn)的杠桿值劫灶。
- 一個(gè)觀測(cè)點(diǎn)是強(qiáng)影響點(diǎn)(influential observation)裸违,表明它對(duì)模型參數(shù)的估計(jì)產(chǎn)生的影響過(guò)大,非常不成比例本昏。強(qiáng)影響點(diǎn)可以通過(guò)Cook距離即Cook’s D統(tǒng)計(jì)量來(lái)鑒別供汛。
-
改進(jìn)的方法
-
car包提供了大量函數(shù),大大增強(qiáng)了擬合和評(píng)價(jià)回歸模型的能力
函 數(shù) 目 的 qqPlot() 分位數(shù)比較圖 durbinWatsonTest() 對(duì)誤差自相關(guān)性做Durbin-Watson檢驗(yàn) crPlots() 成分與殘差圖 ncvTest() 對(duì)非恒定的誤差方差做得分檢驗(yàn) spreadLevelPlot() 分散水平檢驗(yàn) outlierTest() Bonferroni離群點(diǎn)檢驗(yàn) avPlots() 添加的變量圖形 inluencePlot() 回歸影響圖 scatterplot() 增強(qiáng)的散點(diǎn)圖 scatterplotMatrix() 增強(qiáng)的散點(diǎn)圖矩陣 vif() 方差膨脹因子 - 正態(tài)性:與基礎(chǔ)包中的plot()函數(shù)相比涌穆,qqPlot()函數(shù)提供了更為精確的正態(tài)假設(shè)檢驗(yàn)方法紊馏,它畫(huà)出了在n-p-1個(gè)自由度的t分布下的學(xué)生化殘差(studentized residual,也稱學(xué)生化刪除殘差或折疊化殘差)圖形蒲犬,其中n是樣本大小朱监,p是回歸參數(shù)的數(shù)目(包括截距項(xiàng))。
- 獨(dú)立性:Durbin-Watson檢驗(yàn)函數(shù)原叮,能夠檢測(cè)誤差的序列相關(guān)性赫编,該檢驗(yàn)適用于時(shí)間獨(dú)立的數(shù)據(jù),對(duì)于非聚集型的數(shù)據(jù)并不適用
- 線性:通過(guò)成分殘差圖(component plus residual plot)也稱偏殘差圖(partial residual plot)奋隶,可以看看因變量與自變量之間是否呈非線性關(guān)系擂送,也可以看看是否有不同于已設(shè)定線性模型的系統(tǒng)偏差,圖形可用car包中的crPlots()函數(shù)繪制唯欣。若圖形存在非線性嘹吨,則說(shuō)明你可能對(duì)預(yù)測(cè)變量的函數(shù)形式建模不夠充分捡絮,那么就需要添加一些曲線成分尸疆,比如多項(xiàng)式項(xiàng)枉氮,或?qū)σ粋€(gè)或多個(gè)變量進(jìn)行變換(如用log(X)代替X)允粤,或用其他回歸變體形式而不是線性回歸
- 同方差性:car包提供了兩個(gè)有用的函數(shù),可以判斷誤差方差是否恒定问芬。ncvTest()函數(shù)生成一個(gè)計(jì)分檢驗(yàn)悦析,零假設(shè)為誤差方差不變,備擇假設(shè)為誤差方差隨著擬合值水平的變化而變化此衅。若檢驗(yàn)顯著强戴,則說(shuō)明存在異方差性(誤差方差不恒定)。spreadLevelPlot()函數(shù)創(chuàng)建一個(gè)添加了最佳擬合曲線的散點(diǎn)圖挡鞍,展示標(biāo)準(zhǔn)化殘差絕對(duì)值與擬合值的關(guān)系骑歹;通過(guò)分布水平圖看到這一點(diǎn),其中的點(diǎn)在水平的最佳擬合曲線周圍呈水平隨機(jī)分布墨微。若違反了該假設(shè)你將會(huì)看到一個(gè)非水平的曲線道媚。代碼結(jié)果建議冪次變換(suggested power transformation)的含義是,經(jīng)過(guò)p次冪(Y p)變換欢嘿,非恒定的誤差方差將會(huì)平穩(wěn)衰琐。例如也糊,若圖形顯示出了非水平趨勢(shì)炼蹦,建議冪次轉(zhuǎn)換為0.5,在回歸等式中用 Y 代替Y狸剃,可能會(huì)使模型滿足同方差性掐隐。若建議冪次為0,則使用對(duì)數(shù)變換钞馁。
-
-
線性模型假設(shè)的綜合驗(yàn)證
library(gvlma) gvmodel <- gvlma(fit2) summary(gvmodel)
多重共線性
多重共線性可用統(tǒng)計(jì)量VIF(Variance Inflation Factor虑省,方差膨脹因子)進(jìn)行檢測(cè)。VIF的平方根表示變量回歸參數(shù)的置信區(qū)間能膨脹為與模型無(wú)關(guān)的預(yù)測(cè)變量的程度(因此而得名)僧凰。car包中的vif()函數(shù)提供VIF值探颈。一般原則下, vif的平方根>2就表明存在多重共線性問(wèn)題
-
-
異常觀測(cè)值
一個(gè)全面的回歸分析要覆蓋對(duì)異常值的分析训措,包括離群點(diǎn)伪节、高杠桿值點(diǎn)和強(qiáng)影響點(diǎn)
-
離群點(diǎn)
- 離群點(diǎn)是指那些模型預(yù)測(cè)效果不佳的觀測(cè)點(diǎn)。它們通常有很大的绩鸣、或正或負(fù)的殘差怀大。正的殘差說(shuō)明模型低估了響應(yīng)值,負(fù)的殘差則說(shuō)明高估了響應(yīng)值呀闻。
- 鑒別離群點(diǎn)的方法:的Q-Q圖化借,落在置信區(qū)間帶外的點(diǎn)即可被認(rèn)為是離群點(diǎn)。另外一個(gè)粗糙的判斷準(zhǔn)則:標(biāo)準(zhǔn)化殘差值大于2或者小于?2的點(diǎn)可能是離群點(diǎn)捡多,需要特別關(guān)注蓖康。
- car包提供了一種離群點(diǎn)的統(tǒng)計(jì)檢驗(yàn)方法铐炫。outlierTest()函數(shù)可以求得最大標(biāo)準(zhǔn)化殘差絕對(duì)值Bonferroni調(diào)整后的p值.該函數(shù)只是根據(jù)單個(gè)最大(或正或負(fù))殘差值的顯著性來(lái)判斷是否有離群點(diǎn)。若不顯著钓瞭,則說(shuō)明數(shù)據(jù)集中沒(méi)有離群點(diǎn)驳遵;若顯著,則你必須刪除該離群點(diǎn)山涡,然后再檢驗(yàn)是否還有其他離群點(diǎn)存在
-
高杠桿值點(diǎn)
- 高杠桿值觀測(cè)點(diǎn)堤结,即是與其他預(yù)測(cè)變量有關(guān)的離群點(diǎn)。換句話說(shuō)鸭丛,它們是由許多異常的預(yù)測(cè)變量值組合起來(lái)的竞穷,與響應(yīng)變量值沒(méi)有關(guān)系。
- 高杠桿值的觀測(cè)點(diǎn)可通過(guò)帽子統(tǒng)計(jì)量(hat statistic)判斷鳞溉。對(duì)于一個(gè)給定的數(shù)據(jù)集瘾带,帽子均值為p/n,其中p 是模型估計(jì)的參數(shù)數(shù)目(包含截距項(xiàng))熟菲,n是樣本量看政。一般來(lái)說(shuō),若觀測(cè)點(diǎn)的帽子值大于帽子均值的2或3倍抄罕,即可以認(rèn)定為高杠桿值點(diǎn)允蚣。
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))) } hat.plot(fit)
-
強(qiáng)影響點(diǎn)
強(qiáng)影響點(diǎn),即對(duì)模型參數(shù)估計(jì)值影響有些比例失衡的點(diǎn)呆贿。
有兩種方法可以檢測(cè)強(qiáng)影響點(diǎn):Cook距離嚷兔,或稱D統(tǒng)計(jì)量,以及變量添加圖(added variableplot)做入。
-
Cook’s D值大于4/(n-k-1)冒晰,則表明它是強(qiáng)影響點(diǎn),其中n 為樣本量大小竟块,k 是預(yù)測(cè)變量數(shù)目
# Cooks Distance D # identify D values > 4/(n-k-1) cutoff <- 4/(nrow(states)-length(fit$coefficients)-2) plot(fit, which=4, cook.levels=cutoff) abline(h=cutoff, lty=2, col="red")
-
Cook’s D圖有助于鑒別強(qiáng)影響點(diǎn)壶运,但是并不提供關(guān)于這些點(diǎn)如何影響模型的信息。變量添加圖彌補(bǔ)了這個(gè)缺陷浪秘。對(duì)于一個(gè)響應(yīng)變量和k個(gè)預(yù)測(cè)變量蒋情,你可以創(chuàng)建k個(gè)變量添加圖。所謂變量添加圖秫逝,即對(duì)于每個(gè)預(yù)測(cè)變量Xk恕出,繪制Xk 在其他k-1個(gè)預(yù)測(cè)變量上回歸的殘差值相對(duì)于響應(yīng)變量在其他k-1個(gè)預(yù)測(cè)變量上回歸的殘差值的關(guān)系圖。car包中的avPlots()函數(shù)可提供變量添加圖
avPlots(fit, ask=FALSE, id.method="identify")
-
利用car包中的influencePlot()函數(shù)违帆,你還可以將離群點(diǎn)浙巫、杠桿值和強(qiáng)影響點(diǎn)的信息整合到一幅圖形中
influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportial to Cook's Distance" )
-
-
改進(jìn)措施
有四種方法可以處理違背回歸假設(shè)的問(wèn)題:
-
刪除觀測(cè)點(diǎn):
刪除離群點(diǎn)通常可以提高數(shù)據(jù)集對(duì)于正態(tài)假設(shè)的擬合度,而強(qiáng)影響點(diǎn)會(huì)干擾結(jié)果的畴,通常也會(huì)被刪除渊抄。刪除最大的離群點(diǎn)或者強(qiáng)影響點(diǎn)后,模型需要重新擬合丧裁。若離群點(diǎn)或強(qiáng)影響點(diǎn)仍然存在护桦,重復(fù)以上過(guò)程直至獲得比較滿意的擬合。
-
變量變換煎娇;
當(dāng)模型不符合正態(tài)性二庵、線性或者同方差性假設(shè)時(shí),一個(gè)或多個(gè)變量的變換通郴呵海可以改善或調(diào)整模型效果
-
常見(jiàn)的變換
imagesummary(powerTransform(states$Murder)) bcPower Transformation to Normality Est Power Rounded Pwr Wald Lwr bnd Wald Upr Bnd states$Murder 0.6055 1 0.0884 1.1227 Likelihood ratio tests about transformation parameters LRT df pval LR test, lambda = (0) 5.665991 1 0.01729694 LR test, lambda = (1) 2.122763 1 0.14512456 ## 結(jié)果表明催享,你可以用Murder0.6來(lái)正態(tài)化變量Murder。 ## 由于0.6很接近0.5哟绊,你可以嘗試用平方根變換來(lái)提高模型正態(tài)性的符合程度因妙。 ## 但在本例中,λ= 1的假設(shè)也無(wú)法拒絕(p=0.145)票髓,因此沒(méi)有強(qiáng)有力的證據(jù)表明本例需要變量變換
-
當(dāng)違反了線性假設(shè)時(shí)攀涵,對(duì)預(yù)測(cè)變量進(jìn)行變換常常會(huì)比較有用。car包中的boxTidwell()函數(shù)通過(guò)獲得預(yù)測(cè)變量?jī)鐢?shù)的最大似然估計(jì)來(lái)改善線性關(guān)系
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 ## 結(jié)果顯示洽沟,使用變換Population0.87和Illiteracy1.36能夠大大改善線性關(guān)系 ## 但是對(duì)Population(p=0.75)和Illiteracy(p=0.54)的計(jì)分檢驗(yàn)又表明變量并不需要變換以故。
-
添加或刪除變量
- 最常見(jiàn)的方法就是刪除某個(gè)存在多重共線性的變量(某個(gè)變量 sqrt(vif)> 2 )。
- 另外一個(gè)可用的方法便是嶺回歸——多元回歸的變體玲躯,專門(mén)用來(lái)處理多重共線性問(wèn)題
使用其他回歸方法据德。
-
-
選擇“最佳”的回歸模型
- 模型比較
- 用基礎(chǔ)安裝中的anova()函數(shù)可以比較兩個(gè)嵌套模型的擬合優(yōu)度鳄乏。所謂嵌套模型跷车,即它的一些項(xiàng)完全包含在另一個(gè)模型中
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) fit2 <- lm(Murder ~ Population + Illiteracy, data=states) anova(fit2, fit1)
- AIC(Akaike Information Criterion,赤池信息準(zhǔn)則)也可以用來(lái)比較模型橱野,它考慮了模型的統(tǒng)計(jì)擬合度以及用來(lái)擬合的參數(shù)數(shù)目朽缴。AIC值越小的模型要優(yōu)先選擇,它說(shuō)明模型用較少的參數(shù)獲得了足夠的擬合度
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) fit2 <- lm(Murder ~ Population + Illiteracy, data=states) AIC(fit1,fit2)
- 用基礎(chǔ)安裝中的anova()函數(shù)可以比較兩個(gè)嵌套模型的擬合優(yōu)度鳄乏。所謂嵌套模型跷车,即它的一些項(xiàng)完全包含在另一個(gè)模型中
- 變量選擇
-
逐步回歸法(stepwise method):逐步回歸中水援,模型會(huì)一次添加或者刪除一個(gè)變量密强,直到達(dá)到某個(gè)判停準(zhǔn)則為止。
- 向前逐步回歸(forward stepwise)每次添加一個(gè)預(yù)測(cè)變量到模型中蜗元,直到添加變量不會(huì)使模型有所改進(jìn)為止或渤。
- 向后逐步回歸(backward stepwise)從模型包含所有預(yù)測(cè)變量開(kāi)始,一次刪除一個(gè)變量直到會(huì)降低模型質(zhì)量為止奕扣。
- 向前向后逐步回歸(stepwise stepwise薪鹦,通常稱作逐步回歸),結(jié)合了向前逐步回歸和向后逐步回歸的方法,變量每次進(jìn)入一個(gè)池磁,但是每一步中奔害,變量都會(huì)被重新評(píng)價(jià),對(duì)模型沒(méi)有貢獻(xiàn)的變量將會(huì)被刪除地熄,預(yù)測(cè)變量可能會(huì)被添加华临、刪除好幾次,直到獲得最優(yōu)模型為止端考。
- MASS包中的stepAIC()函數(shù)可以實(shí)現(xiàn)逐步回歸模型(向前雅潭、向后和向前向后),依據(jù)的是精確AIC準(zhǔn)則
library(MASS) states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) stepAIC(fit, direction="backward") 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 <none> 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 <none> 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 <none> 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 ## 逐步回歸法其實(shí)存在爭(zhēng)議却特,雖然它可能會(huì)找到一個(gè)好的模型寻馏,但是不能保證模型就是最佳模型,因?yàn)椴皇敲恳粋€(gè)可能的模型都被評(píng)價(jià)了
-
全子集回歸(all-subsets regression)核偿,即所有可能的模型都會(huì)被檢驗(yàn)诚欠,可以選擇展示所有可能的結(jié)果,也可以展示n 個(gè)不同子集大醒馈(一個(gè)轰绵、兩個(gè)或多個(gè)預(yù)測(cè)變量)的最佳模型
- 全子集回歸可用leaps包中的regsubsets()函數(shù)實(shí)現(xiàn)。你能通過(guò)R平方尼荆、調(diào)整R平方或Mallows Cp統(tǒng)計(jì)量等準(zhǔn)則來(lái)選擇“最佳”模型
- R平方含義是預(yù)測(cè)變量解釋響應(yīng)變量的程度
- Mallows Cp統(tǒng)計(jì)量也用來(lái)作為逐步回歸的判停規(guī)則左腔。廣泛研究表明,對(duì)于一個(gè)好的模型捅儒,它的Cp統(tǒng)計(jì)量非常接近于模型的參數(shù)數(shù)目(包括截距項(xiàng))
- 可用leaps包中的plot()函數(shù)繪制液样,或者用car包中的subsets()函數(shù)繪制全子集回歸圖
library(leaps) states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) leaps <-regsubsets(Murder ~ Population + Illiteracy + Income + Frost, data=states, nbest=4) plot(leaps, scale="adjr2") library(car) subsets(leaps, statistic="cp", main="Cp Plot for All Subsets Regression") abline(1,1,lty=2,col="red")
-
- 模型比較
-
深層次分析
-
交叉驗(yàn)證,評(píng)價(jià)回歸方程的泛化能力
- 所謂交叉驗(yàn)證巧还,即將一定比例的數(shù)據(jù)挑選出來(lái)作為訓(xùn)練樣本鞭莽,另外的樣本作保留樣本,先在訓(xùn)練樣本上獲取回歸方程麸祷,然后在保留樣本上做預(yù)測(cè)澎怒。由于保留樣本不涉及模型參數(shù)的選擇,該樣本可獲得比新數(shù)據(jù)更為精確的估計(jì)阶牍。
- k 重交叉驗(yàn)證喷面,在k 重交叉驗(yàn)證中,樣本被分為k個(gè)子樣本走孽,輪流將k?1個(gè)子樣本組合作為訓(xùn)練集惧辈,另外1個(gè)子樣本作為保留集。這樣會(huì)獲得k 個(gè)預(yù)測(cè)方程磕瓷,記錄k個(gè)保留樣本的預(yù)測(cè)表現(xiàn)結(jié)果盒齿,然后求其平均值。
- bootstrap 包中的 crossval() 函數(shù)可以實(shí)現(xiàn) k 重交叉驗(yàn)證。
- 通過(guò)選擇有更好泛化能力的模型县昂,還可以用交叉驗(yàn)證來(lái)挑選變量肮柜。
- 下方示例,通過(guò)定義shrinkage()函數(shù)對(duì)模型的R平方統(tǒng)計(jì)量做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 # raw R2 r2cv <- cor(y,results$cv.fit)**2 # cross-validated R2 cat("Original R-square =", r2, "\n") cat(k, "Fold Cross-Validated R-square =", r2cv, "\n") cat("Change =", r2-r2cv, "\n") } # using it states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit <- lm(Murder ~ Population + Income + Illiteracy + Frost, data=states) shrinkage(fit) fit2 <- lm(Murder~Population+Illiteracy,data=states) shrinkage(fit2)
-
相對(duì)重要性
- 最簡(jiǎn)單的莫過(guò)于比較標(biāo)準(zhǔn)化的回歸系數(shù)倒彰,它表示當(dāng)其他預(yù)測(cè)變量不變時(shí)审洞,該預(yù)測(cè)變量一個(gè)標(biāo)準(zhǔn)差的變化可引起的響應(yīng)變量的預(yù)期變化(以標(biāo)準(zhǔn)差單位度量)。在進(jìn)行回歸分析前待讳,可用scale()函數(shù)將數(shù)據(jù)標(biāo)準(zhǔn)化為均值為0芒澜、標(biāo)準(zhǔn)差為1的數(shù)據(jù)集,這樣用R回歸即可獲得標(biāo)準(zhǔn)化的回歸系數(shù)创淡。(注意痴晦,scale()函數(shù)返回的是一個(gè)矩陣,而lm()函數(shù)要求一個(gè)數(shù)據(jù)框琳彩,你需要用一個(gè)中間步驟來(lái)轉(zhuǎn)換一下誊酌。)
states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) zstates <- as.data.frame(scale(states)) zfit <- lm(Murder~Population + Income + Illiteracy + Frost, data=zstates) coef(zfit) (Intercept) Population Income Illiteracy Frost -2.054026e-16 2.705095e-01 1.072372e-02 6.840496e-01 8.185407e-03 ## 當(dāng)其他因素不變時(shí),文盲率一個(gè)標(biāo)準(zhǔn)差的變化將增加0.68個(gè)標(biāo)準(zhǔn)差的謀殺率露乏。根據(jù)標(biāo)準(zhǔn)化的回歸系數(shù)碧浊,我們可認(rèn)為Illiteracy是最重要的預(yù)測(cè)變量,而Frost是最不重要的
- 相對(duì)權(quán)重(relative weight)是一種比較有前景的新方法瘟仿,它是對(duì)所有可能子模型添加一個(gè)預(yù)測(cè)變量引起的R平方平均增加量的一個(gè)近似值
- 下方示例提供了一個(gè)生成相對(duì)權(quán)重的函數(shù)
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)) lambda <- evec %*% delta %*% t(evec) lambdasq <- lambda ^ 2 beta <- solve(lambda) %*% rxy rsquare <- colSums(beta ^ 2) rawwgt <- lambdasq %*% beta ^ 2 import <- (rawwgt / rsquare) * 100 import <- as.data.frame(import) row.names(import) <- names(fit$model[2:nvar]) names(import) <- "Weights" import <- import[order(import),1, drop=FALSE] dotchart(import$Weights, labels=row.names(import), xlab="% of R-Square", pch=19, main="Relative Importance of Predictor Variables", sub=paste("Total R-Square=", round(rsquare, digits=3)), ...) return(import) } # Applying the relweights function states <- as.data.frame(state.x77[,c("Murder", "Population", "Illiteracy", "Income", "Frost")]) fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) relweights(fit, col="blue") ## 對(duì)結(jié)果進(jìn)行查看箱锐,可以看到各個(gè)預(yù)測(cè)變量對(duì)模型方差的解釋程度(R平方=0.567),Illiteracy解釋了59%的R平方劳较,F(xiàn)rost解釋了20.79%驹止。根據(jù)相對(duì)權(quán)重法,Illiteracy有最大的相對(duì)重要性观蜗,余下依次是Frost臊恋、Population和Income
-