Как запустить модели 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"

Заранее спасибо за помощь. Буду признателен за любые исправления о том, как правильно размещать вопросы. Lina

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 для более подробной информации и способов извлечения информации.

Надеюсь, это поможет.

Answer 2

Иметь дело с:

Аргумент 'na.action' для global.model не задан, а для параметров (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)