R 如何使用 lm() 函数计算回归系数

Question

我想复制 R 对以下数据的 regression equation 估计的计算：

set.seed(1)
Vec = rnorm(1000, 100, 3)
DF = data.frame(X1 = Vec[-1], X2 = Vec[-length(Vec)])

下面 R 报告系数估计值

coef(lm(X1~X2, DF))  ### slope =  -0.03871511

然后我手动估计斜率的回归估计

(sum(DF[,1]*DF[,2])*nrow(DF) - sum(DF[,1])*sum(DF[,2])) / (nrow(DF) * sum(DF[,1]^2) - (sum(DF[,1])^2)) ### -0.03871178

它们很接近，但仍然不完全匹配。

你能帮我理解我在这里错过了什么吗？

任何指针都会很有帮助。

Answer 1

这本身不是 Whosebug 问题，而是 sister site.

的统计问题

狭义的答案是您可以查看 R 源；它通常转移到 LAPACK 和 BLAS，但回归计算的一个关键部分专门用于处理（以统计而非数字方式）低等级案例。

无论如何，在这里，我相信你 'merely' 没有正确调整自由度，'almost but not quite' 当你使用 1000 个观察值时，自由度就会消失。下面是一个更简单的案例，以及计算系数 'by hand' 的 'simpler' 方法，它也具有匹配的优点：

> set.seed(1)
> Vec <- rnorm(5,100,3)
> DF <- data.frame(X1=Vec[-1], X2=Vec[-length(Vec)])
> coef(lm(X1 ~ X2, DF))[2]
       X2 
-0.322898 
> cov(DF$X1, DF$X2) / var(DF$X2)
[1] -0.322898
>

Answer 2

coef(lm(X1~X2, DF)) 
#     (Intercept)           X2 
#    103.83714016  -0.03871511

您可以按以下 OLS 矩阵形式应用 formula 系数。

X = cbind(1,DF[,2])
solve(t(X) %*% (X)) %*% t(X)%*% as.matrix(DF[,1])

给予，

#            [,1]
#[1,] 103.83714016
#[2,]  -0.03871511

与 lm() 输出相同。

数据：

set.seed(1)
Vec = rnorm(1000, 100, 3)
DF = data.frame(X1 = Vec[-1], X2 = Vec[-length(Vec)])

Answer 3

问题是X1和X2在lm中相对于长公式进行了切换

背景

lm(y ~ x) 的斜率公式如下，其中 x 和 y 的长度均为 n，x 是 x[i] 的缩写，y 是 y[i] 的缩写，求和结束我 = 1, 2, ..., n.

问题的根源

因此问题中的长公式，也显示在下面的 (1) 中，对应于 lm(X2 ~ X1, DF) 而不是 lm(X1 ~ X2, DF)。要么像下面的 (1) 那样更改 lm 模型中的公式，要么像下面的 (2) 那样通过将分母中每次出现的 DF[ 1] 替换为 DF[ 2] 来更改答案中的长公式。

# (1)

coef(lm(X2 ~ X1, DF))[[2]]
## [1] -0.03871178

(sum(DF[,1]*DF[,2])*nrow(DF) - sum(DF[,1])*sum(DF[,2])) / 
  (nrow(DF) * sum(DF[,1]^2) - (sum(DF[,1])^2))  # as in question
## [1] -0.03871178

# (2)

coef(lm(X1 ~ X2, DF))[[2]]  # as in question
## [1] -0.03871511

(sum(DF[,1]*DF[,2])*nrow(DF) - sum(DF[,1])*sum(DF[,2])) / 
  (nrow(DF) * sum(DF[,2]^2) - (sum(DF[,2])^2))
## [1] -0.03871511

R 如何使用 lm() 函数计算回归系数

How R calculates the Regression coefficients using lm() function

r

linear-regression

背景

问题的根源