fishlakes.csv

lakes.df2.csv

model selection DRUG .R

# Model selection and multi-model inference
# A not-particularly good introduction by Rosemary Hartman
# (yes, I am just doing this for the flying monkey)

# So, let's say you want to find out where things are
# and why they are there.

# But there are a lot of reasons someone might be somewhere

# Let's say these things are fishermen,
# And we aren't sure what the reasons are yet.

# Start by making some good guesses. Maybe they all go to the 
# big lakes. Maybe they all go to the lakes that CDFW stocks with fish.
# Maybe they go to the lakes with lots of frogs (not likely, but hey, why not try?)
# Once you have all your guesses, you need to figure out which one is right.
# There are two basic methods to use (or two that I have been exploring)

# Method #1: The kitchen sink

# load package "glmulti"
library(glmulti)

# create a model that has all the predictor variables you would like to test
global.model <-  glm(fish~ # number of fishermen
                       stocked+ # CDFW fish stocking y/n
                       cows+ # cattle grazing in area y/n
                       herp.rich+ # amphibian species richness
                       fish.rich+ # fish species richness
                       visitors+ # total number of visitors
                       logsize, # log of the area
                     data=lakes.df2) # data frame we got this from

# Now we will use the "glmulti" function to find the best model;
# this goes through every possible model and finds the best one.
  
fish.model <- glmulti(global.model, # use the model with built as a starting point
                      level = 1,  #  just look at main effects
                      crit="aicc") # use AICc because it works better than AIC for small sample sizes

summary(fish.model)

# That showed us the best model, now lets look at some of the others
weightable(fish.model)

# So this is the best model
  f <- glm(fish~1+cows+herp.rich+visitors+logsize, data=lakes.df2)
summary(f)

# I hate trying to interpret models based on tables of coeficients.
# Let's look at some graphs
library(visreg)
visreg(f)

# But according to theory, models with AIC within two points of each other are basically
# equal. So what about the other models? Should we totally throw them out? Looking
# at the table of aicc weights, there is a pretty big jump between model 8 and model 9.
# So lets try averaging the top 8 models if we want to use this to make predictions,
# evaluate variable importance, etc.

library(MuMIn)

# run the top 8 models
f2 <- glm(fish ~ 1 + cows + herp.rich + fish.rich + visitors + logsize, data=lakes.df2)
f3 <- glm(fish ~ 1 + stocked + cows + herp.rich + visitors + logsize, data=lakes.df2) 
f4 <- glm(fish ~ 1 + cows + visitors + logsize, data=lakes.df2)  
f5 <- glm(fish ~ 1 + cows + fish.rich + visitors + logsize, data=lakes.df2)  
f6 <- glm(fish ~ 1 + stocked + cows + visitors + logsize, data=lakes.df2)  
f7 <- glm(fish ~ 1 + stocked + cows + herp.rich + fish.rich + visitors + logsize, data=lakes.df2)  
f8 <- glm(fish ~ 1 + stocked + cows + fish.rich + visitors + logsize, data=lakes.df2)  

# average the models together
f.ave <- model.avg(f, f2, f3, f4, f5, f6, f7, f8)
summary(f.ave)

# Now we know the relative variable importance, all of the avereged coefficents,
# etc, and can use the "predict" funciton to predict how many fishermen will
# be at the next lake. Graphing these are more difficult, so we'll skip this for now. 
# Ask me later if you really want to know.

# Model selection method #2: Use your brain
# We often can discard (or choose) some models a priori based on our knowlege of
# the system. It's usually better to do it this way if you have several hundered
# possible combination of variables, or want to put in some interaction terms.
# I used this method for my frog data.

# load package bbmle
library(bbmle)

# Decide on which set of models you want to use. This is hard. A statistician
# who knows a lot more than I do told me so. I spent a long time playing around
# with different transformations, predictor variables, and combinations
# before I came up with a set of hypotheses (models) that I was happy with.

m1 <- glm(treatment ~ 1+ logveg + bank.slope, family=binomial("logit"), data= fishlakes)
m2 <- glm(treatment ~ 1+ logveg + silt.total +bank.slope,  family=binomial("logit"), data= fishlakes) 
m3 <- glm(treatment ~ 1+ raca.basin + logveg + bank.slope,  family=binomial("logit"), data= fishlakes) 
m4 <- glm(treatment ~ 1+ BUBO.breeding + logveg + bank.slope,   family=binomial("logit"), data= fishlakes) 
m5 <- glm(treatment ~ 1+ BUBO.breeding + logveg + bank.slope +  raca.basin, family=binomial("logit"), data= fishlakes) 
m6 <- glm(treatment ~ 1+ raca.basin*BUBO.breeding, family=binomial("logit"), data= fishlakes)
m7 <- glm(treatment ~ 1+ herbacious +  raca.basin + lakes.basin, family=binomial("logit"), data= fishlakes) 
m8 <- glm(treatment ~ 1+ BUBO.breeding, family=binomial("logit"), data= fishlakes) 
m9 <- glm(treatment ~ BUBO.breeding + herbacious + bank.slope +  logveg, family=binomial("logit"), data= fishlakes)
m10 <- glm(treatment ~ BUBO.breeding*raca.basin + herbacious+ lakes.basin, family=binomial("logit"), data= fishlakes) 
m11 <- glm(treatment ~ 1 + bank.slope + logveg + BUBO.breeding:raca.basin, family = binomial("logit"), data = fishlakes)

# Now let's rank them via AICc
AICctab(m1,m2,m3,m4,m5,m6,m7,m8,m9,m10,m11, base=T, weights=T, nobs=length(fishlakes))

# Model m11 comes out ahead by quite a bit, but we'll average the top two models, just to show
# you how its done.

m.ave <- model.avg(m4, m11)
summary(m.ave)

#OK, there is a predictive model, but how good is it?
# Let's try cross-validation first. If this was a single model, we could try using
# the cv.glm function from the "boot" library like this:

library("boot")
cost <- function(r, pi = 0) mean(abs(r-pi) > 0.5) ## cost function necessary for binomial data
m11.cv <- cv.glm(data=fishlakes, m11, cost, K=42) # use leave-one-out cross validation (can use K-fold cross validation for larger data sets)

# Now lets see what our error rate was:
m11.cv$delta
# That's not too bad.

# IF we want to check the error rate of an averaged model, you need to get more creative.
# I've written some code to do this for averaged models that only have two component models.
# It shouldn't be too hard to adapt this for more models.

# function for leave one out cross validation of averaged models

Cross <- function(model1, model2, data, cost) {
  library(MuMIn)
  nobs <- nrow(data)
  model.ave <- model.avg(model1, model2)
  values <- matrix(NA, nrow=nobs, ncol=5)
  values[,1] <- data$treatment
  values[,2] <- predict(model.ave, type="response")
  CV=0
  for (i in 1:nobs) {
    data2 <- data[-i,] # leave out one observasion
    model12 <- glm(model1$formula, family=model1$family, data=data2)
    model22 <- glm(model2$formula, family=model2$family, data=data2)
    model.ave2 <- model.avg(model12, model22)
    values[i,3] <- predict(model.ave2, newdata=data[i,], type="response")
    values[i,4] <- round(values[i,3])
    if (values[i,4]==values[i,1]) (values[i,5]=1) else values[i,5]=0
    
  }
  cv = mean(abs(values[,3]-values[,1])>0.5)
  return(cv)
}


# Use the function on the component models
Cross(m11, m4, fishlakes, cost)
# That looks reasonable

# Another method to test model accurace is Area Under the Reciever Operater Curve (AUC)
# This is baisically a plot of true presences versus false presences in a 
# presence-absense model.

# Load the library "pROC"
library(pROC)

# Make your reciever-operater curve 
m.roc <-roc(fishlakes$treatment, predict(m.ave, backtransform=TRUE))
plot(m.roc)
# Looks like a pretty good fit. Not too bad for the small size of the data set.

# And that's all I got. Hopefully you will find it helpful. Let me know if there is anything
# else I have forgotton or done wrong. ~ Rosemary rosehartman@gmail.com