R中的线性插值（lm），奇怪的行为答案

【问题标题】：Linear interpolation (lm) in R, weird behaviorR中的线性插值（lm），奇怪的行为
【发布时间】：2015-12-21 12:17:05
【问题描述】：

使用 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

知道为什么吗？好像哪里有溢出？

谢谢！

【问题讨论】：

标签： r lm linear-interpolation

【解决方案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'X https://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。

【讨论】：

是的，来自biglm 包的biglm 函数给出了相同的答案。
非常感谢您的快速回答，而且它确实适用于 tol 选项。问题是我对 y 的任何一组值都遇到了同样的问题，即test<-data.frame(dt=c(36996616, 36996620, 36996623, 36996626), value=c(4655145,2,-51434,215))。在这种情况下，x 和 y 的相关性较低，但存在相同的 NA 问题。
不客气。我中间提到的计算奇异误差来源于设计矩阵x，具体来说是solve(t(x) %*% x)。运行qr(t(x) %*% x)$rank，你会得到1，这意味着一个奇异矩阵。但是如果你运行qr(t(x) %*% x, tol=1e-31)$rank，你会得到2，这是一个非奇异矩阵。也就是说，我能够使用您更新的数据获得两个参数估计，所以我认为问题是高相关性和设计矩阵的组合，这仅取决于dt
@Yves 我做了一些额外的研究，并认为我有一些更好的信息。已相应更新我的答案

【解决方案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

【讨论】：