setwd('C:\\Users\\v\\Documents\\My Documents full backup\\ODU\\Modern topics seminar\\Day 04') library(AICcmodavg) data(dry.frog) ##setup a subset of models of Table 1 in Mazerolle 2006 (it's not really the same analysis, just an example) Cand.models <- list( ) Cand.models[[1]] <- lm(log_Mass_lost ~ Shade + Substrate + cent_Initial_mass + Initial_mass2, data = dry.frog) Cand.models[[2]] <- lm(log_Mass_lost ~ Shade + Substrate + cent_Initial_mass + Initial_mass2 + Shade:Substrate, data = dry.frog) Cand.models[[3]] <- lm(log_Mass_lost ~ cent_Initial_mass + Initial_mass2, data = dry.frog) Cand.models[[4]] <- lm(log_Mass_lost ~ Shade + cent_Initial_mass + Initial_mass2, data = dry.frog) Cand.models[[5]] <- lm(log_Mass_lost ~ Substrate + cent_Initial_mass + Initial_mass2, data = dry.frog) Cand.models[[6]] <- lm(log_Mass_lost ~ 1, data = dry.frog) ##create a vector of names to trace back models in set Modnames <- paste("mod", 1:length(Cand.models), sep = " ") ##generate AICc table aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE) ##Another way to generate a table, with rounding to 4 digits after decimal point and give log-likelihood print(aictab(cand.set = Cand.models, modnames = Modnames, sort = TRUE), digits = 4, LL = TRUE) ##compute model-averaged value and unconditional SE of predicted log of ##mass lost for frogs of average mass in shade for each substrate type ##first create data set to use for predictions new.dat <- data.frame(Shade = c(1, 1, 1), cent_Initial_mass = c(0, 0, 0), Initial_mass2 = c(0, 0, 0), Substrate = c("SOIL", "SPHAGNUM", "PEAT")) ##compute unconditional SE's modavgPred(cand.set = Cand.models, modnames = Modnames, newdata = new.dat, type = "response", uncond.se = "revised") ##round to 4 digits after decimal point print(modavgPred(cand.set = Cand.models, modnames = Modnames, newdata = new.dat, type = "response", uncond.se = "revised"), digits = 4) ############################################################################################# #model averaged estimates ##beware of variables that have similar names data1=read.csv("example resp data.csv") ##set up models in list Cand <- list( ) Cand[[1]] <- lm(resp ~ size + agecorr, data = data1) Cand[[2]] <- lm(resp ~ size + mass + agecorr, data = data1) Cand[[3]] <- lm(resp ~ age + mass, data = data1) Cand[[4]] <- lm(resp ~ age + mass + mass2, data = data1) Cand[[5]] <- lm(resp ~ mass + mass2 + size, data = data1) Cand[[6]] <- lm(resp ~ mass + mass2 + sizecat, data = data1) Cand[[7]] <- lm(resp ~ sizecat, data = data1) Cand[[8]] <- lm(resp ~ sizecat + mass + sizecat:mass, data = data1) Cand[[9]] <- lm(resp ~ agecorr + sizecat + mass + sizecat:mass, data = data1) ##create vector of model names Modnames <- paste("mod", 1:length(Cand), sep = "") aictab(cand.set = Cand, modnames = Modnames, sort = TRUE) #correct #estimate model averaged age ##no warning issued, because "age" and "agecorr" never appear in same model modavg(cand.set = Cand, parm = "age", modnames = Modnames) ##as expected, issues warning as mass occurs sometimes with "mass2" or ##"sizecatab:mass" in some of the models modavg(cand.set = Cand, parm = "mass", modnames = Modnames) ##correctly excludes models with quadratic term and interaction term ##results are CORRECT modavg(cand.set = Cand, parm = "mass", modnames = Modnames, exclude = list("mass2", "sizecat:mass"))