Data

Recall that we are interested in modeling the number of species \(S_i\) in island \(i\) as a function of island size \(A_i\). We will use the data “islands.csv”. Here is how these data look like:

setwd('U:\\uf\\courses\\bayesian course\\group activities\\7 example Poisson_Logistic reg')
dat=read.csv('islands.csv',as.is=T,fileEncoding="latin1")
dat1=dat[,c('Island','Native.bird.species','island.area')]
head(dat1)
##      Island Native.bird.species island.area
## 1 Alegranza                  21        13.4
## 2   Annobon                  10          17
## 3 Ascension                  12          97
## 4   Bermuda                  14        39.3
## 5  Boavista                  18       634.1
## 6    Branco                   9           .
dat1$island.area=as.numeric(dat1$island.area)
## Warning: NAs introduced by coercion
plot(log(dat1$island.area),dat1$Native.bird.species)

Statistical model

Our goal here is to fit the following model:

\[S_i \sim Poisson(exp(\beta_0+\beta_1 log(A_i)))\] where we assume that the priors are given by:

\[\beta_0 \sim N(0,10)\] \[\beta_1 \sim N(0,10)\]

We can succintly write the posterior distribution as:

\[p(\beta_0,\beta_1|...)\propto [\prod_i Poisson(y_i|exp(\beta_0 + \beta_1 log(A_i)))]\times N(\beta_0|0,10)N(\beta_1|0,10)\]

Analysis steps

  1. We will create JAGS code to fit the model describe above. This will be useful to compare the results from JAGS with those from our customized R code.

  2. We will create our customized R code using Metropolis-Hastings within Gibbs. To this end, we will need to:

  • create a function that evaluates the target distribution for any \(\beta_0\) and \(\beta_1\)
  • create a function that samples \(\beta_0\)
  • create a function that samples \(\beta_1\)
  • develop the full Gibbs sampler
  1. Do our results using customized R code and JAGS match? What can we conclude about our original scaling factor \(z\) (i.e., \(\beta_1\))?



Back to main menu

Comments?

Send me an email at