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

3. 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