在R中使用相同的设计矩阵拟合许多线性模型

Question

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

目前,我只是将它放在for循环中,因此它正在进行数十万次调用lm.似乎必须有一个更好的方法.

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

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

无论如何还有我的蛋糕和吃它？也就是说,为了获得lm的友好性,但是使用这种能力来计算有效地拟合具有相同设计矩阵的许多模型？

Answer 1

是的,还有更好的方法.我们一直在编写fastLm()基于使用来自Armadillo,GSL和Eigen的外部C/C++代码的示例替换函数,包括RcppArmadillo,RcppGSL和RcppEigen.

到目前为止,花费最多的时间来设置模型矩阵并使公式去除.你可以阅读的来源lm(),也许我们在fastLm(),看看如何做到这一点的解析只有一次.保持右侧,然后循环您的不同y向量.您使用哪种拟合功能更少.我喜欢fastLm()RcppArmadillo,但是嘿,我也写了:)

Answer 2

从帮助页面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 )

Answer 3

我不知道有更好的使用方式lm; 但你可能想考虑使用功能lsfit.虽然简单,具有较少花里胡哨,语法lsfit(X,y)允许y为不只是一个响应变量的值向量,也是一个矩阵.然后 通过在同一设计矩阵上回归它们来调整lsfit所有列的单个调用.相当快速和方便.yX