加权最小二乘（Weighted Least Squares，WLS）回归及R语言计算

zrunrun

已于 2022-07-03 22:18:46 修改

阅读量8k

点赞数 3

CC 4.0 BY-SA版权

文章标签： r语言回归

于 2022-06-24 21:31:42 首次发布

本文链接：https://blog.csdn.net/weixin_61593847/article/details/125447278

该博客探讨了在数据存在异方差性的情况下，如何使用加权最小二乘（WLS）回归来改进模型拟合。作者首先展示了数据的分布，并通过图示证明了数据符合幂函数模型。接着，通过对数据取对数进行线性化，使用OLS回归得到初步模型。然后，注意到残差与流速的关系，即流速慢的观测值残差较大，作者提出了WLS回归，通过定义权重来处理这种异方差性。尽管WLS的R方值更高，但参数与OLS接近，且拟合曲线几乎重合。最后，通过计算残差平方和，验证了WLS计算R方的正确公式。博客强调，相关系数不是评估权重选择优劣的合适指标，而WLS有助于使拟合更贴近残差小的观测值。

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

与OLS回归相比，加权最小二乘（Weighted Least Squares，WLS）回归是当残差中方差不变的最小二乘假设被违背（异方差性）时可以考虑的一种方法，即different observations have different residuals。

> #部分示例数据
> head(dat_mod)
      time pace diameter
1 168.5000   40     5.28
2 179.4762   40     5.28
3 155.8077   50     5.28
4 166.5278   50     5.28
5 148.2500   50     5.28
6 143.5714   60     5.28
#做图查看数据分布
dat_mod$diameter=as.factor(dat_mod$diameter)
ggplot(data=dat_mod,aes(x=pace,y=time,color=diameter,shape=diameter))+
  geom_point()+
  scale_y_continuous(limits = c(60,190),n.breaks = 14)+
  scale_x_continuous(breaks = c(30,40,50,60,90,120,150,180,210,240,270,300,330,360,390))+
  scale_color_manual(values=c('#228B22','#FFFF00','#FF0000'))

在这里插入图片描述
可以看出数据比较符合幂函数模型，所以思路是对y和x都取对数，将非线性模型线性化，然后再用最小二乘法拟合。

#幂函数模型
y3=log(dat_mod$time)
mod3=lm(y3~log(dat_mod$pace))
summary(mod3)
> summary(mod3)

Call:
lm(formula = y3 ~ log(dat_mod$pace))

Residuals:
      Min        1Q    Median        3Q       Max 
-0.186191 -0.024459  0.007359  0.031775  0.147288 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        6.37989    0.05151  123.86   <2e-16 ***
log(dat_mod$pace) -0.36094    0.01032  -34.97   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.06502 on 63 degrees of freedom
Multiple R-squared:  0.951,	Adjusted R-squared:  0.9502 
F-statistic:  1223 on 1 and 63 DF,  p-value: < 2.2e-16

pace3=seq(20,400,length.out=1000)
time_pred3=(pace3^(-0.36094))*exp(6.37989)
#或者mod3$fitted.values
dat_pred3=data.frame(pace3,time_pred3)


ggplot()+
  geom_point(data=dat_mod,aes(x=pace,y=time,color=diameter,shape=diameter),size=3)+
  scale_y_continuous(limits = c(60,200),n.breaks = 15)+
  scale_x_continuous(breaks = c(20,30,40,50,60,90,120,150,180,210,240,270,300,330,360,390,400))+
  scale_color_manual(values=c('#228B22','#1E90FF','#FF0000'))+
  theme(panel.grid=element_blank(),panel.background=element_rect(fill='white', color='black',size=0.5),axis.text.x=element_text(size=18),axis.title.x=element_text(size=22,face='bold'),axis.text.y=element_text(size=18),axis.title.y=element_text(size=22,face='bold'),legend.title = element_text(size=22,face='bold'),legend.text = element_text(size=18))+
  labs(x='pace(ug/min)',y='time(s)',color='diameter(cm)',shape='diameter(cm)')+
  geom_line(data=dat_pred3,aes(x=pace3,y=time_pred3),color='#000000',size=1)

在这里插入图片描述
可以看出流速慢的时候残差大，流速越快残差越小，符合WLS使用条件异方差性
OLS是minimize sum(residuals^2),而WLS是minimize sum(w*residuals^2),即将权数与残差平方相乘后再求和，所以要先定义权重。

#定义权重
wt=1/(mod3$residuals)^2

#加权最小二乘（WLS）回归，通过参数 weight 指定加权的变量
fit.wls=lm(y3~log(dat_mod$pace),weights = wt)
summary(fit.wls)
###################################################

Call:
  lm(formula = y3 ~ log(dat_mod$pace), weights = wt)

Weighted Residuals:
  Min      1Q  Median      3Q     Max 
-0.9949 -0.9613  1.0009  1.0294  1.4474 

Coefficients:
  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        6.372919   0.006742   945.3   <2e-16 ***
  log(dat_mod$pace) -0.359768   0.001201  -299.7   <2e-16 ***
  ---
  Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.9959 on 63 degrees of freedom
Multiple R-squared:  0.9993,	Adjusted R-squared:  0.9993 
F-statistic: 8.98e+04 on 1 and 63 DF,  p-value: < 2.2e-16

###################################################

pace4=seq(20,400,length.out=1000)
time_pred4=(pace4^(-0.359768))*exp(6.372919)
#或者predict()
dat_pred4=data.frame(pace4,time_pred4)


ggplot()+
  geom_point(data=dat_mod,aes(x=pace,y=time,color=diameter,shape=diameter),size=3)+
  scale_y_continuous(limits = c(60,200),n.breaks = 15)+
  scale_x_continuous(breaks = c(20,30,40,50,60,90,120,150,180,210,240,270,300,330,360,390,400))+
  scale_color_manual(values=c('#228B22','#1E90FF','#FF0000'))+
  theme(panel.grid=element_blank(),panel.background=element_rect(fill='white', color='black',size=0.5),axis.text.x=element_text(size=18),axis.title.x=element_text(size=22,face='bold'),axis.text.y=element_text(size=18),axis.title.y=element_text(size=22,face='bold'),legend.title = element_text(size=22,face='bold'),legend.text = element_text(size=18))+
  labs(x='pace(ug/min)',y='time(s)',color='diameter(cm)',shape='diameter(cm)')+
  geom_line(data=dat_pred4,aes(x=pace4,y=time_pred4),color='#8B008B',size=1)+
  geom_line(data=dat_pred3,aes(x=pace3,y=time_pred3),color='#000000',size=1)

在这里插入图片描述
这是使用了WLS后得到的拟合曲线，虽然r方=0.9993较之前r方=0.9502增加了很多，但是参数和OLS得到的参数基本一样，得到的拟合曲线基本重合，为什么会这样呢？
我们先比较下两种算法得到的残差平方和：

> sum(mod3$residuals^2)
[1] 0.2663267
> sum(fit.wls$residuals^2)
[1] 0.2664723

可以看到由WLS计算得到的残差平方和要大于由OLS计算得到的残差平方和，如果按照皮尔森相关系数的计算公式1-（RSS/TSS）那么WLS计算得到的R方肯定会更小，所以我猜测WLS计算R方的公式应该是
1-（sum(w×residuals∧2) / sum(w×(y-mean(y))∧2)），下面来验证一下

> tss=(wt*(y3-mean(y3))^2)%>%sum()
> rss=(wt*fit.wls$residuals^2)%>%sum()
> 1-rss/tss
[1] 0.9999152

总结

在加权回归分析中, 相关系数不能作为评判权重选择是否最优的指标,因为计算得到的是加权相关系数而不是皮尔森相关系数，但是采用加权最小二乘法可以使拟合曲线更加靠近残差小的数据点。