1

我正在为系统发育广义线性模型的所有可能模型运行代码。我遇到的问题是为每个模型提取和保存 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
4

2 回答 2

4

如何使用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，当您运行越来越多的模拟时，有时会不断增加内存使用量。在循环结束之前，R 有时不会将早期模型中的对象标记为垃圾，lapply因此有时会出现内存分配错误。使用for循环，我发现 R 将在必要时重用分配给循环上一次迭代的内存，因此如果您可以运行循环一次，则可以多次运行它。

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

于 2013-09-11T11:54:54.650 回答
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
于 2013-09-11T09:31:37.113 回答