残差平方和ssr的计算公式为_手把手教你用R语言做回归后的残差分析(附代码)...

本文详细介绍了如何在R语言中进行残差分析,探讨了残差的本质和重要性,并提供了计算和可视化残差的方法。通过分析残差的分布、偏态和峰度,以及利用Durbin-Watson检验检测残差间的相关性,有助于评估线性回归模型的可靠性。此外,文章还展示了如何识别并处理异常值,以及理解权重在回归分析中的作用。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

cdfb5bea0bf02c80ad296b7ca555d4d3.png
作者:Abhijit Telang翻译:张睿毅校对:丁楠雅本文约 2600字,建议阅读 10分钟。

本文介绍了做残差分析的方法及其重要性,以及利用R语言实现残差分析。

在这篇文章中,我们通过探索残差分析和用R可视化结果,深入研究了R语言。

残差本质上是当一个给定的模型(在文中是线性回归)不完全符合给定的观测值时留下的gap。

医学中的病理学发现的残留分析是一个形象的比喻。人们通常用代谢残留水平来作为衡量药物吸收的指标。

残差是用于建模的原始值与作为模型结果的对于原始值的估计之间的差异。

残差=y-y-hat,其中y是初始值,y-hat是计算值。

期望这个错误尽可能接近于零,并且通过残差找到任何异常值。

找到异常值的一个快速方法是使用标准化残差。第一种方法是简单地求出残差与其标准差的比值,因此,任何超过3个标准差的情况都可以被视为异常值。

## 标准化残差-相对于其标准偏差的比例残差

residueStandard

df[residueStandard>3,]

以下是得到的结果:

days.instant days.atemp days.hum days.windspeed days.casual

442 442 0.505046 0.755833 0.110704 3155

456 456 0.421708 0.738333 0.250617 2301

463 463 0.426129 0.254167 0.274871 3252

470 470 0.487996 0.502917 0.190917 2795

471 471 0.573875 0.507917 0.225129 2846

505 505 0.566908 0.456250 0.083975 3410

512 512 0.642696 0.732500 0.198992 2855

513 513 0.641425 0.697083 0.215171 3283

533 533 0.594708 0.504167 0.166667 2963

624 624 0.585867 0.501667 0.247521 3160

645 645 0.538521 0.664167 0.268025 3031

659 659 0.472842 0.572917 0.117537 2806

当然,我希望我的模型是无偏的,至少我想这样。因此回归线两边的任何残差,如果没有在这条线上,都是随机的,也就是说,没有任何特定的模式。

也就是说,我希望我的剩余误差分布遵循一个普通的正态分布。

使用R语言,只需两行代码就可以优雅地完成这项工作。

  • 绘制残差柱状图;
  • 添加一个分位数图,其中有一条线穿过,即第一个和第三个分位数。

hist(lmfit$residuals)

qqnorm(lmfit$residuals);qqline(lmfit$residuals)

ff607dc8d0d0e52d33d17d33ea2d179a.png

于是,我们知道这个图偏离了正常值(正常值用直线表示)。

但这种非黑即白的信息一般是不够的。因此,我们应该检查偏态峰度,以了解分布的分散性。

首先,我们将计算偏态;我们将使用一个简单的高尔顿偏态(Galton’s skewness)公式。

## 分布对称性检验:偏态

summary

Q1

Q2

Q3

skewness

skewness

[1] 0.2665384

这量化了分布的方向(分布于右侧或左侧,或者是否完全对称或均匀)。

另一个度量角度是峰度,它显示分布是朝向中心(+ve值)还是远离中心(-ve值)。峰度是一个量化离群值可能性的指标。正值表示离群的存在。

kurtosis(lmfit$residuals)

[1] 2.454145

接下来,除了可视化残差的分布,我还想找出残差是否相互关联。

对于一个模型来说,为了解释观测值的所有变化,残差必须随机发生,并且彼此不相关。

Durbin-Watson测试允许检验残差彼此之间的独立性。

dwt(lm(df$days.windspeed~df$days.atemp+df$days.hum),method='resample', alternative='two.sided')

lag Autocorrelation D-W Statistic p-value

1 0.2616792 1.475681 0

Alternative hypothesis: rho != 0

对于给定的自由度和观测次数,需要将统计值与临界值表确定的下限和上限进行比较。文中案例的值域是[1.55,1.67]。

由于计算的D-W统计值低于该范围的较低值,我们拒绝了残差不相关的零假设。因此取而代之的假设是残差之间可能存在相关性。

直观地看,这个假设可以通过研究模型在试图捕获原始Y值的增加值时失败的原因来了解。当捕获增加值时,随着y的增加,残差与y成正比。

0569ca22fcfc308272cc68f751397245.png

将其与绘制拟合y-hat值与y值进行比较。当y-hat值趋于落后时,残差似乎与y共同增长,故此,因为过去的残值似乎继续沿着固定的坡度值运行,过去的残值似乎是当前值的更好预测因子。同时,在达尔文-沃森检验(Darwin-Watson tests)中在残差与先前值之间的差的平方和,与所有观测的给定残差之和的比较和对比中,发现了相关性。

因此,残余误差之间的相关性就像心率测试失灵一样。如果你的模型不能跟上Y的快速变化,它会越来越和一个检验玩命的跑步者的情况类似,没记录上的步数似乎比跑步者实际的步数更相关。

afacdf75394d5ec82019a3c717dec277.png

重构你的线性回归模型中的下意识影响。

有一点很重要:当对任何给定的观测集进行线性回归时,因变量(符号上表示为y)的计算估计量(符号上表示为y-hat)的每个值不仅依赖于当前值(例如,观察值),还依赖于每次观测。

毕竟,你的多元线性回归模型只在最小化因变量y的实际值和所有观测值y(y-hat)的计算估计值之间的误差后,才计算出每个影响因素的系数。

在数学上,这可以用简单的加权平均模型表示,如下所示。

Yhat(j)= w1j*Y(1)+w2j*Y(2)+w3j*Y(3)+w4j*Y(4)+.....+wnj*Y(n);

where n is number of observations and j is one particular estimated value.

如果你想用矩阵表示法来考虑这个问题,矩阵w的维数是[n,n],其中n是观察到的数目。同样,以矩阵形式表示的y具有尺寸[n,1],因此y-hat[j]表示为矩阵,因为它是两个矩阵(w和y)的乘积,并且具有与y相同的尺寸[n,1]。

因此,计算出的估计量的每一个值都可以表示为所有初始y值的加权平均值,表明给定的计算估计量在多大程度上受到每个初始y值的影响。

更简单地说,如果我看到一堆分散的观察结果,我必须以尽可能小的偏差,画一条直线穿过它们,也就是说,如果我必须遗漏一些观察结果,我必须在这条回归线的两侧平等地遗漏。

bc6b4b46dbd4a7c402c65e16cff22bd5.png

但是在这样做的同时,我(或选择的回归程序)不管来源何处,下意识地估计了每次观察对我回归线的斜率的影响程度。

我现在所做的就是重新计算这些影响。为了改变我的回归线,将它们包括在内或排除在外,每一个观察值的权重是多少?

权重的统计术语是hat values,因为它们连接了计算的估计量和它们的原始对应的估计量。

以下是用R语言计算的方法:

modelmatrix

hatvalues

首先,我们得到一个矩阵形式的模型。然后我们计算hat值。

或者,可以使用以下函数获得类似的结果。

hatvalues

让我们考虑一下可以施加在每个权重上的限制。

显然,权重的最小可能值等于所有原始Y值贡献相等的可能性(因为它们必须为线性回归程序贡献一些东西,通过对所有观测进行工作和优化来估计系数)。

在这种情况下,其值域的下限为1/n,其中n是观测总数。在任何情况下,因为n总是+ve,所以权重总是<1

现在试着将hat值加和,你会看到一个有趣的结果…

sum(hatvalues)

[1] 4

它们等于线性回归模型为计算考虑的影响因素数量+1。

例如,在示例数据集中,我们有三个因素,即温度、湿度和风速。

c17463c90038ce1512cf46ba9e46b6fb.png

接下来,我们如何找到最重要或最有影响的观察结果?

一种优雅的方式是:

  • 将hat值切分为四分位数。
  • 应用95%标准过滤最异常值。
  • 将该过滤标准应用于观察结果。

R语言允许你一步完成这些操作!

df[hatvalues>quantile(hatvalues,0.95),]

我们可以得到满足这一标准的具体观测结果。

days.instant days.atemp days.hum days.windspeed days.casual

9 9 0.1161750 0.434167 0.361950 54

21 21 0.1578330 0.457083 0.353242 75

22 22 0.0790696 0.400000 0.171970 93

26 26 0.2036000 0.862500 0.293850 34

32 32 0.2345300 0.829565 0.053213 47

36 36 0.2430580 0.929167 0.161079 100

45 45 0.3983500 0.375833 0.417908 208

50 50 0.3914040 0.187917 0.507463 532

65 65 0.3662520 0.948261 0.343287 114

69 69 0.3856680 0.000000 0.261877 46

90 90 0.2575750 0.918333 0.217646 179

94 94 0.5429290 0.426250 0.385571 734

95 95 0.3983500 0.642083 0.388067 167

106 106 0.4254920 0.888333 0.340808 121

153 153 0.6439420 0.305000 0.292287 736

239 239 0.6355560 0.850000 0.375617 226

249 249 0.5152000 0.886957 0.343943 204

293 293 0.4665250 0.636250 0.422275 471

302 302 0.2279130 0.882500 0.351371 57

341 341 0.4002460 0.970417 0.266175 50

368 368 0.1262750 0.441250 0.365671 89

383 383 0.2752540 0.443333 0.415429 109

388 388 0.2430580 0.911250 0.110708 145

408 408 0.1016580 0.464583 0.409212 73

421 421 0.2556750 0.395833 0.421642 317

433 433 0.5246040 0.567500 0.441563 486

434 434 0.3970830 0.407083 0.414800 447

463 463 0.4261290 0.254167 0.274871 3252

465 465 0.4766380 0.317500 0.358196 905

478 478 0.3895040 0.835417 0.344546 120

543 543 0.5947040 0.373333 0.347642 1077

627 627 0.5650670 0.872500 0.357587 371

667 667 0.4677710 0.694583 0.398008 998

668 668 0.4394000 0.880000 0.358200 2

694 694 0.2487420 0.404583 0.376871 532

722 722 0.2361130 0.441250 0.407346 205

726 726 0.2203330 0.823333 0.316546 9

最终,你可以用颜色编码使他们具象化。

##具象化

plot(df$days.casual,lmfit$fitted.values,col=ifelse(hatvalues>quantile(hatvalues,0.95),'red','blue'))

1cf5453b9df8d42af426c1847d05df87.png

就是这样!这是进行残差分析和其重要性的简单概述。

原文标题:

Doing Residual Analysis Post Regression in R

原文链接:

https://dzone.com/articles/doing-residual-analysis-post-regression-in-r

编辑:王菁

校对:龚力

译者简介

9ede5fb231f72df98360810e3323a2fe.png

张睿毅,北京邮电大学大二物联网在读。我是一个爱自由的人。在邮电大学读第一年书我就四处跑去蹭课,折腾整一年惊觉,与其在当下焦虑,不如在前辈中沉淀。于是在大二以来,坚持读书,不敢稍歇。资本主义国家的科学观不断刷新我的认知框架,同时因为出国考试很早出分,也更早地感受到自己才是那个一直被束缚着的人。太多真英雄在社会上各自闪耀着光芒。这才开始,立志终身向遇到的每一个人学习。做一个纯粹的计算机科学里面的小学生。

— 完 —

关注清华-青岛数据科学研究院官方微信公众平台“THU数据派”及姊妹号“数据派THU”获取更多讲座福利及优质内容。

<think>好的,我需要帮助用户计算sm.OLS.fit()结果的残差方差。首先,我得回忆一下statsmodels库的相关知识。用户提到了引用中的内容,特别是引用[2]中的回归结果参数输出。残差方差通常指的是残差的方差,或者是误差项的估计值,比如σ²(sigma squared)。 首先,我应该确定残差方差的定义。在线性回归中,残差方差通常用均方误差(MSE)来表示,计算公式残差平方和除以自由度,自由度是样本数减去自变量数目减一。也就是MSE = SSR / (n - k - 1),其中SSR残差平方和,n是样本量,k是自变量个数。 接下来,我需要查看statsmodels的OLS结果对象有哪些属性。根据引用[2],回归结果存储在model变量中,调用model.summary()可以输出结果。可能model对象有resid属性,即残差,然后计算残差平方和,再除以自由度。或者,是否已经有现成的属性可以直接获取残差方差? 查阅statsmodels文档,OLSResults对象有一个属性叫mse_resid,这正是残差方差的无偏估计,即MSE。所以用户可以直接用model.mse_resid来获取。但为了确认,可能用户需要手动计算,比如用残差平方和除以自由度。 此外,可能需要展示两种方法:一种是直接调用属性,另一种是手动计算。比如,残差可以通过model.resid获得,然后计算平方和,再除以自由度。自由度可以通过model.df_resid得到,即残差自由度n - k - 1。 所以步骤应该是: 1. 拟合模型,得到model。 2. 获取残差:model.resid 3. 计算残差平方和ssr = np.sum(model.resid ** 2) 4. 自由度:df = model.df_resid 5. 残差方差:ssr / df 同时指出也可以用model.mse_resid直接获取。 需要确保用户正确导入库,比如numpy,并且正确拟合了模型。例如,在引用[2]中,回归是通过sm.OLS(Y,X)然后fit()得到的。所以代码示例应该包括这些步骤,并展示两种方法。 另外,注意残差方差和残差标准误差的区别,后者是方差的平方根,即model.scale的平方根?或者可能model.scale就是残差方差?需要核实。根据statsmodels文档,model.scale是残差方差的无偏估计,等于ssr / df_resid,而model.mse_resid和model.scale是同一个值。所以用户可以使用model.scale或者model.mse_resid。 所以最终回答应该分步骤,说明两种方法,并指出相关属性,同时提醒注意自由度的计算是否正确。然后生成相关问题,比如如何计算残差标准误差,如何获取残差,或者如何手动验证结果。</think>要计算statsmodels中OLS回归结果的残差方差(即误差项的无偏估计量$\hat{\sigma}^2$),可按照以下步骤操作: --- ### **步骤说明** 1. **拟合OLS模型** ```python import statsmodels.api as sm import numpy as np # 假设已有自变量X和因变量Y(需为DataFrame或数组) X = sm.add_constant(X) # 添加常数项 model = sm.OLS(Y, X).fit() # 拟合模型 ``` 2. **方法一:直接调用属性** ```python residual_variance = model.mse_resid # 直接获取残差方差 ``` 或使用: ```python residual_variance = model.scale # 与mse_resid等价[^2] ``` 3. **方法二:手动计算** ```python residuals = model.resid # 获取残差序列 ssr = np.sum(residuals**2) # 残差平方和SSR df_resid = model.df_resid # 自由度 = 样本数 - 自变量数 - 1 residual_variance = ssr / df_resid # 计算残差方差 ``` --- ### **关键公式** 残差方差定义为: $$ \hat{\sigma}^2 = \frac{SSR}{n - k - 1} $$ 其中: - $SSR = \sum_{i=1}^n (y_i - \hat{y}_i)^2$(残差平方和) - $n$为样本量 - $k$为自变量个数(不含常数项) --- ### **验证结果** 可通过以下代码验证两种方法的一致性: ```python print("方法一结果:", model.mse_resid) print("方法二结果:", residual_variance) ``` ---
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值