2

对于神经影像学应用程序，我试图通过 R 中的最小二乘法拟合许多线性模型（标准调用lm）。想象一下，我有一个设计矩阵 X。这个设计矩阵在所有模型中都是相同的。正在拟合的数据 (Y) 将发生变化，因此所有拟合参数（例如 beta、p 值、残差等）也会发生变化。

目前，我只是把它放在一个 for 循环中，所以它对lm. 似乎必须有更好的方法。

我相信计算成本最高的部分是矩阵求逆。看起来这可以通过 lm.fit 中的 Fortran 调用来处理。

如果我手动进行此回归，我会进行矩阵求逆，然后将其乘以各种数据集。事实上，当我有良好的设计矩阵（例如，所有连续赋值的协变量）时，我已经编写了一个函数来做到这一点。但是，我真的很喜欢所做的所有工作lm，比如适当地重新编码我的因素等，并且输出lm也非常好。

无论如何，我的蛋糕也可以吃吗？即，为了获得 lm 的友好性，但使用这种能力在计算上有效地拟合许多具有相同设计矩阵的模型？

4

3 回答 3

11

是的，有更好的方法。我们一直fastLm()在使用包 RcppArmadillo、RcppGSL 和 RcppEigen 中的 Armadillo、GSL 和 Eigen 的外部 C/C++ 代码编写示例替换函数。

到目前为止，花费最多的时间是建立模型矩阵和解析公式。您可以阅读 的源代码lm()，或者可能是我们的源代码fastLm()，并了解如何只进行一次解析。保留右侧，然后循环遍历不同的y向量。您使用哪种拟合功能不太重要。我喜欢fastLm()来自 RcppArmadillo，但是，嘿，我也写了 :)

于 2013-01-28T21:22:31.903 回答
8

从帮助页面lm：

如果“响应”是矩阵，则线性模型通过最小二乘法分别拟合到矩阵的每一列。

因此，一种简单的方法似乎是将所有不同的 y 向量组合成一个矩阵，并将其作为响应传递给lm. 例如：

(fit <- lm( cbind(Sepal.Width,Sepal.Length) ~ Petal.Width+Petal.Length+Species, data=iris))
summary(fit)
summary(fit)[2]
coef(summary(fit)[2])
coef(summary(fit))[2]
sapply( summary(fit), function(x) x$r.squared )
于 2013-01-28T22:13:16.717 回答
5

我不知道更好的使用方法lm；但您可能需要考虑使用 function lsfit。虽然更简单且花里胡哨的东西更少，但语法lsfit(X,y)允许y不仅是具有响应变量值的向量，而且是矩阵。然后 通过在同一设计矩阵上对它们进行回归来lsfit拟合所有列。相当快速和方便。yX

于 2013-01-28T21:07:03.517 回答