20180404（從有道遷移）
回歸
-
回歸的多面性
回歸分析的各種變體
回歸類型 用 途 簡單線性 用一個量化的解釋變量預(yù)測一個量化的響應(yīng)變量 多項式 用一個量化的解釋變量預(yù)測一個量化的響應(yīng)變量，模型的關(guān)系是n階多項式 多元線性 用兩個或多個量化的解釋變量預(yù)測一個量化的響應(yīng)變量 多變量 用一個或多個解釋變量預(yù)測多個響應(yīng)變量 Logistic 用一個或多個解釋變量預(yù)測一個類別型響應(yīng)變量 泊松 用一個或多個解釋變量預(yù)測一個代表頻數(shù)的響應(yīng)變量 Cox比例風(fēng)險 用一個或多個解釋變量預(yù)測一個事件（死亡、失敗或舊病復(fù)發(fā)）發(fā)生的時間 時間序列 對誤差項相關(guān)的時間序列數(shù)據(jù)建模 非線性 用一個或多個量化的解釋變量預(yù)測一個量化的響應(yīng)變量，不過模型是非線性的 非參數(shù) 用一個或多個量化的解釋變量預(yù)測一個量化的響應(yīng)變量，模型的形式源自數(shù)據(jù)形式，不事先設(shè)定 穩(wěn)健 用一個或多個量化的解釋變量預(yù)測一個量化的響應(yīng)變量，能抵御強影響點的干擾 -
OLS 回歸
-
普通最小二乘（OLS）回歸法的適用情境
OLS回歸是通過預(yù)測變量的加權(quán)和來預(yù)測量化的因變量，其中權(quán)重是通過數(shù)據(jù)估計而得的參數(shù)
-
OLS回歸擬合模型的形式,其中，n 為觀測的數(shù)目，k 為預(yù)測變量的數(shù)目。
\widehat{y}_{i}=\widehat{\beta}_{0}+\widehat{\beta}_{1}X_{1i}+...+\widehat{\beta}_{k}X_{ki} i=1...nimage 為了能夠恰當(dāng)?shù)亟忉孫LS模型的系數(shù)，數(shù)據(jù)必須滿足以下統(tǒng)計假設(shè)。
正態(tài)性 對于固定的自變量值，因變量值成正態(tài)分布。
獨立性 Yi值之間相互獨立。
線性 因變量與自變量之間為線性相關(guān)。
-
同方差性 因變量的方差不隨自變量的水平不同而變化。也可稱作不變方差，但是說同方差性感覺上更犀利。
如果違背了以上假設(shè)，你的統(tǒng)計顯著性檢驗結(jié)果和所得的置信區(qū)間很可能就不精確。注意，OLS回歸還假定自變量是固定的且測量無誤差，但在實踐中通常都放松了這個假設(shè)
-
用lm()擬合回歸模型
-
在R中，擬合線性模型最基本的函數(shù)就是lm()，格式為：
myfit <- lm(formula,data),其中:- formula指要擬合的模型形式
- data是一個數(shù)據(jù)框，包含了用于擬合模型的數(shù)據(jù)。
- 結(jié)果對象（本例中是myfit）存儲在一個列表中，包含了所擬合模型的大量信息。
-
表達(dá)式（formula）形式如下：
y ~ x1+x2+x3+...+xk- ~左邊為響應(yīng)變量
- 右邊為各個預(yù)測變量，預(yù)測變量之間用+符號分隔
-
R表達(dá)式中常用的符號
符號 用 途 ~ 分隔符號，左邊為響應(yīng)變量，右邊為解釋變量。</br>例如，要通過x、z和w預(yù)測y，代碼為y ~ x + z + w + 分隔預(yù)測變量 : 表示預(yù)測變量的交互項。</br>例如，要通過x、z及x與z的交互項預(yù)測y，代碼為y ~ x + z + x:z * 表示所有可能交互項的簡潔方式。</br>代碼y~ x * z * w可展開為y ~ x + z + w + x:z + x:w + z:w +x:z:w ^ 表示交互項達(dá)到某個次數(shù)。</br>代碼y ~ (x + z + w)^2可展開為y ~ x + z + w + x:z + x:w + z:w . 表示包含除因變量外的所有變量。</br>例如，若一個數(shù)據(jù)框包含變量x、y、z和w，代碼y ~ .可展開為y ~ x +z + w - 減號，表示從等式中移除某個變量。</br>例如，y ~ (x + z + w)^2 – x:w可展開為y ~ x + z + w + x:z + z:w -1 刪除截距項。</br>例如，表達(dá)式y(tǒng) ~ x - 1擬合y在x上的回歸，并強制直線通過原點 I() 從算術(shù)的角度來解釋括號中的元素。</br>例如，y ~ x + (z + w)^2將展開為y ~ x + z + w + z:w。</br>相反, 代碼y~ x + I((z + w)^2)將展開為y ~ x + h, h是一個由z和w的平方和創(chuàng)建的新變量 function 可以在表達(dá)式中用的數(shù)學(xué)函數(shù)。</br>例如，log(y) ~ x + z + w表示通過x、z和w來預(yù)測log(y) -
對擬合線性模型非常有用的其他函數(shù)
函 數(shù) 用 途 summary() 展示擬合模型的詳細(xì)結(jié)果 coefficients() 列出擬合模型的模型參數(shù)（截距項和斜率） confint() 提供模型參數(shù)的置信區(qū)間（默認(rèn)95%） fitted() 列出擬合模型的預(yù)測值 residuals() 列出擬合模型的殘差值 anova() 生成一個擬合模型的方差分析表，或者比較兩個或更多擬合模型的方差分析表 vcov() 列出模型參數(shù)的協(xié)方差矩陣 AIC() 輸出赤池信息統(tǒng)計量 plot() 生成評價擬合模型的診斷圖 predict() 用擬合模型對新的數(shù)據(jù)集預(yù)測響應(yīng)變量值
-
-
簡單線性回歸：當(dāng)回歸模型包含一個因變量和一個自變量時，我們稱為簡單線性回歸
- 示例：基礎(chǔ)安裝中的數(shù)據(jù)集women提供了15個年齡在30~39歲間女性的身高和體重信息，我們想通過身高來預(yù)測體重，獲得一個等式可以幫助我們分辨出那些過重或過瘦的個體
- 代碼：
## 擬合回歸模型 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平方項（0.991）表明模型可以解釋體重99.1%的方差，它也是實際和預(yù)測值之間的相關(guān)系數(shù)（R2 = r2?Y）。 ## 殘差標(biāo)準(zhǔn)誤（1.53 lbs）則可認(rèn)為是模型用身高預(yù)測體重的平均誤差。 ## F統(tǒng)計量檢驗所有的預(yù)測變量預(yù)測響應(yīng)變量是否都在某個幾率水平之上。由于簡單回歸只有一個預(yù)測變量，此處F檢驗等同于身高回歸系數(shù)的t檢驗 women$weight ## 列出擬合模型的預(yù)測值 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)
-
多項式回歸,。當(dāng)只有一個預(yù)測變量，但同時包含變量的冪（比如，X、X2、X3）時，我們稱之為多項式回歸
針對上訴身高問題，你可以通過添加一個二次項（即
X2）來提高回歸的預(yù)測精度，因此可以擬合含二次項的等式：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%。 ## 二次項的顯著性（t = 13.89，p<0.001）表明包含二次項提高了模型的擬合度。 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)系圖,==這個功能加強的圖形，既提供了身高與體重的散點圖、線性擬合曲線和平滑擬合（loess）曲線，還在相應(yīng)邊界展示了每個變量的箱線圖==
- spread=FALSE選項刪除了殘差正負(fù)均方根在平滑曲線上的展開和非對稱信息。
- lty.smooth=2選項設(shè)置loess擬合曲線為虛線。
- pch=19選項設(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)有不止一個預(yù)測變量時，則稱為多元線性回歸
多元回歸分析中，第一步最好檢查一下變量間的相關(guān)性。cor()函數(shù)提供了二變量之間的相關(guān)系數(shù)，car包中scatterplotMatrix()函數(shù)則會生成散點圖矩陣
-
使用lm()函數(shù)擬合多元線性回歸模型
- 當(dāng)預(yù)測變量不止一個時，回歸系數(shù)的含義為，一個預(yù)測變量增加一個單位，其他預(yù)測變量保持不變時，因變量將要增加的數(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)
-
有交互項的多元線性回歸
在初步判斷多個變量之間存在互相作用時，可以將變量之間的交互項加入擬合模型中，通過對交互項在模型中的pvalue驗證是否相關(guān)性顯著
-
可以通過effects包中的effect()函數(shù)，你可以用圖形展示交互項的結(jié)果
- 格式：plot(effect(term,mod,xlevels),multiple=TRUE)
- term即模型要畫的項，mod為通過lm()擬合的模型，xlevels是一個列表，指定變量要設(shè)定的常量值，multiline=TRUE選項表示添加相應(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)方法：最常見的方法就是對lm()函數(shù)返回的對象使用plot()函數(shù)，可以生成評價模型擬合情況的四幅圖形
image正態(tài)性：當(dāng)預(yù)測變量值固定時，因變量成正態(tài)分布，則殘差值也應(yīng)該是一個均值為0的正態(tài)分布。正態(tài)Q-Q圖（Normal Q-Q，右上）是在正態(tài)分布對應(yīng)的值下，標(biāo)準(zhǔn)化殘差的概率圖。若滿足正態(tài)假設(shè)，那么圖上的點應(yīng)該落在呈45度角的直線上；若不是如此，那么就違反了正態(tài)性的假設(shè)。
獨立性：從這些圖中分辨出因變量值是否相互獨立，只能從收集的數(shù)據(jù)中來驗證。
線性：若因變量與自變量線性相關(guān)，那么殘差值與預(yù)測（擬合）值就沒有任何系統(tǒng)關(guān)聯(lián)。換句話說，除了白噪聲，模型應(yīng)該包含數(shù)據(jù)中所有的系統(tǒng)方差。在“殘差圖與擬合圖”（Residuals vs Fitted，左上）中可以清楚的看到一個曲線關(guān)系，這暗示著你可能需要對回歸模型加上一個二次項。
同方差性：若滿足不變方差假設(shè)，那么在位置尺度圖（Scale-Location Graph，左下）中，水平線周圍的點應(yīng)該隨機分布。該圖似乎滿足此假設(shè)。
-
最后一幅“殘差與杠桿圖”（Residuals vs Leverage，右下）提供了你可能關(guān)注的單個觀測點的信息。從圖形可以鑒別出離群點、高杠桿值點和強影響點。
- 一個觀測點是離群點，表明擬合回歸模型對其預(yù)測效果不佳（產(chǎn)生了巨大的或正或負(fù)的殘差）。
- 一個觀測點有很高的杠桿值，表明它是一個異常的預(yù)測變量值的組合。也就是說，在預(yù)測變量空間中，它是一個離群點。因變量值不參與計算一個觀測點的杠桿值。
- 一個觀測點是強影響點（influential observation），表明它對模型參數(shù)的估計產(chǎn)生的影響過大，非常不成比例。強影響點可以通過Cook距離即Cook’s D統(tǒng)計量來鑒別。
-
改進(jìn)的方法
-
car包提供了大量函數(shù)，大大增強了擬合和評價回歸模型的能力
函 數(shù) 目 的 qqPlot() 分位數(shù)比較圖 durbinWatsonTest() 對誤差自相關(guān)性做Durbin-Watson檢驗 crPlots() 成分與殘差圖 ncvTest() 對非恒定的誤差方差做得分檢驗 spreadLevelPlot() 分散水平檢驗 outlierTest() Bonferroni離群點檢驗 avPlots() 添加的變量圖形 inluencePlot() 回歸影響圖 scatterplot() 增強的散點圖 scatterplotMatrix() 增強的散點圖矩陣 vif() 方差膨脹因子 - 正態(tài)性：與基礎(chǔ)包中的plot()函數(shù)相比，qqPlot()函數(shù)提供了更為精確的正態(tài)假設(shè)檢驗方法，它畫出了在n-p-1個自由度的t分布下的學(xué)生化殘差（studentized residual，也稱學(xué)生化刪除殘差或折疊化殘差）圖形，其中n是樣本大小，p是回歸參數(shù)的數(shù)目（包括截距項）。
- 獨立性：Durbin-Watson檢驗函數(shù)，能夠檢測誤差的序列相關(guān)性，該檢驗適用于時間獨立的數(shù)據(jù)，對于非聚集型的數(shù)據(jù)并不適用
- 線性：通過成分殘差圖（component plus residual plot）也稱偏殘差圖（partial residual plot），可以看看因變量與自變量之間是否呈非線性關(guān)系，也可以看看是否有不同于已設(shè)定線性模型的系統(tǒng)偏差，圖形可用car包中的crPlots()函數(shù)繪制。若圖形存在非線性，則說明你可能對預(yù)測變量的函數(shù)形式建模不夠充分，那么就需要添加一些曲線成分，比如多項式項，或?qū)σ粋€或多個變量進(jìn)行變換（如用log(X)代替X），或用其他回歸變體形式而不是線性回歸
- 同方差性：car包提供了兩個有用的函數(shù)，可以判斷誤差方差是否恒定。ncvTest()函數(shù)生成一個計分檢驗，零假設(shè)為誤差方差不變，備擇假設(shè)為誤差方差隨著擬合值水平的變化而變化。若檢驗顯著，則說明存在異方差性（誤差方差不恒定）。spreadLevelPlot()函數(shù)創(chuàng)建一個添加了最佳擬合曲線的散點圖，展示標(biāo)準(zhǔn)化殘差絕對值與擬合值的關(guān)系；通過分布水平圖看到這一點，其中的點在水平的最佳擬合曲線周圍呈水平隨機分布。若違反了該假設(shè)你將會看到一個非水平的曲線。代碼結(jié)果建議冪次變換（suggested power transformation）的含義是，經(jīng)過p次冪（Y p）變換，非恒定的誤差方差將會平穩(wěn)。例如，若圖形顯示出了非水平趨勢，建議冪次轉(zhuǎn)換為0.5，在回歸等式中用 Y 代替Y，可能會使模型滿足同方差性。若建議冪次為0，則使用對數(shù)變換。
-
-
線性模型假設(shè)的綜合驗證
library(gvlma) gvmodel <- gvlma(fit2) summary(gvmodel) 多重共線性
多重共線性可用統(tǒng)計量VIF（Variance Inflation Factor，方差膨脹因子）進(jìn)行檢測。VIF的平方根表示變量回歸參數(shù)的置信區(qū)間能膨脹為與模型無關(guān)的預(yù)測變量的程度（因此而得名）。car包中的vif()函數(shù)提供VIF值。一般原則下， vif的平方根>2就表明存在多重共線性問題
-
-
異常觀測值
一個全面的回歸分析要覆蓋對異常值的分析，包括離群點、高杠桿值點和強影響點
-
離群點
- 離群點是指那些模型預(yù)測效果不佳的觀測點。它們通常有很大的、或正或負(fù)的殘差。正的殘差說明模型低估了響應(yīng)值，負(fù)的殘差則說明高估了響應(yīng)值。
- 鑒別離群點的方法：的Q-Q圖，落在置信區(qū)間帶外的點即可被認(rèn)為是離群點。另外一個粗糙的判斷準(zhǔn)則：標(biāo)準(zhǔn)化殘差值大于2或者小于?2的點可能是離群點，需要特別關(guān)注。
- car包提供了一種離群點的統(tǒng)計檢驗方法。outlierTest()函數(shù)可以求得最大標(biāo)準(zhǔn)化殘差絕對值Bonferroni調(diào)整后的p值.該函數(shù)只是根據(jù)單個最大（或正或負(fù)）殘差值的顯著性來判斷是否有離群點。若不顯著，則說明數(shù)據(jù)集中沒有離群點；若顯著，則你必須刪除該離群點，然后再檢驗是否還有其他離群點存在
-
高杠桿值點
- 高杠桿值觀測點，即是與其他預(yù)測變量有關(guān)的離群點。換句話說，它們是由許多異常的預(yù)測變量值組合起來的，與響應(yīng)變量值沒有關(guān)系。
- 高杠桿值的觀測點可通過帽子統(tǒng)計量（hat statistic）判斷。對于一個給定的數(shù)據(jù)集，帽子均值為p/n，其中p 是模型估計的參數(shù)數(shù)目（包含截距項），n是樣本量。一般來說，若觀測點的帽子值大于帽子均值的2或3倍，即可以認(rè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) -
強影響點
強影響點，即對模型參數(shù)估計值影響有些比例失衡的點。
有兩種方法可以檢測強影響點：Cook距離，或稱D統(tǒng)計量，以及變量添加圖（added variableplot）。
-
Cook’s D值大于4/(n-k-1)，則表明它是強影響點，其中n 為樣本量大小，k 是預(yù)測變量數(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圖有助于鑒別強影響點，但是并不提供關(guān)于這些點如何影響模型的信息。變量添加圖彌補了這個缺陷。對于一個響應(yīng)變量和k個預(yù)測變量，你可以創(chuàng)建k個變量添加圖。所謂變量添加圖，即對于每個預(yù)測變量Xk，繪制Xk 在其他k-1個預(yù)測變量上回歸的殘差值相對于響應(yīng)變量在其他k-1個預(yù)測變量上回歸的殘差值的關(guān)系圖。car包中的avPlots()函數(shù)可提供變量添加圖
avPlots(fit, ask=FALSE, id.method="identify")
-
利用car包中的influencePlot()函數(shù)，你還可以將離群點、杠桿值和強影響點的信息整合到一幅圖形中
influencePlot(fit, id.method="identify", main="Influence Plot", sub="Circle size is proportial to Cook's Distance" )
-
-
改進(jìn)措施
有四種方法可以處理違背回歸假設(shè)的問題：
-
刪除觀測點:
刪除離群點通?？梢蕴岣邤?shù)據(jù)集對于正態(tài)假設(shè)的擬合度，而強影響點會干擾結(jié)果，通常也會被刪除。刪除最大的離群點或者強影響點后，模型需要重新擬合。若離群點或強影響點仍然存在，重復(fù)以上過程直至獲得比較滿意的擬合。
-
變量變換；
當(dāng)模型不符合正態(tài)性、線性或者同方差性假設(shè)時，一個或多個變量的變換通?？梢愿纳苹蛘{(diào)整模型效果
-
常見的變換
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來正態(tài)化變量Murder。 ## 由于0.6很接近0.5，你可以嘗試用平方根變換來提高模型正態(tài)性的符合程度。 ## 但在本例中，λ= 1的假設(shè)也無法拒絕（p=0.145），因此沒有強有力的證據(jù)表明本例需要變量變換 -
當(dāng)違反了線性假設(shè)時，對預(yù)測變量進(jìn)行變換常常會比較有用。car包中的boxTidwell()函數(shù)通過獲得預(yù)測變量冪數(shù)的最大似然估計來改善線性關(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)系 ## 但是對Population（p=0.75）和Illiteracy（p=0.54）的計分檢驗又表明變量并不需要變換。
-
添加或刪除變量
- 最常見的方法就是刪除某個存在多重共線性的變量（某個變量 sqrt(vif)> 2 ）。
- 另外一個可用的方法便是嶺回歸——多元回歸的變體，專門用來處理多重共線性問題
使用其他回歸方法。
-
-
選擇“最佳”的回歸模型
- 模型比較
- 用基礎(chǔ)安裝中的anova()函數(shù)可以比較兩個嵌套模型的擬合優(yōu)度。所謂嵌套模型，即它的一些項完全包含在另一個模型中
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)則）也可以用來比較模型，它考慮了模型的統(tǒng)計擬合度以及用來擬合的參數(shù)數(shù)目。AIC值越小的模型要優(yōu)先選擇，它說明模型用較少的參數(shù)獲得了足夠的擬合度
fit1 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states) fit2 <- lm(Murder ~ Population + Illiteracy, data=states) AIC(fit1,fit2)
- 用基礎(chǔ)安裝中的anova()函數(shù)可以比較兩個嵌套模型的擬合優(yōu)度。所謂嵌套模型，即它的一些項完全包含在另一個模型中
- 變量選擇
-
逐步回歸法（stepwise method）：逐步回歸中，模型會一次添加或者刪除一個變量，直到達(dá)到某個判停準(zhǔn)則為止。
- 向前逐步回歸（forward stepwise）每次添加一個預(yù)測變量到模型中，直到添加變量不會使模型有所改進(jìn)為止。
- 向后逐步回歸（backward stepwise）從模型包含所有預(yù)測變量開始，一次刪除一個變量直到會降低模型質(zhì)量為止。
- 向前向后逐步回歸（stepwise stepwise，通常稱作逐步回歸），結(jié)合了向前逐步回歸和向后逐步回歸的方法，變量每次進(jìn)入一個，但是每一步中，變量都會被重新評價，對模型沒有貢獻(xiàn)的變量將會被刪除，預(yù)測變量可能會被添加、刪除好幾次，直到獲得最優(yōu)模型為止。
- MASS包中的stepAIC()函數(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 ## 逐步回歸法其實存在爭議，雖然它可能會找到一個好的模型，但是不能保證模型就是最佳模型，因為不是每一個可能的模型都被評價了 -
全子集回歸（all-subsets regression），即所有可能的模型都會被檢驗，可以選擇展示所有可能的結(jié)果，也可以展示n 個不同子集大?。ㄒ粋€、兩個或多個預(yù)測變量）的最佳模型
- 全子集回歸可用leaps包中的regsubsets()函數(shù)實現(xiàn)。你能通過R平方、調(diào)整R平方或Mallows Cp統(tǒng)計量等準(zhǔn)則來選擇“最佳”模型
- R平方含義是預(yù)測變量解釋響應(yīng)變量的程度
- Mallows Cp統(tǒng)計量也用來作為逐步回歸的判停規(guī)則。廣泛研究表明，對于一個好的模型，它的Cp統(tǒng)計量非常接近于模型的參數(shù)數(shù)目（包括截距項）
- 可用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")
-
- 模型比較
-
深層次分析
-
交叉驗證，評價回歸方程的泛化能力
- 所謂交叉驗證，即將一定比例的數(shù)據(jù)挑選出來作為訓(xùn)練樣本，另外的樣本作保留樣本，先在訓(xùn)練樣本上獲取回歸方程，然后在保留樣本上做預(yù)測。由于保留樣本不涉及模型參數(shù)的選擇，該樣本可獲得比新數(shù)據(jù)更為精確的估計。
- k 重交叉驗證，在k 重交叉驗證中，樣本被分為k個子樣本，輪流將k?1個子樣本組合作為訓(xùn)練集，另外1個子樣本作為保留集。這樣會獲得k 個預(yù)測方程，記錄k個保留樣本的預(yù)測表現(xiàn)結(jié)果，然后求其平均值。
- bootstrap 包中的 crossval() 函數(shù)可以實現(xiàn) k 重交叉驗證。
- 通過選擇有更好泛化能力的模型，還可以用交叉驗證來挑選變量。
- 下方示例，通過定義shrinkage()函數(shù)對模型的R平方統(tǒng)計量做k 重交叉驗證
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) -
相對重要性
- 最簡單的莫過于比較標(biāo)準(zhǔn)化的回歸系數(shù)，它表示當(dāng)其他預(yù)測變量不變時，該預(yù)測變量一個標(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ù)返回的是一個矩陣，而lm()函數(shù)要求一個數(shù)據(jù)框，你需要用一個中間步驟來轉(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)其他因素不變時，文盲率一個標(biāo)準(zhǔn)差的變化將增加0.68個標(biāo)準(zhǔn)差的謀殺率。根據(jù)標(biāo)準(zhǔn)化的回歸系數(shù)，我們可認(rèn)為Illiteracy是最重要的預(yù)測變量，而Frost是最不重要的- 相對權(quán)重（relative weight）是一種比較有前景的新方法，它是對所有可能子模型添加一個預(yù)測變量引起的R平方平均增加量的一個近似值
- 下方示例提供了一個生成相對權(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") ## 對結(jié)果進(jìn)行查看，可以看到各個預(yù)測變量對模型方差的解釋程度（R平方=0.567），Illiteracy解釋了59%的R平方，F(xiàn)rost解釋了20.79%。根據(jù)相對權(quán)重法，Illiteracy有最大的相對重要性，余下依次是Frost、Population和Income
-
