在 WLS 的情况下与 R 的 "Residual Standard Error"（以 lm 为单位）不一致

Question

我正在尝试使用 R 在 Excel 中重现加权最小二乘法 (WLS) 进行确认。我使用（琐碎但可重现的）以下数据集来执行双重检查：

x<-c(1,2,3,4,5,6)
y<-c(9,23,30,42,54,66)
w<-1/x

当我使用 lm 和权重参数计算 WLS 时，如下所示：

WLS<-lm(y~x, weights = w)
summary(WLS)

输出为：

> summary(WLS)

Call:
lm(formula = y ~ x, weights = w)

Weighted Residuals:
       1        2        3        4        5        6 
-0.50162  1.67280 -1.02017 -0.44984 -0.01447  0.34087 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.6311     1.2241  -1.333    0.254    
x            11.1327     0.4181  26.627 1.18e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.05 on 4 degrees of freedom
Multiple R-squared:  0.9944,    Adjusted R-squared:  0.993 
F-statistic:   709 on 1 and 4 DF,  p-value: 1.182e-05

我已阅读 here 可以使用以下行（为方便起见适应上述模型）手动计算从 R 计算的残差标准误差：

k=length(WLS$coefficients)-1 #Subtract one to ignore intercept
SSE=sum(WLS$residuals**2)
n=length(WLS$residuals)
sqrt(SSE/(n-(1+k))) #Residual Standard Error

这个计算和我在很多书上看到的公式是一致的(e.g. here)。但是，当运行这种手工计算时，返回的结果是1.618487（即不是1.05）。

我发现 here WLS 也可以通过将 OLS 应用于转换后的变量（矩阵符号中的模型：Y'=X'B+e'）使用以下转换来执行：Y=W^( 1/2)Y; X'=W^(1/2)X ; e'=W^(1/2)e。我使用以下代码在 R 中执行了它：

v<-w^(1/2)
x2<-x*v
y2<-y*v
WLS2<-lm(y2~0+v+x2)

即截距为零的模型，v 代表新的截距。这样做，我得到以下输出：

> summary(WLS2)

Call:
lm(formula = y2 ~ 0 + v + x2)

Residuals:
       1        2        3        4        5        6 
-0.50162  1.67280 -1.02017 -0.44984 -0.01447  0.34087 

Coefficients:
   Estimate Std. Error t value Pr(>|t|)    
v   -1.6311     1.2241  -1.333    0.254    
x2  11.1327     0.4181  26.627 1.18e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.05 on 4 degrees of freedom
Multiple R-squared:  0.9982,    Adjusted R-squared:  0.9972 
F-statistic:  1085 on 2 and 4 DF,  p-value: 3.388e-06

请注意，回归系数和残差标准误差相同，但 R²、F 统计量和残差不同。此外，当我使用该模型 (WLS2) 的残差计算残差标准误差时，我确实得到 1.049927。

我的问题：有人可以解释为什么尽管残差不同，但 R 返回的残差标准误差对于两个模型是相同的吗？对于第一个模型（没有数据转换），残差标准误差应该是 1.618487（手动计算）是否正确？ R 内部计算 WLS 的方式有问题吗？似乎 R 在计算残差标准误差之前忽略了对残差的反变换。

谢谢！

Answer 1

Is it correct that the Residual standard error should be 1.618487 (as calculated manually) for the first model (without data transformation) ?

否，因为您为模型拟合了权重，但随后在计算 sigma 时忘记了权重。

x <- c(1, 2, 3, 4, 5, 6)
y <- c(9, 23, 30, 42, 54, 66)
w <- 1 / x

WLS <- lm(y ~ x, weights = w)
summary(WLS)$sigma
#> [1] 1.049927

# You computed sigma for lm(y ~ x)
k <- length(WLS$coefficients) - 1 # Subtract one to ignore intercept
SSE <- sum(WLS$residuals**2)
n <- length(WLS$residuals)
sqrt(SSE / (n - (1 + k))) # Residual Standard Error, without weighting
#> [1] 1.618487

# But what you really wanted is to compute sigma for lm (y ~ x | weights = w)
SSE <- sum(w * (WLS$residuals)**2)
sqrt(SSE / (n - (1 + k))) # Residual Standard Error, with weighting
#> [1] 1.049927

^{由 reprex package (v2.0.1)}

于 2022 年 3 月 20 日创建

Why the Residual standard error returned by R are the same for the two models despite having different residuals ?

这是因为第二个模型明确包含权重：

WLS: y ~ 1 + x 权重 = w
WLS2: sqrt(w) * y ~ sqrt(w) + sqrt(w) * x 权重 = 1

在 WLS 的情况下与 R 的 "Residual Standard Error"（以 lm 为单位）不一致

Inconsistence with R's "Residual Standard Error" (in lm) in case of WLS

regression

r

transformation

weighted

glm