R中所有可能的回归：将系数保存在矩阵中

Question

我正在为系统发育广义线性模型的所有可能模型运行代码。 我遇到的问题是为每个模型提取和保存 beta 系数。

我想将系数保存到矩阵中，其中列对应于特定变量，行对应于公式。问题出现是因为每个模型的变量都不同。因此，不能简单地将系数行绑定到矩阵。

下面的例子显示了我在问题上的解决方法：

y = rnorm(10)
inpdv = matrix(c(rnorm(10), runif(10), rpois(10, 1)), ncol = 3)
colnames(inpdv) = c("A", "B", "C")
data = cbind(y, inpdv)

model.mat = expand.grid(c(TRUE,FALSE), c(TRUE,FALSE), c(TRUE,FALSE))
names(model.mat) = colnames(inpdv)

formula = apply(model.mat, 1, function(x)
                       paste(colnames(model.mat)[x], collapse=" + "))
formula = paste("y", formula, sep = " ~ ")
formula[8] = paste(formula[8], 1, sep = "")

beta = matrix(NA, nrow = length(formula), ncol = 3)

for(i in 1:length(formula)){
   fit = lm(formula(formula), data)
   ## Here I want to extract the beta coeffecients here into a n * k matrix
   ## However, I cannot find a way to assign the value to the right cell in the matrix

}

所以我想每个系数都需要放入相应的单元格中，但我想不出一种快速有效的方法。

真正的分析将发生大约 30,000 次，因此任何有关效率的提示也将不胜感激。

编辑： 例如，y ~ a + c 模型的输出需要采用以下形式：

a NA b 

其中字母代表该模型的系数。 下一个模型可能是 y ~ b + c ，然后将其添加到底部。所以结果现在看起来像

a  NA b
NA b c

Answer 1

如何使用names 和%in% 对右列进行子集化。使用coef 提取系数值。像这样：

beta = matrix(NA, nrow = length(formula), ncol = 3)
colnames(beta) <- colnames(inpdv)

for(i in 1:length(formula)){
   fit = lm(formula(formula[i]), data)
    coefs <- coef(fit)
    beta[ i , colnames(beta) %in% names( coefs ) ] <- coefs[ names( coefs ) %in% colnames( beta ) ]
}
#              A          B         C
#[1,] -0.4229837 -0.0519900 0.3787666
#[2,]         NA  0.7015679 0.0555350
#[3,] -0.4165834         NA 0.3692974
#[4,]         NA         NA 0.1346726
#[5,] -0.2035173  0.7049951        NA
#[6,]         NA  0.7978726        NA
#[7,] -0.2229959         NA        NA
#[8,]         NA         NA        NA

我认为在这里使用for 循环是完全可以接受的。对于使用lapply 之类的初学者来说，当您运行越来越多的模拟时，有时会不断增加内存使用量。在lapply 循环完成之前，R 有时不会将早期模型中的对象标记为垃圾，因此有时会出现内存分配错误。使用for 循环，我发现R 将在必要时重用分配给循环上一次迭代的内存，因此如果您可以运行循环一次，则可以多次运行它。

不使用for 循环的另一个原因是速度，但我认为迭代时间与拟合模型的时间相比可以忽略不计，所以我会使用它。

Answer 2

for(i in 1:length(formula)){
    fit = lm(formula(formula), data)
     beta[i, 1:length(fit$coefficients)] <- fit$coefficients
}

更新

想法：以系数命名列，并按名称为列分配值。

这只是一个虚拟示例，但应该可以帮助您： 创建你的输出矩阵：

beta <- matrix(NA,  nrow=7, ncol=4)
colnames(beta) <- c("(Intercept)", 'A', 'B', 'C')

创建一些虚拟数据：

 A <- rnorm(10)
 B <- rpois(10, 1)
 C <- rnorm(10, 2)
 Y <- rnorm(10, -1)

现在你可以这样做了：

fit <- lm(Y ~ A + B + C)
beta[1, names(fit$coefficients)] <- fit$coefficients

fit <- lm(Y ~ A + B)
beta[2, names(fit$coefficients)] <- fit$coefficients

fit <- lm(Y ~ A + C)
beta[3, names(fit$coefficients)] <- fit$coefficients

fit <- lm(Y ~ B + C)
beta[4, names(fit$coefficients)] <- fit$coefficients

fit <- lm(Y ~ A)
beta[5, names(fit$coefficients)] <- fit$coefficients

fit <- lm(Y ~ B)
beta[6, names(fit$coefficients)] <- fit$coefficients

fit <- lm(Y ~ C)
beta[7, names(fit$coefficients)] <- fit$coefficients