#LOAD PACKAGES library(lme4) library(bblme) #LOAD DATA FROM DRYAD AND SPECIFY VARIABLES THAT ARE FACTORS #ALL OTHER VARIABLES SHOULD BE NUMERICAL data<-read.csv("Squirrel Stress Plasticity Dryad.csv") data$year.f<-as.factor(data$year.f) data$breed_status<-as.factor(data$breed_status) data$assay_id<-as.factor(data$assay_id) data$squirrel_id<-as.factor(data$squirrel_id) #MODEL 1: NO RANDOM SLOPES OR RANDOM INTERCEPTS mod1<-lm(log(cort)~assay_id+breed_status+scale(poly(juldate,2))+scale(loc_density)+year.f, data,na.action=na.omit) #MODEL 2: NO RANDOM SLOPES BUT WITH RANDOM INTERCEPTS mod2<-lmer(log(cort)~assay_id+breed_status+scale(poly(juldate,2))+scale(loc_density)+year.f +(1|squirrel_id), data,na.action=na.omit,REML = F) #MODEL 3: RANDOM SLOPES AND RANDOM INTERCEPTS BUT NO COVARIANCE BETWEEN SLOPES AND INTERCEPTS mod3<-lmer(log(cort)~assay_id+breed_status+scale(poly(juldate,2))+scale(loc_density)+year.f +(1|squirrel_id)+(0+scale(loc_density)|squirrel_id), data,na.action=na.omit,REML = F) #MODEL 4: RANDOM SLOPES AND RANDOM INTERCEPTS, WITH COVARIANCE BETWEEN SLOPES AND INTERCEPT mod4<-lmer(log(cort)~assay_id+breed_status+scale(poly(juldate,2))+scale(loc_density)+year.f +(1+scale(loc_density)|squirrel_id), data,na.action=na.omit,REML = F) #GENERATE AICc TABLE bbmle::AICctab(mod1,mod2,mod3,mod4,base = T,weights = T) #GENERATE MODEL OUTPUT AND CONFIDENCE INTERVALS VIA BOOTSTRAP FOR THE TWO TOP MODELS #(IE MODEL 3 AND MODEL 4) summary(mod3) confint.merMod(mod3,method=c("boot")) summary(mod4) confint.merMod(mod4,method=c("boot"))