r - 如何使用多个变量和因子的所有可能组合运行 lm 模型

Question

这不是我的主题，所以如果我的问题问得不好或数据不完整，我很抱歉。我正在尝试运行 31 个具有单个响应变量 (VELOC) 的线性模型，并且作为预测变量有一个因子(TRAT，具有 2 个级别，A 和 B)和五个连续变量。

我有一个用于 gls 的循环，但只有连续预测变量，所以我认为它可以工作。但它没有，我相信这个问题一定是一件愚蠢的事情。

我不知道如何包含数据，但它看起来像:

   TRAT    VELOC        l        b        h        t        m

1     A  0.02490 -0.05283  0.06752  0.03435 -0.03343  0.10088

2     A  0.01196 -0.01126  0.10604 -0.01440 -0.08675  0.18547

3     A -0.06381  0.00804  0.06248 -0.04467 -0.04058 -0.04890

4     A  0.07440  0.04800  0.05250 -0.01867 -0.08034  0.08049

5     A  0.07695  0.06373 -0.00365 -0.07319 -0.02579  0.06989

6     B -0.03860 -0.01909  0.04960  0.09184 -0.06948  0.17950

7     B  0.00187 -0.02076 -0.05899 -0.12245  0.12391 -0.25616

8     B -0.07032 -0.02354 -0.05741  0.03189  0.05967 -0.06380

9     B -0.09047 -0.06176 -0.17759  0.15136  0.13997  0.09663

10    B -0.01787  0.01665 -0.08228 -0.02875  0.07486 -0.14252

现在，我使用的脚本是:

pred.vars = c("TRAT","l", "b", "h","t","m") #define predictor variables

m.mat = permutations(n = 2, r = 6, v = c(F, T), repeats.allowed = T)# I run        all possible combinations of pred.vars
models = apply(cbind(T, m.mat), 1, function(xrow) {paste(c("1", pred.vars)
[xrow], collapse = "+")})# fill the models 
models = paste("VELOC", models, sep = "~")#fill the left side
all.aic = rep(NA, length(models))# AIC of models
m.mat = cbind(1, m.mat)# Which predictors are estimated in the models beside
#the intercept

colnames(m.mat) = c("(Intercept)", pred.vars)

n.par = 2 + apply(m.mat,1, sum)# number of parameters estimated in the Models
coefs=m.mat# define an object to store the coefficients 

for (k in 1:length(models)) {
   res = try(lm(as.formula(models[k]), data = xdata))
   if (class(res) != "try-error") {
   all.aic[k] = -2 * logLik(res)[1] + 2 * n.par[k]
   xx = coefficients(res)
   coefs[k, match(names(xx), colnames(m.mat))] = xx
    }
}

我收到此错误: "Error in coefs[k, match(names(xx), colnames(m.mat))] = xx : NAs are not allowed in subscripted assignments"
在此先感谢您的帮助。对于如何正确发布问题的任何更正，我将不胜感激。
丽娜

Answer 1

我怀疑 dredge MuMIn 包中的函数会帮助你。您指定一个包含所有参数的“完整”模型，然后运行 ​​dredge(fullmodel)获取嵌套在完整模型中的所有组合。

然后，您应该能够从结果中获得系数和 AIC 值。

就像是:

require(MuMIn)
data(iris)

globalmodel <- lm(Sepal.Length ~ Petal.Length + Petal.Width + Species, data = iris)

combinations <- dredge(globalmodel)

print(combinations)

要获得所有模型的参数估计值(有点困惑)，您可以使用

coefTable(combinations)

或者要获取特定模型的系数，您可以使用疏浚对象中的行号对其进行索引，例如

coefTable(combinations)[1]

获取模型中第 1 行的系数。这也应该打印因子水平的系数。

见 MuMIn helpfile有关更多详细信息和提取信息的方法。

希望有帮助。