其他分享
首页 > 其他分享> > 回归分析-线性回归-检验-模型

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

作者：互联网

OLS：最小二乘法

简单线性回归

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)

每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同
在这里插入图片描述

回归诊断

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

方法

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

在这里插入图片描述
 一个观测点是离群点，表明拟合回归模型对其预测效果不佳（产生了巨大的或正或负的
残差）。
 一个观测点有很高的杠杆值，表明它是一个异常的预测变量值的组合。也就是说，在预
测变量空间中，它是一个离群点。因变量值不参与计算一个观测点的杠杆值。
 一个观测点是强影响点（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")

在这里插入图片描述

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），或用其他回归变体形式而不是线性回归

同方差性

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（不需要进行变换）

多重共线性

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为强影响点。

改进措施

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准则：

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

在这里插入图片描述

标签：变量,fit,方差,回归,残差,检验,线性,模型,Population
来源： https://blog.csdn.net/qq_40243662/article/details/122726671