如何使用多个变量和因子的所有可能组合运行 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

我怀疑 MuMIn 包中的 dredge 函数会对您有所帮助。您指定一个包含所有要包含的参数的“完整”模型，然后运行 ​​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 了解更多详细信息和提取信息的方法。

希望对您有所帮助。

Answer 2

处理：

'global.model' 的 'na.action' 参数未设置并且 options('na.action') 是 "na.omit"

require(MuMIn)
data(iris)
options(na.action = "na.fail") # change the default "na.omit" to prevent models
# from being fitted to different datasets in
# case of missing values.
globalmodel <- lm(Sepal.Length ~ Petal.Length + Petal.Width + Species, data = iris)
combinations <- dredge(globalmodel)
print(combinations)