A more 'manual' way to do it is this.
If you have all your models named properly and well organized, say Model1,
Model2, ..., Model 47, with a slot with the AIC (glm produces a list with one
component named 'aic') then something like this:
x <- matrix(0,1081,3)
x[,1:2] <- t(combn(47,2))
for(i in 1:n){x[,3] <-
abs(unlist(combn(eval(as.name(paste("Model",i,sep="")))$aic[1,])-unlist(combn(eval(as.name(paste("Model",i,sep="")))$aic,2)[2,]))}
will calculate all the 1081 AIC differences between paired comparisons of the
47 models. It may not be pretty but works for me.
An AIC difference equal or less than 2 is a tie, anything higher is evidence
for the model with the less AIC (Sakamoto et al., Akaike Information Criterion
Statistics, KTK Scientific Publishers, Tokyo).
HTH
Rubén H. Roa-Ureta, Ph. D.
AZTI Tecnalia, Txatxarramendi Ugartea z/g,
Sukarrieta, Bizkaia, SPAIN
-----Original Message-----
From: r-help-boun...@r-project.org on behalf of Milan Bouchet-Valat
Sent: Wed 1/25/2012 10:32 AM
To: Jhope
Cc: r-help@r-project.org
Subject: Re: [R] How do I compare 47 GLM models with 1 to 5 interactions and
unique combinations?
Le mardi 24 janvier 2012 à 20:41 -0800, Jhope a écrit :
> Hi R-listers,
>
> I have developed 47 GLM models with different combinations of interactions
> from 1 variable to 5 variables. I have manually made each model separately
> and put them into individual tables (organized by the number of variables)
> showing the AIC score. I want to compare all of these models.
>
> 1) What is the best way to compare various models with unique combinations
> and different number of variables?
See ?step or ?stepAIC (from package MASS) if you want an automated way
of doing this.
> 2) I am trying to develop the most simplest model ideally. Even though
> adding another variable would lower the AIC, how do I interpret it is worth
> it to include another variable in the model? How do I know when to stop?
This is a general statistical question, not specific to R. As a general
rule, if adding a variable lowers the AIC by a significant margin, then
it's worth including it. You should only stop when a variable increases
the AIC. But this is assuming you consider it a good indicator and you
know what you're doing. There's plenty of literature on this subject.
> Definitions of Variables:
> HTL - distance to high tide line (continuous)
> Veg - distance to vegetation
> Aeventexhumed - Event of exhumation
> Sector - number measurements along the beach
> Rayos - major sections of beach (grouped sectors)
> TotalEggs - nest egg density
>
> Example of how all models were created:
> Model2.glm <- glm(cbind(Shells, TotalEggs-Shells) ~ Aeventexhumed,
> data=data.to.analyze, family=binomial)
> Model7.glm <- glm(cbind(Shells, TotalEggs-Shells) ~ HTL:Veg, family =
> binomial, data.to.analyze)
> Model21.glm <- glm(cbind(Shells, TotalEggs-Shells) ~ HTL:Veg:TotalEggs,
> data.to.analyze, family = binomial)
> Model37.glm <- glm(cbind(Shells, TotalEggs-Shells) ~
> HTL:Veg:TotalEggs:Aeventexhumed, data.to.analyze, family=binomial)
To extract the AICs of all these models, it's easier to put them in a
list and get their AICs like this:
m <- list()
m$model2 <- glm(cbind(Shells, TotalEggs-Shells) ~ Aeventexhumed,
data=data.to.analyze, family=binomial)
m$model3 <- glm(cbind(Shells, TotalEggs-Shells) ~ HTL:Veg, family =
binomial, data.to.analyze)
sapply(m, extractAIC)
Cheers
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
