title: "lesson 8 回歸"
author: "wintryheart"
date: "2019/9/9"
output:
html_document:
toc: TRUE
number_sections: TRUE

knitr::opts_chunk$set(echo = TRUE)

研究的三個(gè)困難

• 發(fā)現(xiàn)一個(gè)有趣的問(wèn)題
• 設(shè)計(jì)一個(gè)有用的、可以測(cè)量的響應(yīng)變量
• 收集合適的數(shù)據(jù)

OLS回歸

統(tǒng)計(jì)假設(shè)：

• 正態(tài)性：對(duì)于固定的自變量值，因變量值成正態(tài)分布；
• 獨(dú)立性：值之間相互獨(dú)立；
• 線性
• 同方差性： 因變量的方差不隨自變量的水平不同而變化。

lm()擬合模型

簡(jiǎn)單線性回歸

fit <- lm(weight ~ height, data=women)

summary(fit)

fitted(fit)

residuals(fit)

plot(women$height, women$weight)

abline(fit)

unnamed-chunk-1-1.png

多項(xiàng)式回歸

fit2 <- lm(weight ~ height + I(height^2), data=women)

summary(fit2)

• car包scatterplot()函數(shù)繪制二元關(guān)系圖

library(car)
scatterplot(weight~height,data=women, spread=F, smoother.args=list(lty=2), pch=19,
            main="Women Age 30-39", xlab="hight", ylab="weight")

unnamed-chunk-3-1.png

spread=FALSE，刪除殘差正負(fù)均方根在平滑曲線上的展開(kāi)和非對(duì)稱信息。

smoother.args=list(lty=2),設(shè)置loess擬合曲線為虛線。

多元回歸

• car包scatterplotMatrix()生成散點(diǎn)圖矩陣。
• scatterplotMatrix()默認(rèn)在非對(duì)角線區(qū)域繪制變量間散點(diǎn)圖，并添加平滑和線性擬合曲線。
• 對(duì)角線區(qū)域繪制每個(gè)變量的密度圖和軸須圖。

states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])

cor(states)

library(car)

scatterplotMatrix(states, spread=FALSE, smoother.args=list(lty=2), main="Scatter Plot Matrix")

fit3 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)

summary(fit3)

unnamed-chunk-4-1.png

交互項(xiàng)

fit4 <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)

summary(fit4)

• effects包專門用于顯示線性模型、一般線性模型和其他許多模型的作用效應(yīng)。
• effects包中的effect()函數(shù)，可以用圖形展示線性模型的交互項(xiàng)的結(jié)果。
• plot(effect(term, mod, ,xlevels), multiline=TRUE)

library(effects)
plot(effect("hp:wt", fit4, , list(wt=c(2.2, 3.2, 4.2))))
plot(effect("hp:wt", fit4, , list(wt=c(2.2, 3.2, 4.2))), multiline=TRUE)

unnamed-chunk-6-1.png

unnamed-chunk-6-2.png

回歸診斷

# 置信區(qū)間
confint(fit3)

結(jié)果表明，文盲率改變1%， 謀殺率在95%的置信區(qū)間[2.38, 5.90]中變化。

因?yàn)镕rost的置信區(qū)間包含0， 所以可以得出結(jié)論：當(dāng)其他變量不變時(shí)，溫度的改變與謀殺率無(wú)關(guān)。

標(biāo)準(zhǔn)方法

對(duì)lm()函數(shù)返回的對(duì)象使用plot()函數(shù)，可以生成評(píng)價(jià)模型擬合情況的四幅圖形。

• 正態(tài)Q-Q圖，檢驗(yàn)正態(tài)性。如果滿足正態(tài)假設(shè)，圖上的點(diǎn)應(yīng)該落在呈45度角的直線上。
• 殘差圖與擬合圖，檢驗(yàn)線性。如果因變量與自變量線性相關(guān)，殘差與預(yù)測(cè)值應(yīng)該沒(méi)有任何系統(tǒng)關(guān)聯(lián)。
• 位置尺度圖（Scale-Location），檢驗(yàn)同方差性。如果滿足不變方差假設(shè)，水平線周圍的點(diǎn)應(yīng)該隨機(jī)分布。
• 殘差與杠桿圖（Residuals VS Leverage），檢驗(yàn)離群點(diǎn)、杠桿值和強(qiáng)影響點(diǎn)。

fit <- lm(weight ~ height, data=women)
par(mfrow=c(2,2))
plot(fit)

#殘差擬合圖顯示一個(gè)曲線關(guān)系，暗示需要對(duì)回歸模型加上一個(gè)二次項(xiàng)。
fit2 <- lm(weight ~ height + I(height^2), data = women)
plot(fit2)

fit3 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
plot(fit3)

unnamed-chunk-8-1.png

unnamed-chunk-8-2.png

unnamed-chunk-8-3.png

改進(jìn)方法

car包中的回歸診斷實(shí)用函數(shù)

正態(tài)性檢驗(yàn)qqPlot()

states <- as.data.frame(state.x77[,c("Murder", "Population","Illiteracy","Income","Frost")])
fit3 <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
par(mfrow=c(1,1))
qqPlot(fit3, labels=row.names(states), simulate=TRUE, main="Q-Q Plot")
states["Nevada",]
fitted(fit3)["Nevada"]
residuals(fit3)["Nevada"]
rstudent(fit3)["Nevada"]

unnamed-chunk-9-1.png

• Default S3 method:

qqPlot(x, distribution="norm", groups, layout,
ylim=range(x, na.rm=TRUE), ylab=deparse(substitute(x)),
xlab=paste(distribution, "quantiles"), glab=deparse(substitute(groups)),
main=NULL, las=par("las"),
envelope=.95, col=carPalette()[1], col.lines=carPalette()[2],
lwd=2, pch=1, cex=par("cex"),
line=c("quartiles", "robust", "none"), id=TRUE, grid=TRUE, ...)

• S3 method for class 'formula'

qqPlot(formula, data, subset, id=TRUE, ylab, glab, ...)

• S3 method for class 'lm'

qqPlot(x, xlab=paste(distribution, "Quantiles"), ylab=paste("Studentized Residuals(",deparse(substitute(x)), ")", sep=""),
main=NULL, distribution=c("t", "norm"),
line=c("robust", "quartiles", "none"), las=par("las"),
simulate=TRUE, envelope=.95, reps=100,
col=carPalette()[1], col.lines=carPalette()[2], lwd=2, pch=1, cex=par("cex"),
id=TRUE, grid=TRUE, ...)

distribution:
root name of comparison distribution – e.g., "norm" for the normal distribution; t for the t-distribution.

envelope:
confidence level for point-wise confidence envelope, or FALSE for no envelope. 置信區(qū)間

id:
controls point identification; if FALSE, no points are identified; can be a list of named arguments to the showLabels function; TRUE is equivalent to list(method="y", n=2, cex=1, col=carPalette()[1], location="lr"), which identifies the 2 points with the 2 points with the most extreme verical values — studentized residuals for the "lm" method. Points labels are by default taken from the names of the variable being plotted is any, else case indices are used. Unlike most graphical functions in car, the default is id=TRUE to include point identification.

simulate:
if TRUE calculate confidence envelope by parametric bootstrap; for lm object only. The method is due to Atkinson (1985). 參數(shù)自助法生成置信區(qū)間

誤差的獨(dú)立性檢驗(yàn)，Durbin-Watson檢驗(yàn)

durbinWatsonTest(fit3)

p值不顯著，說(shuō)明無(wú)自相關(guān)性，誤差項(xiàng)之間獨(dú)立。

線性檢驗(yàn)，crPlots()

通過(guò)成分殘差圖（component plus residual plot），或偏差圖（partial residual plot），看因變量和自變量之間是否呈非線性關(guān)系。

crPlots(fit3)

unnamed-chunk-11-1.png

成分殘差圖說(shuō)明，線性模型形式是合適的。

若圖形存在非線性，說(shuō)明線性建模不夠充分，需要添加一些曲線成分（如多項(xiàng)式），或?qū)ψ兞窟M(jìn)行變換（如log變換），或用其它回歸變體形式而不是線性回歸。

同方差性檢驗(yàn)，判斷誤差方差是否恒定

ncvTest()生成一個(gè)計(jì)分檢驗(yàn)，零假設(shè)為誤差方差不變。若檢驗(yàn)顯著，說(shuō)明存在異方差性（誤差方差不恒定）。

spreadLevelPlot()創(chuàng)建一個(gè)最佳擬合曲線的散點(diǎn)圖，展示標(biāo)準(zhǔn)化殘差絕對(duì)值與擬合值的關(guān)系。

ncvTest(fit3)
spreadLevelPlot(fit3)

unnamed-chunk-12-1.png

ncvTest()計(jì)分檢驗(yàn)不顯著，說(shuō)明 滿足同方差性假設(shè)。

spreadLevelPlot()建議冪次變換（suggested power transformation)，即經(jīng)過(guò)p次冪變換（），非恒定的誤差方差將會(huì)平穩(wěn)。。由于建議冪次接近1，異方差性不明顯，不需要進(jìn)行變換。

線性模型假設(shè)的綜合驗(yàn)證

gvlma包中的gvlma()函數(shù)，對(duì)線性模型假設(shè)進(jìn)行綜合驗(yàn)證，同時(shí)還能做偏斜度、峰度和異方差性的評(píng)價(jià)。

library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

Decision列表明是否違反假設(shè)條件。

多重共線性

多重共線性用統(tǒng)計(jì)量VIF（variance inflation factor， 方差膨脹因子）進(jìn)行檢測(cè)。

VIF的平方根表示，變量回歸參數(shù)的置信區(qū)間能膨脹為與模型無(wú)關(guān)的預(yù)測(cè)變量的程度。

car包中的vif()函數(shù)提供VIF值。一般原則下，就表明存在多重共線性問(wèn)題。

library(car)

vif(fit3)

sqrt(vif(fit3))>2

異常值

離群點(diǎn)

離群點(diǎn)指那些模型預(yù)測(cè)不佳的觀測(cè)點(diǎn)。

• Q-Q圖中，落在置信區(qū)間帶外的點(diǎn)，即可被認(rèn)為是離群點(diǎn)。
• 標(biāo)準(zhǔn)化殘差的絕對(duì)值大于2的點(diǎn)可能是離群點(diǎn)。
• car包的outlierTest()函數(shù)，根據(jù)單個(gè)絕對(duì)值最大的殘差值的顯著性來(lái)判斷是否有離群點(diǎn)。

library(car)
fit5 <- lm(mpg~ wt+qsec, data = mtcars)
outlierTest(fit5)

高杠桿值點(diǎn)

高杠桿值點(diǎn)指那些與其他預(yù)測(cè)變量有關(guān)的離群點(diǎn)，即它們是許多異常的預(yù)測(cè)變量值組合起來(lái)的，與響應(yīng)變量的值無(wú)關(guān)。

高杠桿值點(diǎn)可以通過(guò)帽子統(tǒng)計(jì)量(hat statistic)判斷。帽子均值為，其中p為模型估計(jì)的參數(shù)數(shù)目（包含截距項(xiàng)），n是樣本量。
如果觀測(cè)點(diǎn)的帽子值大于帽子均值的2或3倍，就可以認(rèn)定為高杠桿值點(diǎn)。

hatvalues(fit3)

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(x=1:n, y=hatvalues(fit), labels=names(hatvalues(fit)))
}
hat.plot(fit3)

unnamed-chunk-16-1.png

強(qiáng)影響點(diǎn)

強(qiáng)影響點(diǎn)指對(duì)模型參數(shù)估計(jì)值影響有些比例失衡的點(diǎn)。即移除該觀測(cè)點(diǎn)，模型會(huì)發(fā)生重大改變。

有兩種方法可以檢測(cè)：Cook距離或稱D統(tǒng)計(jì)量，以及變量添加圖（added variable plot）。

Cook距離

一般來(lái)說(shuō)，Cook's D值大于4/(n-k-1)則表明它是強(qiáng)影響點(diǎn)，其中n為樣本量，k為預(yù)測(cè)變量數(shù)。

cutoff <- 4/(nrow(states)-length(fit3$coefficients)-2)
plot(fit3, which = 4, cook.levels = cutoff)
abline(h=cutoff, lty=2, col="red")

unnamed-chunk-17-1.png

其實(shí)以1為分割點(diǎn)，比以4/(n-k-1)為分割點(diǎn)更具有一般性。

變量添加圖

變量添加圖，即對(duì)于每個(gè)預(yù)測(cè)變量，繪制在其他k-1個(gè)預(yù)測(cè)變量上回歸的殘差值相對(duì)于響應(yīng)變量在其他k-1預(yù)測(cè)變量上的回歸殘差值的關(guān)系圖。

car包中的avPlots()函數(shù)，可以作變量添加圖。

library(car)

avPlots(fit3, ask=FALSE, id.method="identify")

unnamed-chunk-18-1.png

car包中的influencePlot()函數(shù)，可以把離群點(diǎn)、杠桿值和強(qiáng)影響點(diǎn)的信息整合到一幅圖形中。

influencePlot(fit3, id.method="identify", main="Influence Plot",
              sub="Circle size is proportional to Cook's distance")

unnamed-chunk-19-1.png

其中，

• 縱坐標(biāo)超過(guò)+2或小于-2的點(diǎn)，可被認(rèn)為離群點(diǎn)。
• 水平坐標(biāo)超過(guò)0.2或0.3的點(diǎn)，可被認(rèn)為高杠桿值。
• 圓圈大小與影響成比例，圓圈很大的點(diǎn)可能是強(qiáng)影響點(diǎn)。

改進(jìn)措施

刪除觀測(cè)點(diǎn)

變量變換

car包

• powerTransform()通過(guò)的最大似然估計(jì)來(lái)正態(tài)化變量。
• boxTidwell()通過(guò)獲得預(yù)測(cè)變量?jī)鐢?shù)的最大似然估計(jì)來(lái)改善線性關(guān)系。
• spreadLevelPlot()提供冪次變換應(yīng)用。

library(car)
summary(powerTransform(states$Murder))

結(jié)果表明，可用來(lái)正態(tài)化變量Murder；不過(guò)的假設(shè)也無(wú)法拒絕。因此沒(méi)有強(qiáng)有力的證據(jù)表明需要變換變量。

boxTidwell(Murder ~ Population + Illiteracy, data=states)

結(jié)果顯示，使用和能夠改善線性關(guān)系。但是檢驗(yàn)結(jié)果又表明變量并不需要變換。

如果變量變換得沒(méi)有意義，應(yīng)該避免這樣做。

增刪變量

嘗試其他方法

• 如果存在離群點(diǎn)或強(qiáng)影響點(diǎn)，可以使用穩(wěn)健回歸；
• 如果違背了正態(tài)性假設(shè)，可以使用非參數(shù)回歸；
• 如果存在顯著的非線性，可以嘗試非線性回歸；
• 如果違背誤差獨(dú)立性假設(shè)，可以嘗試專門研究誤差結(jié)構(gòu)的模型，如時(shí)間序列模型或多層回歸模型。
• 廣義線性回歸模型，適用于許多OLS回歸假設(shè)不成立的情況。

選擇最佳模型

模型比較

anova()函數(shù)可以比較兩個(gè)嵌套模型的擬合優(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)

結(jié)果表明，模型預(yù)測(cè)效果差異并不顯著，Income和Frost兩個(gè)變量沒(méi)必要添加到模型中。

AIC()函數(shù)可以比較兩個(gè)模型的AIC。AIC值越小越好。

AIC(fit1, fit2)

注意，anova()需要嵌套模型，而AIC()不需要。

變量選擇

anova()和AIC()只能比較兩個(gè)模型。
如果有更多的模型，需要用以下方法。

逐步回歸

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")

逐步回歸法可能會(huì)找到一個(gè)好的模型，但是不能保證模型就是最佳的，因?yàn)椴皇敲恳粋€(gè)可能的模型都被評(píng)價(jià)。

全子集回歸

leaps包中的regsubsets()函數(shù)，所有可能的模型都會(huì)被檢驗(yàn)。

通過(guò)R方，調(diào)整R方或Mallows Cp統(tǒng)計(jì)量等準(zhǔn)則選擇最佳模型。

對(duì)于一個(gè)好的模型，它的Cp統(tǒng)計(jì)量非常接近于模型的參數(shù)數(shù)目（包括截距項(xiàng)）。

結(jié)果可用leaps包中的plot()繪制，也可用car包中的subsets()繪制。

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", legend=c(3.5, 50), main="Cp Plot for All Subsets Regression")
abline(1, 1, lty=2, col="red")

subsets(leaps, statistic="bic", legend=c(3.5, 0), main="Cp Plot for All Subsets Regression")
subsets(leaps, statistic="adjr2", legend=c(3.5, 0.1), main="Cp Plot for All Subsets Regression")

unnamed-chunk-25-1.png

unnamed-chunk-25-2.png

unnamed-chunk-25-3.png

unnamed-chunk-25-4.png

在subsets()繪制的Cp圖中，越好的模型離截距項(xiàng)和斜率均為1的直線越近。

深層次分析：相對(duì)重要性

最簡(jiǎn)單的方法是比較標(biāo)準(zhǔn)化的回歸系數(shù)。

利用scale()函數(shù)將數(shù)據(jù)標(biāo)準(zhǔn)化為均值為0，標(biāo)準(zhǔn)差為1的數(shù)據(jù)集。
（注意，scale()返回的是矩陣，而lm()要求的是數(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 + Illiteracy + Income + Frost, data=zstates)

coef(zfit)

round(coef(zfit), digits = 2) #系數(shù)保留兩位小數(shù)

從結(jié)果可以看出Illiteracy的標(biāo)準(zhǔn)化系數(shù)最大（r round(coef(zfit), digits = 2)[3]），可以被認(rèn)為是最重要的變量。

還有其它方法可以定量分析相對(duì)重要性，如可以將相對(duì)重要性看作每個(gè)預(yù)測(cè)變量（本身或與其他預(yù)測(cè) 變量組合）對(duì)R方的貢獻(xiàn)。relaimpo包涵蓋了一些相對(duì)重要性的評(píng)價(jià)方法。

library(relaimpo)

linmod <- lm(Fertility~Agriculture+Examination+Education+Catholic+Infant.Mortality,swiss)
crlm <- calc.relimp(linmod, type = c("lmg", "last", "first", "betasq", "pratt", "genizi", "car"), rela = TRUE )
crlm
plot(crlm)

unnamed-chunk-27-1.png

calc.relimp()計(jì)算線性模型的幾種相對(duì)重要性指標(biāo)。推薦使用lgm。該方法將R2。

• lmg 將R2分解為每個(gè)自變量的貢獻(xiàn)，參見(jiàn) Lindeman, Merenda and Gold 1980, p.119ff ，Chevan and Sutherland (1991)和《線性回歸模型中自變量相對(duì)重要性的衡量》(孫紅衛(wèi)，王玖，羅文海,2012)

• pmvd 是比例邊際方差分解，可以解釋為帶樣本權(quán)重的加權(quán)自變量貢獻(xiàn)。參見(jiàn)Feldman (2005)

• last 是每一個(gè)變量最后被包括進(jìn)來(lái)時(shí)的貢獻(xiàn)，有時(shí)也稱有用度。

• first 是每一個(gè)變量最初被包括進(jìn)來(lái)時(shí)的貢獻(xiàn)，也即是y和該自變量的協(xié)方差平方。

• betasq 是標(biāo)準(zhǔn)化系數(shù)的平方

• pratt 是標(biāo)準(zhǔn)化系數(shù)和相關(guān)系數(shù)的乘積

• genizi 是根據(jù)Genizi 1993的R2分解

• car 是根據(jù)Zuber and Strimmer 2010的R2分解。

calc.relimp 中有五個(gè)指標(biāo)(lmg, pmvd, pratt, genizi and car)分解R2。當(dāng)rela=FALSE，它們的值相加等于R2；當(dāng)rela=TRUE，它們的值相加等于100(%)