Command that generates & evaluates model set from a mixed model?

[Matthew Savoca]

Is there a (relatively) straightforward R command to systematically evaluate a family of models using an information criterion (e.g., AIC or BIC) where the R command systematically removes/adds in different predictors in every possible combination? 

Example:

I want to run a GLMM, with response variable 'a'

I have various fixed and random predictor variables where:

variables b & c are fixed effects and 

variables d, e, & f are random effects. 

So, the full model in this set would be coded as follows:

a ~ b + c + (1|d) + (1|e) + (1|f)

I could do this by hand, but I assume theres a command that does this much faster/efficiently and would also ensure I'm not missing a model. I've looked into doing this with the dredge command in the MuMIn package, but it won't work and the error it's giving me leads me to believe that this command won't work with a full model that includes both fixed and random effects.

Thanks in advance for your help!

--

Matthew S. Savoca

Ph.D. Candidate; Graduate Group in Ecology

1147 Life Sciences Addition

University of California, Davis

office: 530-754-9500
cell: 718-490-9356

[Michael Hannon]

I don't know that this will help, but expand.grid() gives all
combinations of its inputs. See the appended for a simple, numerical
example.

-- Mike

> x <- c(3, 5, 7)
> y <- c(4, 6, 8)
> expand.grid(x, y)
Var1 Var2
1 3 4
2 5 4
3 7 4
4 3 6
5 5 6
6 7 6
7 3 8
8 5 8
9 7 8

> --
> Check out our R resources at
> http://www.noamross.net/davis-r-users-group.html
> ---
> You received this message because you are subscribed to the Google Groups
> "Davis R Users' Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an
> email to davis-rug+...@googlegroups.com.
> Visit this group at http://groups.google.com/group/davis-rug.
> For more options, visit https://groups.google.com/d/optout.

[Jaime Ashander]

Assuming you're not changing RE structure, you could do:

drop1(model) # for the simple model you specify.

For models with many terms, you may need to apply drop1, or similar commands, several times.

Both add1 and drop1 (and MASS equivalents) systematically evaluate models that are one single term away from the model passed, while stepAIC uses these to select a best model in a stepwise way. I'm not sure if it returns a table of all models

Commands to look at:

- stepAIC (from MASS package)
- add1 (in base) or addterm (from MASS)
- drop1 (in base) or dropterm (from MASS)

On Wed, Apr 29, 2015 at 4:19 PM, Matthew Savoca <msav...@gmail.com> wrote:

--

[Bonnie Dixon]

I am not aware of any tool that does that. (That functionality would be useful for me too in some circumstances.)

Here is what I once did to compare a bunch of mixed models programmatically, but it was a simpler situation than what you are describing. Anyway, perhaps this takes a step in the direction you are trying to go.

vars <- 
  names(df[ , setdiff(names(df), c("group.var", "outcome.var"))])

mod.func <- function(var.name) {
  formula <- as.formula(paste("outcome.var", "~", var.name))
  do.call(lme, 
          list(fixed = formula, 
               random = ~ 1| group.var, 
               data = df, 
               method = "ML"))
}

mods <- lapply(vars, mod.func)

mod.fit.stats <- rsquared.glmm(mods)
rownames(mod.fit.stats) <- vars
mod.fit.stats <- 
  mod.fit.stats[order(mod.fit.stats$AIC) , c("Marginal", "Conditional", "AIC")]
names(mod.fit.stats) <- c("Marginal_R2", "Conditional_R2", "AIC")
mod.fit.stats

Bonnie

[Noam Ross]

Rosie gave a presentation on doing this with glmulti a while ago: http://www.noamross.net/blog/2013/2/20/model-selection-drug.html, but only with fixed effects. It might work if you specify lmer rather than glm as the base model fitting function, but I’m not sure.

​

[Michael Koontz]

Wow, that is a really nicely written tutorial. I will definitely use that again.

If I remember correctly, log likelihood based model comparisons aren’t meaningful when both the fixed and random effects are changed. So you have to go through your whole process of choosing your random effects structure and then go through your process of choosing your fixed effects. (Or vice versa, I think). Is this right?

###

In case models with random effects still aren’t working with MuMIn, this code might give you some flexibility to define your random effects structure  (in the as.formula() call, I would think) and preferred modeling function. It always leaves all of the interactions in if there are more than 2 variables.

rm(list=ls())

data <- iris

response.idx <- 1 # index of response variable column

variables <- colnames(data)[-response.idx] # All of the column names except the response variable

response <- colnames(data)[response.idx]

variables

response

parameter.combos <- list()

# The combn function takes a vector (x), and creates all possible combinations of that vector using m elements from it. So a vector of length 3 with an m of 1 will have three possible combinations (each of the 3 elements by itself). A vector of length 2 with an m of 2 will have 3 possible combinations (1st and 2nd elements, 2nd and 3rd elements, 1st and 3rd elements), and so on...

# The combn() function returns a matrix the same mode as the x= argument (in this case character) with each column representing the different combinations of the x= argument with a different element of x= on each row. (So the number of rows is equal to m)

for (i in 1:length(variables))

{

  parameter.combos[[i]] <- combn(x=variables, m=i)  

}

parameter.combos

# For each set of parameter combinations, collapse each column into a single character vector element with a '*' as a separator

list.of.model.parameters <- lapply(parameter.combos, FUN=function(x) apply(x, MARGIN=2, FUN=function(z) paste(z, collapse="*")))

model.parameters <- unlist(list.of.model.parameters) # Unlist to make into a character vector

model.parameters

models <- list() # Create an empty list

# Go through the different parameter combinations, building a model each time, and saving the model result in the list

for (j in 1:length(model.parameters))

{

  formula <- as.formula(paste(response, " ~ ", model.parameters[j], sep="")) # Coerce the character string to be a formula.

  models[[j]] <- lm(formula, data=data) # Assign the jth list element of 'models' to be the lm() call.

}

models # A list of all of the models

AIC.table <- data.frame(model.parameters=model.parameters, AIC=NA, logLik=NA, df=NA)

for (k in 1:length(models))

{

  AIC.table[k, "AIC"] <- AIC(models[[k]])

  AIC.table[k, "logLik"] <- logLik(models[[k]])

  AIC.table[k, "df"] <- attr(logLik(models[[k]]), "df")

}

AIC.table <- AIC.table[order(AIC.table$AIC), ]

AIC.table

--

Michael Koontz — website

Graduate Group in Ecology

Plant and Environmental Sciences, Room 2211

University of California, Davis

Davis, CA 95616

[Rosemary Hartman]

Have you tried glmulti? I don't know that it does mixed models, but it automatically runs all the factor combinations and ranks them with AIC.

It's one of the few R things I know how to do! See here:

http://www.noamross.net/blog/2013/2/20/model-selection-drug.html

On Wed, Apr 29, 2015 at 4:19 PM, Matthew Savoca <msav...@gmail.com> wrote:

--

Check out our R resources at http://www.noamross.net/davis-r-users-group.html
---
You received this message because you are subscribed to the Google Groups "Davis R Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to davis-rug+...@googlegroups.com.
Visit this group at http://groups.google.com/group/davis-rug.
For more options, visit https://groups.google.com/d/optout.

--

Rosemary Hartman
roseh...@gmail.com
301-538-0136

[Rosemary Hartman]

Oh, shoot, Noam beat me to it! Anyway, good luck!

[Laurie Lawyer]

If you're using lmer to estimate your model, you might want to take a look at the step function in the lmerTest package.  It uses least squares not AIC, but the functionality is similar.

[Matthew Savoca]

This group is awesome! I have a lot of things to look into now, thanks a bunch all!

[Yee, Julie]

Hi, everyone,

I've been using a nifty package written by Mark Herzog (mhe...@usgs.gov) to do this very thing with mixed models.  He calls it 'ModelInference' and it's only available through him until he can get it documented more completely for public distribution.  He's graciously been willing to share it to those who inquire.  It works with model types in lme4 and glmmADMB, so you can include mixed effects and generalized distribution structures.  Not only does it permute all additive combinations of variables, but it can also construct all combinations of interaction effects, include quadratic terms, and exclude certain combinations from occurring.  For example, there is a "never.together" option where you can list predictors that should not occur together in the same model.  And a "required.together" option if you never want one to occur without the other.  You can also force certain predictors to occur in all models, and you can limit the complexity of the models by putting a cap on the number of terms or on the size of interaction effects.  The results of all models are collected in a model container, which can be summarized into a complete table of model rankings, sorted by AICc.  And there are functions that will calculate model-averaged coefficients and predictions.  Or, you can export the model container to a list of models and use the AICcmodavg package on it.  

Julie

____________________________________________________

Julie L. Yee, Statistician

U.S. Geological Survey

Western Ecological Research Center

www.werc.usgs.gov

Voice  530 - 669 - 5097

Fax     707 - 678 - 5039