在R中的线性插值(lm)出现奇怪的行为

Question

使用R 3.2.2版本，我在运行一个简单的线性插值时发现了一种奇怪的行为。第一个数据框给出了正确的结果：

test<-data.frame(dt=c(36996616, 36996620, 36996623, 36996626), value=c(1,2,3,4))
lm(value~dt, test)$coefficients

  (Intercept)            dt 
-1.114966e+07  3.013699e-01 

通过增加dt变量，系数现在为NA：

test$dt<-test$dt+1
lm(value~dt, test)$coefficients

(Intercept)          dt 
        2.5          NA 

有任何想法吗？好像有地方溢出了？谢谢！

- Yves

Answer 1

编辑

我找到了一些更好的关于这个问题的信息。

如果预测变量完全相关，你会得到NA系数。但是，由于我们只有一个预测变量，这似乎是一个不寻常的情况。因此，在这种情况下，dt似乎与截距成线性相关。

我们可以使用alias找到线性相关的变量。请参见https://stats.stackexchange.com/questions/112442/what-are-aliased-coefficients

在第一个示例中：

test<-data.frame(dt=c(36996616, 36996620, 36996623, 36996626), value=c(1,2,3,4))
fit1 <- lm(value ~ dt, test)
alias(fit1)
Model :
value ~ dt

没有线性相关的术语。但在第二个例子中。

test$dt <- test$dt + 1
fit2 <- lm(value ~ dt, test)
alias(fit2)
Model :
value ~ dt

Complete :
   [,1]       
dt 147986489/4

这似乎显示了dt和intercept之间呈线性相关。

有关lm如何处理降秩模型的其他信息：https://stat.ethz.ch/pipermail/r-help/2002-February/018512.html。

lm不会反转X'Xhttps://stat.ethz.ch/pipermail/r-help/2008-January/152456.html，但我仍然认为以下内容有助于显示X'X的奇异性。

x <- matrix(c(rep(1, 4), test$dt), ncol=2)
y <- test$value

b <- solve(t(x) %*% x) %*% t(x) %*% y
Error in solve.default(t(x) %*% x) : 
system is computationally singular: reciprocal condition number = 7.35654e-30

lm.fit 的默认 tol 值为 1e-7。这是计算 qr 分解中的线性相关性所需的公差。

qr(t(x) %*% x)$rank
[1] 1

如果你减少这个值，你将得到一个关于dt的参数估计。

# decrease tol in qr
qr(t(x) %*% x, tol = 1e-31)$rank
[1] 2

# and in lm
lm(value~dt, test, tol=1e-31)$coefficients
  (Intercept)            dt 
-1.114966e+07  3.013699e-01 

请参考https://stats.stackexchange.com/questions/86001/simple-linear-regression-fit-manually-via-matrix-equations-does-not-match-lm-o，了解简单线性回归中的矩阵代数的详细信息。

- Whitebeard

---

是的，biglm包中的biglm函数会给出相同的答案。 - user3710546

非常感谢您的快速回复，确实使用 tol 选项效果很好。问题是我在任何一组 y 值中都遇到了相同的问题，即 test<-data.frame(dt=c(36996616, 36996620, 36996623, 36996626), value=c(4655145,2,-51434,215))。在这种情况下，x 和 y 的相关性较低，但存在相同的 NA 问题。 - Yves

不客气。我提到的计算奇异误差源于设计矩阵x，具体来说是solve(t(x) %*% x)。运行qr(t(x) %*% x)$rank，你会得到1，这意味着一个奇异矩阵。但是如果你运行qr(t(x) %*% x, tol=1e-31)$rank，你会得到2，这是一个非奇异矩阵。话虽如此，我能够用你更新的数据得到两个参数估计值，所以我认为问题是高相关性和仅取决于dt的设计矩阵的组合。 - Whitebeard

@Yves 我进行了一些额外的研究，认为我有更好的信息。已经相应地更新了我的答案。 - Whitebeard

Answer 2

biglm 函数来自于 biglm，似乎直接管理这个：

library(biglm)
test <- data.frame(dt=c(36996616, 36996620, 36996623, 36996626), 
                   value=c(1,2,3,4))
test$dt <- test$dt+1

coefficients(biglm(value ~ dt, test))
#   (Intercept)            dt 
# -1.114966e+07  3.013699e-01 

- user3710546