R MuMIn model.avg() - относительная важность не работает с вставленными формулами
Я пытаюсь использовать функцию MuMIn model.avg с модельными формулами, которые вставляются и используют индекс, а не прямой ввод, например:
m1<-gls(as.formula(paste(response,"~",paste(combns[,j], collapse="+"))), data=dat)
'combns' - это двумерный массив, созданный функцией combn(), содержащий комбинации переменных-предикторов. Это производит усредненные по модели коэффициенты и значения AICc, идентичные полученным, если функции gls содержат формулы непосредственно, например:
m1<-gls(median_Ta ~ day_of_season + hour_of_day + pct_grey_cover +
foliage_height_diversity + tree_shannon_diversity + median_patch_size, data=dat)
Однако относительная важность переменной не является вычислением, и я полагаю, что это связано с использованием цикла for или с использованием переменной для доступа к индексу списка, в котором хранятся модели, каким-то образом приводя к тому, что термины модели компонентов не будут "прочитайте" правильно (см. термины коды для моделей):
Component models:
df logLik AICc delta weight
1234567b 7 -233.08 481.43 0.00 0.59
1234567f 3 -237.97 482.21 0.78 0.40
1234567e 4 -241.32 491.08 9.65 0.00
1234567a 9 -241.15 502.39 20.96 0.00
1234567c 6 -248.37 509.68 28.25 0.00
1234567d 5 -250.22 511.11 29.68 0.00
Term codes:
day_of_season foliage_height_diversity hour_of_day
1 2 3
median_patch_size pct_grey_cover tree_shannon_diversity
4 5 6
urban_boundary_distance
7
Это приводит к тому, что относительная важность переменной задается как:
Relative variable importance:
day_of_season foliage_height_diversity hour_of_day
Importance: 1 1 1
N containing models: 6 6 6
median_patch_size pct_grey_cover tree_shannon_diversity
Importance: 1 1 1
N containing models: 6 6 6
urban_boundary_distance
Importance: 1
N containing models: 6
Принимая во внимание, что если я использую model.avg для тех же моделей с формулами, набранными индивидуально, я получаю следующий правильный вывод:
Component models:
df logLik AICc delta weight
23456 7 -233.08 481.43 0.00 0.59
1 3 -237.97 482.21 0.78 0.40
57 4 -241.32 491.08 9.65 0.00
1234567 9 -241.15 502.39 20.96 0.00
1467 6 -248.37 509.68 28.25 0.00
147 5 -250.22 511.11 29.68 0.00
Relative variable importance:
pct_grey_cover median_patch_size tree_shannon_diversity
Importance: 0.6 0.59 0.59
N containing models: 3 4 3
foliage_height_diversity hour_of_day day_of_season
Importance: 0.59 0.59 0.4
N containing models: 2 2 4
urban_boundary_distance
Importance: <0.01
N containing models: 4
Как я могу заставить model.avg правильно читать переменные предиктора в формулах? Я включил только шесть моделей в качестве примера, но я хочу сравнить полный набор из 128 моделей (и у меня есть другие переменные ответа с большим количеством переменных-предикторов), поэтому вывод их по отдельности невозможен.
Заранее спасибо.
Редактировать: воспроизводимый пример
Мне понадобилось время, чтобы сузить проблему. Первый пример, m.ave, показывает проблему в действии с циклом for. Во втором примере m.ave2 показывает, что он работает с типизированными индексами, а не с переменной. Очевидно, что это всего лишь небольшое подмножество переменных-предикторов.
require(nlme)
require(MuMIn)
dat<-data.frame(median_Ta=rnorm(100), day_of_season=runif(100), hour_of_day=runif(100), pct_grey_cover=rnorm(100),
foliage_height_diversity=rnorm(100), urban_boundary_distance=runif(100), tree_shannon_diversity=rnorm(100),
median_patch_size=rnorm(100))
f1<-"median_Ta ~ day_of_season + hour_of_day + pct_grey_cover + foliage_height_diversity +
urban_boundary_distance + tree_shannon_diversity + median_patch_size"
f1<-gsub("\\s", "", f1) # remove whitespace
f1split <- strsplit(f1, split="~") # split predictors and response
response <- f1split[[1]][1]
predictors <- strsplit(f1split[[1]][2], split="+", fixed=TRUE)[[1]]
modelslist<-list()
combns <- combn(predictors, 6)
for (j in 1:7) {
modelslist[[j]]<-gls(as.formula(paste(response,"~",paste(combns[,j], collapse="+"))), data=dat)
}
m.ave<-model.avg(modelslist[[2]], modelslist[[3]], modelslist[[4]],
modelslist[[5]], modelslist[[6]], modelslist[[7]], modelslist[[8]])
summary(m.ave)
#compare....
modelslist2<-list()
modelslist2[[1]]<-gls(as.formula(paste(response,"~",paste(combns[,1], collapse="+"))), data=dat)
modelslist2[[2]]<-gls(as.formula(paste(response,"~",paste(combns[,2], collapse="+"))), data=dat)
modelslist2[[3]]<-gls(as.formula(paste(response,"~",paste(combns[,3], collapse="+"))), data=dat)
modelslist2[[4]]<-gls(as.formula(paste(response,"~",paste(combns[,4], collapse="+"))), data=dat)
modelslist2[[5]]<-gls(as.formula(paste(response,"~",paste(combns[,5], collapse="+"))), data=dat)
modelslist2[[6]]<-gls(as.formula(paste(response,"~",paste(combns[,6], collapse="+"))), data=dat)
modelslist2[[7]]<-gls(as.formula(paste(response,"~",paste(combns[,7], collapse="+"))), data=dat)
m.ave2<-model.avg(modelslist2[[1]], modelslist2[[2]], modelslist2[[3]], modelslist2[[4]],
modelslist2[[5]], modelslist2[[6]], modelslist2[[7]])
summary(m.ave2)
1 ответ
Это ошибка в formula метод для gls (в упаковке nlme). Поскольку фактическая формула нигде не хранится в объекте, она оценивает "model" аргумент в вызове функции. В случае элементов modellist они все одинаковые, например:
modelslist[[1]]$call$model
modelslist[[7]]$call$model
оба возвращаются
> formula(paste(response, "~", paste(combns[, j], collapse = "+")))
который, когда evalиспользовать текущее (последнее) значение jтак что все formula(modellist[[N]]) возвращает последнюю формулу модели.
all.equal(formula(modelslist[[1]]), formula(modelslist[[7]]))
возвращается
> TRUE
Нужно сказать, что все это смущает model.avg который использует формулы для построения таблицы выбора модели (это запасной вариант, потому что gls не хватает terms также).
Изменить: возможные обходные пути
Намного проще получить то, что вы хотите:
model.avg(dredge(..., m.lim = c(6,6)))
или, если вы хотите делать прогнозы:
modellist <- lapply(dredge(..., m.lim = c(6,6), evaluate = FALSE), eval)
Но если вы хотите использовать произвольный набор моделей, замените $call$model элемент в каждом gls модельный объект с правильной формулой, например
combns <- combn(1:7, 6)
modellist <- vector("list", 7)
for (j in 1:7) {
f <- reformulate(predictors[combns[, j]], response = response)
fm <- gls(f, data = dat)
fm$call$model <- f # assign the actual formula
modellist[[j]] <- fm
}