Below is a small subset of my data frame. The actual dataframe has an explicit name for each variable; not just "DepVar1, and DepVar2 (2 response variables)" or "IndVar (1-9)" (9 explanatory variables - 1 categorical and 8 continuous variables).
I'd like to adapt the loop written by Bergan by changing the function glm() to lmer() found in the lme4 package to produce a series of generalized linear mixed models (GLMM) containing ALL possible combinations of explanatory variables (Indvar 1-9) with random effects specified using a (1|IndVarType) syntax to explain variance in the response variable (DepVar1 and DepVar2).
Example of glmm models:
DepVar1 ~ Indvar (1-9) + (1|IndVarType)
DepVar2 ~ Indvar (1-9) + (1|IndVarType)
After running the loop to produce all glmm models, my aim to sort the best glmm models by the lowest AICc values using the function aictab() in the AICcmodavg package to display associated statistics: (1) Delta_AICc; (2) AICcWt; and (3) Cum.Wt.
I have been attempting to adapt Bergans code to incorporate random effects (1|IndVarType) but so far I have been unsuccessful. Any suggestions how to do that? I have done some searches and can only find examples for loops containing the glm() function. Many thanks in advance if anyone has a solution.
Code
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 information
# 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)
Data Structure
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")
