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")