### script for running the Bayesian analyses in Ellison 2004


### import the data
datum=read.csv("EllisData.csv")

head(datum)


# run analysis with 'frequentist methods'
results=glm(Species~Habitat+Latitude+Elevation,data=datum,family=poisson)
summary(results)
### Note that results are pseudoreplicated with respect to Latitude and Elevation

#run analysis with random blocking effect for site
library(lme4)
resultsWRandom=lmer(Species~Habitat+Latitude+Elevation+(1|Site),data=datum,
	family=poisson)
summary(resultsWRandom)

### Run analysis with non-informative prior
library(MCMCpack)
help(MCMCpoisson)
results2=MCMCpoisson(Species~Habitat+Latitude+Elevation,data=datum,mcmc=25000,
	burnin=25000,b0=0,B0=0.0001)
	### burnin and MCMC (iterations) values were taken from Table 4; b0 and B0 (mean and 1/sd of prior) was taken
	### from page 514, top of second column (in this case, all priors non-informative, so only need to least once)
summary(results2)
plot(results2) #plots the posterior distributions

### Run analysis with informative prior
results3=MCMCpoisson(Species~Habitat+Latitude+Elevation,data=datum,mcmc=25000,
	burnin=25000,b0=c(0,0.37,-0.17,-0.002),B0=c(0.0001,1,1/0.04,1/0.0003))
	### burnin and MCMC (iterations) values were taken from Table 4; b0 and B0 (mean and 1/sd of priors) was taken
	### from page 514, bottom of second column - note that b0 and B0 are lists of means and sd's
summary(results3)
plot(results3) ###plots the posterior distributions

###Run analysis with informative prior that is way off the data, but with strong certainty (low se)
results4=MCMCpoisson(Species~Habitat+Latitude+Elevation,data=datum,mcmc=25000,
	burnin=25000,b0=c(0,-5,-2,0.02),B0=c(0.0001,10,1/0.004,1/0.00003))
	### burnin and MCMC (iterations) values were taken from Table 4; b0 and B0 (mean and 1/sd of priors) was taken
	### from page 514, bottom of second column - note that b0 and B0 are lists of means and sd's
summary(results4)
plot(results4) ###plots the posterior distributions


### Lets compare a simpler model
resultsReduced=glm(Species~Habitat+Latitude,data=datum,family=poisson)
summary(resultsReduced)
#### compare AIC values and relative weights

###Now compare Bayesian analyses
results2=MCMCpoisson(Species~Habitat+Latitude+Elevation,data=datum,mcmc=25000,
	burnin=25000,b0=0,B0=0.0001,marginal.likelihood="Laplace")

results2Reduced=MCMCpoisson(Species~Habitat+Latitude,data=datum,
	mcmc=25000,burnin=25000,b0=0,B0=0.0001,marginal.likelihood="Laplace")

help(BayesFactor)

Compare=BayesFactor(results2,results2Reduced)
Compare