randomeffects.R

library(lme4)

bats <- 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")

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
glm_mixed <- lapply(formula_vec, function(f)   {
      form <- formula(paste(f, "(1 | IndvarType)", sep="+"))
      fit1 <- lmer(form, data = bats)
      fit1@call$fixed <- form
      return(fit1)
    }
)
names(glm_mixed) <- formula_vec




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)