6

考虑以下 R 代码（我认为它最终会调用一些 Fortran）：

X <- 1:1000
Y <- rep(1,1000)
summary(lm(Y~X))

为什么摘要返回值？由于 Y 没有变化，这个模型不应该不适合吗？更重要的是，为什么模型 R^2 ~= .5？

编辑

我跟踪了从 lm 到 lm.fit 的代码，可以看到这个调用：

z <- .Fortran("dqrls", qr = x, n = n, p = p, y = y, ny = ny,
   tol = as.double(tol), coefficients = mat.or.vec(p, ny), residuals = y,
   effects = y, rank = integer(1L), pivot = 1L:p, qraux = double(p),
   work = double(2 * p), PACKAGE = "base")

这就是真正适合的地方。查看http://svn.r-project.org/R/trunk/src/appl/dqrls.f）并没有帮助我理解发生了什么，因为我不知道 fortran。

4

3 回答 3

5

从统计学上讲，我们应该预期什么（我想说“预期”，但这是一个非常具体的术语 ;-)）？系数应为 (0,1)，而不是“无法拟合”。假设 (X,Y) 的协方差与 X 的方差成正比，而不是相反。由于 X 具有非零方差，因此没有问题。由于协方差为 0，X 的估计系数应为 0。因此，在机器容差范围内，这就是您得到的答案。

这里没有统计异常。可能存在统计上的误解。还有机器公差的问题，但考虑到预测变量和响应值的规模，1E-19 数量级的系数可以忽略不计。

更新 1：可以在此 Wikipedia 页面上找到对简单线性回归的快速回顾。需要注意的关键Var(x)是在分母中，Cov(x,y)在分子中。在这种情况下，分子为 0，分母为非零，因此没有理由期望 aNaN或NA。但是，有人可能会问为什么 a 的结果系数不是x，0这与 QR 分解的数值精度问题有关。

于 2012-02-12T01:44:08.890 回答
2

我相信这仅仅是因为 QR 分解是用浮点运算实现的。

该singular.ok参数实际上是指设计矩阵（即仅X）。尝试

lm.fit(cbind(X, X), Y)

对比

lm.fit(cbind(X, X), Y, singular.ok=F)
于 2012-02-12T00:26:18.663 回答
2

我同意问题可能出在浮点数上。但我不认为是奇点。

如果您检查使用solve(t(x1)%*%x1)%*%(t(x1)%*%Y)而不是QR，(t(x1)%*%x1)则不是单数

使用x1 = cbind(rep(1,1000,X)因为lm(Y~X)包括拦截。

于 2012-02-12T01:33:44.263 回答