Generating AIC table with model and variable names

aic_table.R

require(survival)
require(kimisc) # has the nlist function to create a named list
require(AICcmodavg) # has the aictab function
require(dplyr)
require(ggplot2)
require(reshape2)

dat <- read.csv("data/obs_matched_pr.csv")

dat$strat_var <- factor(dat$strat_var)

dat_long <- melt(dat, id.vars=c("plot", "strat_var", "pres"))

# Code for facet (small multiples) plot to see all variable at once
ggplot(dat_long, aes(x=value, y = pres)) + 
  facet_wrap(~ variable) + # add scales="free_x" if you want to use the whole x axis
  geom_point() + 
  geom_smooth(method = "glm", family="binomial") # the 


# Generate some models
mod1 <- clogit(pres ~ sand + grindellia + strata(strat_var), data = dat)
mod2 <- clogit(pres ~ sand + strata(strat_var), data = dat)
mod3 <- clogit(pres ~ grindellia + strata(strat_var), data = dat)

# Put the models all together in a named list using nlist function from kimisc package
model_list <- nlist(mod1, mod2, mod3)

# Compare models with AIC table
aic_table <- aictab(model_list)

# Add a column to the table with variable names
model_vars <- sapply(model_list, function(x) gsub(":", "*", paste(names(x$means), collapse = " + ")))
aic_table <- cbind(Modnames = names(model_vars), model_vars) %>%
  merge(aic_table, by="Modnames") %>%
  arrange(AICc)

# http://r.789695.n4.nabble.com/Re-prediction-based-on-conditional-logistic-regression-clogit-td4692244.html

# Example plotting for model 1
# Get odds ratios and confidence intervals and put them together in a dataframe
OR <- exp(coef(mod1))
CI <- exp(confint(mod1, level = .95))
ggdf <- data.frame(var = names(OR), OR, LCI = CI[,1], UCI = CI[,2], row.names = NULL)

ggplot(ggdf, aes(x = var, y = OR)) + 
  geom_point(size = 5) + 
  geom_errorbar(aes(ymin=LCI, ymax=UCI), size=1, width=0) + 
  labs(title="Odds ratios for factors in model 1", 
       x = "Factor", y = "Odds Ratio (95% CI)") + 
  coord_flip()

# compute model averages

vars_to_avg <- c("sand", "grindellia")

avgs <- sapply(vars_to_avg, function(x) modavg(model_list, x), simplify = FALSE)

avg_parms <- lapply(avgs, function(x) {
  var = x$Parameter
  beta = x$Mod.avg.beta
  uncond_se = x$Uncond.SE
  LCI = x$Lower.CL
  UCI = x$Upper.CL
  data.frame(var = var, beta=beta, uncond_se=uncond_se, OR = exp(beta), 
             LCI=exp(LCI), UCI=exp(UCI))
}) %>% rbind_all()

# Plot model averaged results
ggplot(avg_parms, aes(x = var, y = OR)) + 
  geom_point(size = 5) + 
  geom_errorbar(aes(ymin=LCI, ymax=UCI), size=1, width=0) + 
  labs(title="Odds ratios for averaged models 1:3", 
       x = "Factor", y = "Odds Ratio (95% CI)") + 
  coord_flip()