Цикл для включения случайных эффектов в обобщенные линейные смешанные модели с использованием функции lmer()
Ниже приведено небольшое подмножество моего фрейма данных. Фактический фрейм данных имеет явное имя для каждой переменной; не только "DepVar1, а DepVar2 (2 переменные отклика)" или "IndVar (1-9)" (9 объясняющих переменных - 1 категориальная и 8 непрерывных переменных).
Я хотел бы адаптировать цикл, написанный Берганом, изменив функцию glm () на lmer (), найденную в пакете lme4, чтобы получить серию обобщенных линейных смешанных моделей (GLMM), содержащих ВСЕ возможные комбинации объясняющих переменных (Indvar 1-9) со случайными эффектами, указанными с использованием синтаксиса (1|IndVarType) для объяснения дисперсии переменной отклика (DepVar1 и DepVar2).
Example of glmm models:
DepVar1 ~ Indvar (1-9) + (1|IndVarType)
DepVar2 ~ Indvar (1-9) + (1|IndVarType)
После запуска цикла для создания всех моделей glmm моя цель состоит в том, чтобы отсортировать лучшие модели glmm по самым низким значениям AICc, используя функцию aictab () в пакете AICcmodavg для отображения соответствующей статистики: (1) Delta_AICc; (2) AICcWt; и (3) Cum.Wt.
Я пытался адаптировать код Берганса для включения случайных эффектов (1|IndVarType), но до сих пор мне это не удалось. Любые предложения, как это сделать? Я провел некоторые поиски и могу найти только примеры для циклов, содержащих функцию glm (). Большое спасибо заранее, если у кого-то есть решение.
Код
library(lme4)
ind_vars <- c("Indvar1",
"Indvar2",
"Indvar3",
"Indvar4",
"Indvar4",
"Indvar5",
"Indvar6",
"Indvar7",
"Indvar8",
"Indvar9",
"IndvarType")
dep_vars <- c("Depvar1", "DepVar2")
# create all combinations of ind_vars
ind_vars_comb <-
unlist(sapply( seq_len(length(ind_vars)),
function(i) {
apply( combn(ind_vars,i), 2, function(x) paste(x, collapse = "+"))
}))
# pair with dep_vars:
var_comb <- expand.grid(dep_vars, ind_vars_comb )
# formulas for all combinations
formula_vec <- sprintf("%s ~ %s", var_comb$Var1, var_comb$Var2)
# create models
# create models
glm_mixed <- lapply(formula_vec, function(f) {
fit1 <- lmer(f, (1|IndvarType), data = bats)
fit1$coefficients <- coef(summary(fit1))
return(fit1)
})
names(glm_mixed) <- formula_vec
##Error: No random effects terms specified in formula
# Model selection
# Installed AICcmodavg package for AICc values into R
Информация AICc
# R code from Mazerolle (2014)
library(AICcmodavg)
mydata.aov <- glm_mixed # list of models
mydata.model.names <- formula_vec # list of model names
# generates AICc values # sort models into order of AIC value
aictab(mydata.aov, mydata.model.names, second.ord = TRUE, sort = TRUE)
Структура данных
structure(list(Indvar1 = c(0, 5, 10, 19, 30, 33, 39, 44, 54,
63, 68, 72, 81, 87, 93, 100, 105, 110, 119, 127, 134, 141, 149,
155, 115, 120, 125, 0, 5, 9, 17, 22, 29, 35, 39, 44, 45, 50,
55, 63), IndvarType = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L), .Label = c("CONTROL", "LED", "Metal Halide", "SOX"), class = "factor"),
`IndvarCat ` = c(26.9, 25.16, 39, 29.81, 21.83, 20.22, 2.9,
2.1, 0.85, 0.62, 0.39, 0.26, 24.7, 21.99, 20.46, 26.32, 0,
0, 0.43, 0.02, 0.02, 0.03, 0.02, 0.03, 2.62, 0.43, 0.44,
25.16, 39, 29.81, 21.83, 20.22, 20.88, 0.63, 0.56, 0.56,
86.63, 87.97, 88.59, 0.31), Indvar2 = c(10.34, 12.56, 15.76,
10.35, 11.15, 14.6, 15.05, 12.54, 15.29, 19.5, 17.12, 17.62,
13.92, 12.7, 12.55, 17.86, 18.86, 18.23, 19.65, 19.59, 18.11,
19.04, 16.92, 18.39, 18.97, 18.96, 17.72, 7.65, 8.61, 8.98,
8.68, 12.25, 11.71, 16.19, 15.73, 16.02, 13.62, 14.89, 14.98,
17.14), Indvar3 = c(7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
5, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2, 11, 11, 11, 11, 11, 11,
11, 11, 11, 13, 13, 13, 13, 13, 8, 8), Indvar4 = structure(c(1L,
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 1L, 1L, 1L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L), .Label = c("Full Moon",
"Waning Gibbous", "Waxing Crescent", "Waxing Gibbous"), class = "factor"),
Indvar5 = c(32.2, 32.2, 32.2, 32.2, 32.2, 32.2, 32.2, 32.2,
32.2, 32.2, 32.2, 32.2, 32.2, 32.2, 32.2, 32.2, 32.2, 32.2,
32.2, 32.9, 32.9, 32.9, 32.9, 32.9, 41.4, 41.4, 41.4, 41.4,
41.4, 41.4, 41.4, 41.4, 41.4, 41.1, 41.1, 41.1, 41.1, 41.1,
42.2, 42.2), Indvar6 = c(2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,
2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3,
3, 3, 3, 3, 3, 3, 3, 3, 3, 3), Indvar7 = c(18, 18, 18, 18,
18, 18, 18, 18, 18, 18, 18, 18, 18, 14, 14, 14, 14, 14, 14,
14, 14, 13, 13, 13, 14.3, 14.3, 14.3, 14.3, 14.3, 14.3, 14.3,
14.3, 14.3, 15.5, 15.5, 15.5, 15.5, 15.5, 14.6, 14.6), Indvar8 = c(51,
51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 51, 69, 69, 69,
69, 69, 69, 77, 77, 77, 77, 77, 62, 62, 62, 62, 62, 62, 62,
62, 62, 57, 57, 57, 57, 57, 61, 61), Indvar9 = c(0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0), Depvar1 = c(3,
2, 5, 6, 15, 2, 10, 12, 17, 2, 0, 0, 15, 7, 17, 0, 1, 0,
14, 10, 12, 7, 4, 1, 5, 4, 2, 9, 7, 7, 9, 5, 4, 3, 0, 0,
12, 11, 9, 1), DepVar2 = c(0.444444444, 0, 0, 0.027777778,
0, 0, 0.25, 0, 0.08650519, 0, 0, 0, 0.111111111, 0, 0.124567474,
0, 0, 0, 0.25, 0.01, 0.111111111, 0.081632653, 0, 0, 0.04,
0.25, 0.25, 0.790123457, 0.510204082, 2.040816327, 1.777777778,
0, 2.25, 0.111111111, 0, 0, 0.027777778, 0.074380165, 0.012345679,
0)), .Names = c("Indvar1", "IndvarType", "IndvarCat ", "Indvar2",
"Indvar3", "Indvar4", "Indvar5", "Indvar6", "Indvar7", "Indvar8",
"Indvar9", "Depvar1", "DepVar2"), row.names = c(NA, 40L), class = "data.frame")