回归分析-线性回归-检验-模型

2023-11-19 22:17:58

OLS：最小二乘法

  • 通过预测变量的加权和来预测量化的因变量，其中权重是通过数据估计而得的参数
    数据特征：
     正态性 对于固定的自变量值，因变量值成正态分布。
     独立性 Yi值之间相互独立。
     线性 因变量与自变量之间为线性相关。
     同方差性 因变量的方差不随自变量的水平不同而变化。也可称作不变方差

  • 回归模型包含一个因变量和一个自变量时，我们称为简单线性回归

  • 当只有一个预测变量，但同时包含变量的幂（比如，X、X 2、X 3）时，我们称之为多项式回归

  • 当有不止一个预测变量时，则称为多元线性回归

简单线性回归

fit <- lm(weight ~ height, data = women)
summary(fit)
women$weight
fitted(fit)
residuals(fit)

plot(women$height, women$weight, main = "Women Age 30-39", xlab = "Height (in inches)", ylab = "Weight (in pounds)")
abline(fit)

回归分析-线性回归-检验-模型

多项式回归

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

plot(women$height, women$weight, main = "Women Age 30-39", 
    xlab = "Height (in inches)", ylab = "Weight (in lbs)")
lines(women$height, fitted(fit2))

回归分析-线性回归-检验-模型
线性模型： Y ∼ log ⁡ ( x 1 ) + sin ⁡ ( x 2 ) Y \sim \log (x_1)+\sin (x_2) Y∼log(x1​)+sin(x2​)
一般来说，n次多项式生成一个n-1个弯曲的曲线

install.packages("carData")
library(carData)
library(car)
scatterplot(weight ~ height, data = women, spread = FALSE, 
    lty.smooth = 2, pch = 19, main = "Women Age 30-39", xlab = "Height (inches)", 
    ylab = "Weight (lbs.)")

回归分析-线性回归-检验-模型
spread=FALSE选项删除了残差正负均方根在平滑曲线上的展开和非对称信息。lty.smooth=2选项设置loess拟合曲线为虚线。pch=19选项设置点为实心圆（默认为空心圆）。

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

library(car)
scatterplotMatrix(states, spread = FALSE, lty.smooth = 2, 
    main = "Scatterplot Matrix")

回归分析-线性回归-检验-模型

  • 谋杀率是双峰的曲线，每个预测变量都一定程度上出现了偏斜。谋杀率随着人口和文盲率的增加而增加，随着收入水平和结霜天数增加而下降。同时，越冷的州府文盲率越低，收入水平越高。
fit <- lm(Murder ~ Population + Illiteracy + Income + 
    Frost, data = states)

Call:
lm(formula = Murder ~ Population + Illiteracy + Income + Frost, 
    data = states)

Residuals:
    Min      1Q  Median      3Q     Max 
-4.7960 -1.6495 -0.0811  1.4815  7.6210 

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    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

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

回归系数的含义为，一个预测变量增加一个单位，其他预测变量保持不变时，因变量将要增加的数量。例如本例中，文盲率的回归系数为4.14，表示控制人口、收入和温度不变时，文盲率上升1%，谋杀率将会上升4.14%，它的系数在p<0.001的水平下显著不为0。相反，Frost的系数没有显著不为0（p=0.954），表明当控制其他变量不变时，Frost与Murder
不呈线性相关

fit <- lm(mpg ~ hp + wt + hp:wt, data = mtcars)
summary(fit)
install.packages("effects")
library(carData)
library(effects)
plot(effect("hp:wt"),fit,multiline=TRUE)
plot(effect("hp:wt", fit, xlevels=list(wt = c(2.2, 3.2, 4.2))), 
    multiline = TRUE)

每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同
回归分析-线性回归-检验-模型

  • 从图中可以很清晰地看出，随着车重的增加，马力与每加仑汽油行驶英里数的关系减弱了。当wt = 4.2时，直线几乎是水平的，表明随着hp的增加，mpg不会发生改变

回归诊断

  • 回归诊断技术向你提供了评价回归模型适用性的必要工具，它能帮助发现并纠正问题
fit <- lm(Murder ~ Population + Illiteracy + Income +
    Frost, data=states)
confint(fit)

2.5 % 97.5 %
(Intercept) -6.552191e+00 9.0213182149
Population 4.136397e-05 0.0004059867
Illiteracy 2.381799e+00 5.9038743192
Income -1.312611e-03 0.0014414600
Frost -1.966781e-02 0.0208304170

  • 因为Frost的置信区间包含0，可以得出结论说，当其他变量不变时，温度的改变与谋杀率无关

方法

fit <- lm(weight ~ height, data = women)
par(mfrow = c(2, 2))
plot(fit)
par(opar)
  • 正态性 当预测变量值固定时，因变量成正态分布，则残差值也应该是一个均值为0的正态分布。正态Q-Q图（Normal Q-Q，右上）是在正态分布对应的值下，标准化残差的概率图。若满足正态假设，那么图上的点应该落在呈45度角的直线上；若不是如此，那么就违反了正态性的假设。
  • 独立性 你无法从这些图中分辨出因变量值是否相互独立，只能从收集的数据中来验证。
    上面的例子中，没有任何先验的理由去相信一位女性的体重会影响另外一位女性的体重。假若你发现数据是从一个家庭抽样得来的，那么可能必须要调整模型独立性的假设。
  • 线性 若因变量与自变量线性相关，那么残差值与预测（拟合）值就没有任何系统关联。换句话说，除了白噪声，模型应该包含数据中所有的系统方差。在“残差图与拟合图”（Residuals vs Fitted，左上）中可以清楚的看到一个曲线关系，这暗示着你可能需要对回归模型加上一个二次项。
  • 同方差性 若满足不变方差假设，那么在位置尺度图（Scale-Location Graph，左下）中，水平线周围的点应该随机分布。该图似乎满足此假设。
  • “残差与杠杆图”（Residuals vs Leverage，右下）提供了你可能关注的单个观测点

回归分析-线性回归-检验-模型
 一个观测点是离群点，表明拟合回归模型对其预测效果不佳（产生了巨大的或正或负的
残差）。
 一个观测点有很高的杠杆值，表明它是一个异常的预测变量值的组合。也就是说，在预
测变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。
 一个观测点是强影响点（influential observation），表明它对模型参数的估计产生的影响过
大，非常不成比例。强影响点可以通过Cook距离即Cook’s D统计量来鉴别。

newfit <- lm(weight ~ height + I(height^2), data = women)
par(mfrow = c(2, 2))
plot(newfit)

回归分析-线性回归-检验-模型
对于删除数据，要非常小心，因为本应是你的模型去匹配数据，而不是反过来。

除此之外，可以选择更好的检验方法！

正态性

library(car)
fit <- lm(Murder ~ Population + Illiteracy + Income + 
    Frost, data = states)
qqPlot(fit, labels = FALSE, simulate = TRUE, main = "Q-Q Plot")

回归分析-线性回归-检验-模型

  • qqPlot()函数提供了更为精确的正态假设检验方法，它画
    出了在n-p-1个*度的t分布下的学生化残差（studentized residual，也称学生化删除残差或折叠化残差）图形，其中n是样本大小，p是回归参数的数目（包括截距项）

  • id.method = "identify"选项能够交互式绘图——
    待图形绘制后，用鼠标单击图形内的点，将会标注函数中label选项的设定值。敲击Esc键，从图形下拉菜单中选择Stop，或者在图形上右击，都将关闭这种交互模式。此处，我已经鉴定出了Nevada异常。当simulate=TRUE时，95%的置信区间将会用参数自助法

  • 除了Nevada，所有的点都离直线很近，并都落在置信区间内，这表明正态性假设符合得很好。但是你也必须关注下Nevada，它有一个很大的正残差值（真实值—预测值），表明模型低估了该州的谋杀率

residplot <- function(fit, nbreaks=10) {
    z <- rstudent(fit)
    hist(z, breaks=nbreaks, freq=FALSE,
    xlab="Studentized Residual",
    main="Distribution of Errors")
    rug(jitter(z), col="brown")
    curve(dnorm(x, mean=mean(z), sd=sd(z)),
        add=TRUE, col="blue", lwd=2)
    lines(density(z)$x, density(z)$y,
        col="red", lwd=2, lty=2)
    legend("topright",
        legend = c( "Normal Curve", "Kernel Density Curve"),
        lty=1:2, col=c("blue","red"), cex=.7)
}

residplot(fit)

回归分析-线性回归-检验-模型
正如你所看到的，除了一个很明显的离群点，误差很好地服从了正态分布。虽然Q-Q图已经蕴藏了很多信息，但我总觉得从一个柱状图或者密度图测量分布的斜度比使用概率图更容易

误差的独立性

判断因变量值（或残差）是否相互独立，最好的方法是依据收集数据方式的先验知识，时间序列数据通常呈现自相关性——相隔时间越近的观测相关性大于相隔越远的观测。
car包提供了一个可做Durbin-Watson检验的函数，能够检测误差的序列相关性

library(carData)
library(car)
residplot(fit)
durbinWatsonTest(fit)

lag Autocorrelation D-W Statistic p-value
1 -0.2006929 2.317691 0.27
Alternative hypothesis: rho != 0

p值不显著（p=0.282）说明无自相关性，误差项之间独立。滞后项（lag=1）表明数据集中每个数据都是与其后一个数据进行比较的。该检验适用于时间独立的数据，对于非聚集型的数据并不适用。注意，durbinWatsonTest()函数使用自助法来导出p值，如果添加了选项simulate=FALSE，则每次运行测试时获得的结果都将略有不同。

成分残差图

通过成分残差图（component plus residual plot）也称偏残差图（partial residual plot），你可以看看因变量与自变量之间是否呈非线性关系，也可以看看是否有不同于已设定线性模型的系统偏差，图形可用car包中的crPlots()函数绘制。

crPlots(fit, one.page = TRUE, ask = FALSE)

回归分析-线性回归-检验-模型
若图形存在非线性，则说明你可能对预测变量的函数形式建模不够充分，那么就需要添加一些曲线成分，比如多项式项，或对一个或多个变量进行变换（如用log(X)代替X），或用其他回归变体形式而不是线性回归

同方差性

  • ncvTest()函数生成一个计分检验，零假设为误差方差不变，备择假设为误差方差随着拟合值水平的变化而变化。若检验显著，则说明存在异方差性（误差方差不恒定）
  • spreadLevelPlot()函数创建一个添加了最佳拟合曲线的散点图，展示标准化残差绝对值与拟合值的关系

Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 1.746514, Df = 1, p = 0.18632

回归分析-线性回归-检验-模型
计分检验不显著（p=0.19），说明满足方差不变假设。你还可以通过分布水平图看到这一点，其中的点在水平的最佳拟合曲线周围呈水平随机分布。若违反了该假设，你将会看到一个非水平的曲线。代码结果建议幂次变换（suggested power transformation）的含义是，经过p次幂（Y p）变换，非恒定的误差方差将会平稳。例如，若图形显示出了非水平趋势，
建议幂次转换为0.5，在回归等式中用Y 代替Y，可能会使模型满足同方差性。若建议幂次为0，则使用对数变换。对于当前例子，异方差性很不明显，因此建议幂次接近1（不需要进行变换）

多重共线性

  • 那么此处就相当于假定年龄不变，然后测量握力与年龄的关系，这种问题就称作多重共线性（multicollinearity）。它会导致模型参数的置信区间过大，使单个系数解释起来很困难。
  • 多重共线性可用统计量VIF（Variance Inflation Factor，方差膨胀因子）进行检测。VIF的平方根表示变量回归参数的置信区间能膨胀为与模型无关的预测变量的程度。
  • car包中的vif()函数提供VIF值。一般原则下， vif >2就表明存在多重共线性问题。

vif(fit)
Population Illiteracy Income Frost
1.245282 2.165848 1.345822 2.082547
sqrt(vif(fit)) > 2
Population Illiteracy Income Frost
FALSE FALSE FALSE FALSE

异常值

一个全面的回归分析要覆盖对异常值的分析，包括离群点、高杠杆值点和强影响点。
car包也提供了一种离群点的统计检验方法。outlierTest()函数可以求得最大标准化残差绝对值Bonferroni调整后的p值：
library(car)
outlierTest(fit)
rstudent unadjusted p-value Bonferroni p
Nevada 3.542929 0.00095088 0.047544

你可以看到Nevada被判定为离群点（p=0.048）。注意，该函数只是根据单个最大（或正或负）残差值的显著性来判断是否有离群点。若不显著，则说明数据集中没有离群点；若显著，则你必须删除该离群点，然后再检验是否还有其他离群点存在。

高杠杆值观测点

即是与其他预测变量有关的离群点。换句话说，它们是由许多异常的预测变量值组合起来的，与响应变量值没有关系

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)

水平线标注的即帽子均值2倍和3倍的位置。定位函数（locator function）能以交互模式绘图：单击感兴趣的点，然后进行标注，停止交互时，用户可按Esc键退出，或从图形下拉菜单中选择Stop

强影响点

即对模型参数估计值影响有些比例失衡的点
有两种方法可以检测强影响点：Cook距离，或称D统计量，以及变量添加图（added variable plot）。一般来说，Cook’s D值大于4/(n-k -1)，则表明它是强影响点，其中n 为样本量大小，k 是预测变量数目

Cook’s D图有助于鉴别强影响点，但是并不提供关于这些点如何影响模型的信息。变量添加图弥补了这个缺陷。对于一个响应变量和k 个预测变量，你可以如下图创建k 个变量添加图。
回归分析-线性回归-检验-模型

avPlots(fit, ask = FALSE, onepage = TRUE, id.method = "identify")

回归分析-线性回归-检验-模型
利用car包中的influencePlot()函数，你还可以将离群点、杠杆值和强影响点的信息整合到一幅图形中：

influencePlot(fit, id.method = "identify", main = "Influence Plot", 
sub = "Circle size is proportial to Cook's Distance")

回归分析-线性回归-检验-模型
Nevada和Rhode Island是离群点，New York、California、Hawaii和Washington有高杠杆值，Nevada、Alaska和Hawaii为强影响点。

改进措施

  • 不过，我对删除观测点持谨慎态度。若是因为数据记录错误，或是没有遵守规程，或是受试对象误解了指导说明，这种情况下的点可以判断为离群点，删除它们是十分合理的。不过在其他情况下，所收集数据中的异常点可能是最有趣的东西。发掘为何该观测点不同于其他点，有助于你更深刻地理解研究的主题，或者发现其他你可能没有想过的问题。

  • 若Y是比例数，通常使用logit变换[ln (Y/1Y )]。当模型违反了正态假设时，通常可以对响应变量尝试某种变换。car包中的powerTransform()函数通过λ 的最大似然估计来正态化变量X λ

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

MLE of lambda Score Statistic (z)
Population 0.86939 -0.3228
Illiteracy 1.35812 0.6194
Pr(>|z|)
Population 0.7468
Illiteracy 0.5357
iterations = 19

使用变换Population0.87和Illiteracy1.36能够大大改善线性关系。但是对Population（p=0.75）和Illiteracy（p=0.54）的计分检验又表明变量并不需要变换。这些结果与图8-11的成分
残差图是一致的。

如果变换得有意义，比如收入的对数变换、距离的逆变换，
解释起来就会容易得多。但是若变换得没有意义，你就应该避免这样做

增删变量

  • 最常见的方法就是删除某个存在多重共线性的变量
  • 多元回归的变体，专门用来处理多重共线性问题

模型选择

用基础安装中的anova()函数可以比较两个嵌套模型的拟合优度。所谓嵌套模型，即它的一些项完全包含在另一个模型中。

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

anova()函数同时还对是否应该添加Income和Frost到线性模型
中进行了检验。由于检验不显著（p=0.994），因此我们可以得出结论：不需要将这两个变量添加到线性模型中，可以将它们从模型中删除。

AIC（Akaike Information Criterion，赤池信息准则）也可以用来比较模型，它考虑了模型的统计拟合度以及用来拟合的参数数目。AIC值越小的模型要优先选择，它说明模型用较少的参数获得了足够的拟合度

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

Analysis of Variance Table

Model 1: Murder ~ Population + Illiteracy
Model 2: Murder ~ Population + Illiteracy + Income + Frost
Res.Df RSS Df Sum of Sq F Pr(>F)
1 47 289.25
2 45 289.17 2 0.078505 0.0061 0.9939

大量变量的选择

逐步回归法（stepwisemethod）和全子集回归（all-subsets regression）。

逐步回归法的实现依据增删变量的准则不同而不同。MASS包中的stepAIC()函数可以实现逐步回归模型（向前、向后和向前向后），依据的是精确AIC准则：

  • 逐步回归法其实存在争议，虽然它可能会找到一个好的模型，但是不能保证模型就是最佳模型，因为不是每一个可能的模型都被评价了。为克服这个限制，便有了全子集回归法。
  • 全子集回归可用leaps包中的regsubsets()函数实现。你能通过R平方、调整R平方或Mallows Cp统计量等准则来选择“最佳”模型。
library(leaps)
leaps <- regsubsets(Murder ~ Population + Illiteracy + 
    Income + Frost, data = states, nbest = 4)
plot(leaps, scale = "adjr2")

回归分析-线性回归-检验-模型
第一行中（图底部开始），可以看到含intercept（截距项）和Income的模型调整R平方为0.33，含intercept和Population的模型调整R平方为0.1。跳至第12行，你会看到
含intercept、Population、Illiteracy和Income的模型调整R平方值为0.54，而仅含intercept、Population和Illiteracy的模型调整R平方为0.55。此处，你会发现含预测变量越少的模型调整R平方越大（对于非调整的R平方，这是不可能的）。图形表明，双预测变量模型（Population和Illiteracy）是最佳模型。

回归分析-线性回归-检验-模型

library(car)
subsets(leaps, statistic = "cp", 
    main = "Cp Plot for All Subsets Regression")
abline(1, 1, lty = 2, col = "red")

回归分析-线性回归-检验-模型

  • 你会看到对于不同子集大小，基于Mallows Cp统计量的四个最佳模型。越好的模型离截距项和斜率均为1的直线越近。

  • 变量自动选择应该被看做是对模型选择的一种辅助方法，
    而不是直接方法。拟合效果佳而没有意义的模型对你毫无帮助，主题背景知识的理解才能最终指
    引你获得理想的模型。

上一篇：#7：整数反转


下一篇：LeetCode刷题 7. 整数反转