library(lattice) library(rethinking) ###Fit full first-order model (requires `rethinking' package)### flist <- alist( estrogen ~ dnorm(b.intercept + b.rep*rep.success + b.plastic*plastic.mass + b.cortizol*cortizol + b.age*age + b.fish*fish.mass, sigma), b.intercept ~ dnorm(0,200), b.rep ~ dnorm(0,200), b.plastic ~ dnorm(0,200), b.cortizol ~ dnorm(0,200), b.age ~ dnorm(0,200), b.fish ~ dnorm(0,200), sigma ~ dunif(0,500) ) m1 <- map(flist,data=seabirdData) ###Examine marginal posteriors and summary stats for full first-order model### precis(m1,prob=.95) post <- extract.samples(m1,n=10000) dens(post) AIC(m1) BIC(m1) ###Compute the model residuals### m1.fitted <- coef(m1)['b.intercept'] + coef(m1)['b.rep']*seabirdData$rep.success + coef(m1)['b.plastic']*seabirdData$plastic.mass + coef(m1)['b.cortizol']*seabirdData$cortizol + coef(m1)['b.age']*seabirdData$age + coef(m1)['b.fish']*seabirdData$fish.mass m1.resid <- seabirdData$estrogen - m1.fitted ###Plot residuals vs. fitted values (requires `lattice' package)### xyplot(m1.resid~m1.fitted, type=c("smooth","p"),col="black") ###Plot residuals vs. predictors### xyplot(m1.resid~seabirdData$rep.success, type=c("smooth","p"),col="black") xyplot(m1.resid~seabirdData$plastic.mass, type=c("smooth","p"),col="black") xyplot(m1.resid~seabirdData$cortizol, type=c("smooth","p"),col="black") xyplot(m1.resid~seabirdData$age, type=c("smooth","p"),col="black") xyplot(m1.resid~seabirdData$fish.mass, type=c("smooth","p"),col="black")