在R中为表的每一行计算线性趋势线

Question

是否可以在不使用循环的情况下对数据帧的每一行进行线性回归？趋势线的输出（截距+斜率）应添加到原始数据框中作为新列。

为了更清楚地表达我的意图，我准备了一个非常小的数据示例：

day1 <- c(1,3,1)
day2 <- c(2,2,1)
day3 <- c(3,1,5)
output.intercept <- c(0,4,-1.66667)
output.slope <- c(1,-1,2)
data <- data.frame(day1,day2,day3,output.intercept,output.slope)

输入变量是day1-3；假设这些是连续三天不同商店的销售额。我想要做的是为这三行计算一条线性趋势线，并将输出参数（查看output.intercept + output.slope）作为新列添加到原始表中。

解决方案应该在计算时间方面非常高效，因为真实数据框有成千上万的行。

最好，克里斯多夫

- user2635656

---

响应变量是什么？ - Sven Hohenstein

@SvenHohenstein 响应变量已经显示，协变量被暗示为 1:3（在这种情况下），在更一般的情况下是 seq_len(nrow(dat))。 - Gavin Simpson

Answer 1

design.mat <- cbind(1,1:3)
response.mat <- t(data[,1:3])

reg <- lm.fit(design.mat, response.mat)$coefficients
data <- cbind(data, t(reg))
#  day1 day2 day3 output.intercept output.slope        x1 x2
#1    1    2    3          0.00000            1  0.000000  1
#2    3    2    1          4.00000           -1  4.000000 -1
#3    1    1    5         -1.66667            2 -1.666667  2

然而，如果你有大量的数据，可能由于内存限制需要循环。如果是这种情况，我会使用长格式的data.table并使用该包的by语法进行循环。

- Roland

---

哇，完美运作。非常感谢！我稍后会用大数据集来尝试它。“design.mat”到底是什么？用于模拟x值吗？ - user2635656

如果你不知道设计矩阵是什么，那么你应该学习一本关于回归的教科书。 - Roland

再次感谢，解决方案即使在处理大数据时也非常完美。但是当“数据”包含缺失数据点（以NA形式）时会出现一个问题。（“lm.fit（design.mat，response.mat）中的NA / NaN / Inf in' y'”）有没有办法解决一些缺失数据点的问题？我已经尝试将na.exclude函数包含到lm.fit语句中，但在这种情况下它不起作用。 - user2635656

不幸的是，在一个特殊情况下，删除包含NA的行不是一个选项，因为在我的一个数据表中，几乎每一列都包含NA。难道没有另一种可能性只对可用数据进行回归，并“简单地”不考虑带有NA的值吗？否则，我将不得不事先操纵我的原始数据文件。 - user2635656

考虑一下你在请求什么。 - Roland

Answer 2

我和楼主遇到了同样的问题。这个解决方案适用于带有NAs的数据。之前的所有答案在这种情况下都会生成错误：

slp = function(x) {
  y = t(x)
  y = y[!is.na(y)]
  len = length(y):1
  b = cov(y,len)/var(len)
  return(b)}

reg_slp <- apply(data,1,slp)

只获取斜率，但截距可以很容易地添加。我怀疑这并不特别高效，但在我的情况下是有效的。

- user3055034

Answer 3

使用您的数据，

day1 <- c(1,3,1)
day2 <- c(2,2,1)
day3 <- c(3,1,5)
output.intercept <- c(0,4,-1.66667)
output.slope <- c(1,-1,2)
dat <- data.frame(day1,day2,day3)

我想您需要类似这样的内容：

fits <- lm.fit(cbind(1, seq_len(nrow(dat))), t(dat))
t(coef(fits))

这给出了什么

R> t(coef(fits))
         x1 x2
[1,]  0.000  1
[2,]  4.000 -1
[3,] -1.667  2

这些可以像这样添加到dat中。

dat <- cbind(dat, t(coef(fits)))
names(dat)[-(1:3)] <- c("Intercept","Slope")

R> dat
  day1 day2 day3 Intercept Slope
1    1    2    3     0.000     1
2    3    2    1     4.000    -1
3    1    1    5    -1.667     2

如果您有控制数据初始排列的能力，将数据以列为时间序列存储可能会更容易，而不是以行为基础。这样做可以避免在通过lm.fit()进行拟合时需要转置大矩阵。理想情况下，您希望数据最初按以下方式排列：

     [,1] [,2] [,3]
day1    1    3    1
day2    2    2    1
day3    3    1    5

即，将行作为时间点而不是像现在这样的单独系列。这是因为R期望数据的排列方式。请注意，在lm.fit()调用中，我们必须转置您的dat，这将涉及到一个大对象的复制。因此，如果您可以控制这些数据在进入R之前的排列/提供方式，那么对于大问题会有所帮助。

lm.fit()被用作底层精简代码，用于lm()，但我们避免了解析公式和创建模型矩阵的复杂性。如果您想要更高效的方法，您可能需要自己进行QR分解（代码在lm.fit()中执行此操作），因为lm.fit()会进行一些健全性检查，如果您确定数据不会导致奇异矩阵等问题，则可以省略这些检查。

- Gavin Simpson

---

非常感谢。我意识到在R中还有很多东西需要学习，甚至是基础知识。感谢您提供有关数据结构的提示。由于我之前在R中进行了一些数据准备工作，因此我对数据排列有控制权。我认为这样更有效率，因为我的真实数据文件包含600k行和仅100列。 - user2635656

一个注意事项：我认为语句“fits <- lm.fit(cbind(1, seq_len(nrow(dat))), t(dat))”应该调整为“fits <- lm.fit(cbind(1, seq_len(ncol(dat))), t(dat))”。或者我错了吗？它在示例中起作用是因为ncol(dat) = nrow(dat)。 - user2635656

Answer 4

或者像这样？

day1 <- c(1,3,1)
day2 <- c(2,2,1)
day3 <- c(3,1,5)
data <- data.frame(day1,day2,day3)
y<-1:3

reg<-apply(data,1,function(x) lm(as.numeric(x)~y))
data[,c("intercept","slope")]<-rbind(reg[[1]]$coef,reg[[2]]$coef,reg[[3]]$coef)

- DatamineR

---

这是正确的，但不是很高效。请注意，lm() 必须解析公式 nrow(dat) 次，如果只执行 3 次，则速度很快，如果执行 100K 次则速度慢。此外，这里遗漏了 lm() 的一个特性，即它接受矩阵响应。因此，您根本不需要在这里使用 apply() 或循环；您可以在单个 lm() 调用中拟合所有系列：lm(t(data[, 1:3]) ~ I(1:3))。但是，如果您想要高效的话，最好不要解析公式并生成模型框架和模型矩阵以及 lm() 给出的所有额外内容。对此，请使用 lm.fit() 来进行改进。 - Gavin Simpson