---
title: "Exercise 12.3 Generalized Linear Models"
output: html_document
---

<!----  

Make sure to save this file with a different name and
then delete this comment. 

Note:  this template contains code for fitting and evaluating a 
Poisson regression model. Your job is to modify it to fit a
Negative Binomial regression model.

--->

## Poisson model

Load packages

```{r message=FALSE, warning=FALSE}
library(spatstat) # for the data
library(raster) # for dealing with spatially gridded data
library(MASS) # for glm.nb used to fig neg binomial model in frequentist framework
library(R2jags) # to interface R with JAGS
library(mcmcplots) # for plotting MCMC output
library(MCMCvis) # for summarizing MCMC output  
library(ggplot2) # for plotting
options(width=150)
knitr::opts_chunk$set(fig.align = "center")
```

1. Run the code below and then comment on the appropriateness of the Poisson model in light of the goodness-of-fit test and the prediction intervals.

```{r fig.width=10}
plot(bei$x, bei$y, pch = 19, cex = 0.5,
     main = "Spatial distribution of individuals in the 50-ha Barro Colorado plot",
     xlab = "x coordinate [m]", ylab = "y coordinate [m]", frame = FALSE)
abline(h = 0, col = "grey")
abline(h = 500, col = "grey")
abline(v = 0, col = "grey")
abline(v = 1000, col = "grey")


# Coarsen the predictor data by projecting into a 50 x 50 m grid
# Take the mean of elevation within each of the 5 x 5 m grid cells:
elev <- raster(bei.extra[[1]])

# Crop the data so that they have exactly 500 x 1000 cells
ext <- extent(2.5, 1002.5, 2.5, 1002.5)
elev <- crop(elev, ext)

# Aggregate the elevation data
elev50.raster <- aggregate(elev, fact = 10, fun = mean)

# Summarize the point data by counting the number of observations in each 50 x 50 m grid cell
xy <- data.frame(x = bei$x, y = bei$y)
y.raster <- rasterize(xy, elev50.raster, fun = "count")

# Replace the NA values by 0
y.raster[is.na(y.raster)] <- 0
```

Plot the data

```{r fig.width=10.3, fig.height=2.8}
plot(stack(elev50.raster, y.raster), 
     main = c("Predictor: Mean Elevation in 50x50 m cells",
              "Response: # of Individuals in 50x50 m cells"),
       axes=FALSE)


```
```{r fig.width=6, fig.height=6}
plot(elev50.raster[], y.raster[], cex = 1, pch = 19, col = "grey", 
     ylab = "# of Individuals", xlab = "Mean Elevation [m]")  
```

Form the quadratic term

```{r }
# Scale and center elevation
elev50 <- as.numeric(scale(elev50.raster[]))

# Create elevation squared
pow.elev50 <- elev50^2 

# Form the response = number of trees in each grid cell
y <- y.raster[]
```

<!---  
-----------------------------------------------------------------------------------------

NOTE: in step 2, you will copy and then modify code below this point 


-----------------------------------------------------------------------------------------

--->


Fit the Poisson model using glm

```{r }
m.glm <- glm(y ~ elev50 + pow.elev50, family = "poisson")
summary(m.glm)
```

Look at model fit graphically   

```{r }
elev.seq <- seq(-3, 2, by = 0.05)
new.data <- data.frame(elev50 = elev.seq,
                       pow.elev50 = elev.seq ^ 2)
new.predict.glm <- predict(m.glm, newdata = new.data, type = "response")
```


```{r fig.width=6, fig.height=6}
plot(elev50, y, cex=1, col="lightgrey", pch=19,
       ylab="# of Individuals", xlab="Scaled Mean Elevation")            
lines(elev.seq, new.predict.glm, col="red", lwd=2)
```

# Now, fit the Poisson model using JAGS 

JAGS model

```{r }
model<-function(){
# priors
  beta0 ~ dnorm(0,0.001)
  beta1 ~ dnorm(0,0.001)
  beta2 ~ dnorm(0,0.001)
        
  # likelihood
  for(i in 1:N.cells){
      y[i] ~ dpois(lambda[i])
      log(lambda[i]) <- beta0 + beta1*elev50[i] + beta2*pow(elev50[i], 2)
      
      # this part is here in order to make nice prediction curves:
      prediction[i] ~ dpois(lambda[i])
          
      # Statistics for Bayesian p-value
      Presi[i] <- (y[i] - lambda[i]) / sqrt(lambda[i]) # Pearson residuals
      Presi.new[i]<-(prediction[i] - lambda[i]) / sqrt(lambda[i]) # Pearson residuals
      D[i] <- pow(Presi[i], 2)
      D.new[i] <- pow(Presi.new[i], 2)
  } 
# Add up discrepancy measures
  fit <- sum(D[])
  fit.new <- sum(D.new[])
} 
```

Gather data

```{r }
jags.data <- list(N.cells = length(y),
                    y = y,
                    elev50 = elev50)
```

Gather data nd parameters to save 

```{r }
params <- c("beta0", "beta1", "beta2", "prediction", "Presi","fit", "fit.new", "lambda")
```

Fit the model

```{r }
jm.sample <- jags.parallel(model.file = model, parameters.to.save = params,
                  data = jags.data,
                  n.chains = 3,
                  n.burnin = 2000, n.iter = 4000, DIC = TRUE)
```

Look at summary for regression parameters

```{r }
MCMCsummary(jm.sample, params = c("beta0", "beta1", "beta2"))

```

And, make sure we sampled long enough that the chains converged to the same place:

```{r fig.width=10.3, fig.height=6.8}
denplot(jm.sample, c("beta0","beta1", "beta2"))
traplot(jm.sample,  c("beta0","beta1", "beta2"))
```

Lets look at Pearson Residuals plotted against fitted values:

```{r fig.width=6, fig.height=6}
Presids <- MCMCpstr(jm.sample, params = "Presi", func = mean)
Fitted <- MCMCpstr(jm.sample, params = "lambda", func = mean)
jagsresids <- data.frame(Presids = Presids$Presi, Fitted = Fitted$lambda) 
ggplot(jagsresids, aes(Fitted, Presids)) + geom_point() + 
         geom_hline(yintercept = 0) 
```

Bayesian P-value for an overall goodness-of-fit test:

```{r fig.width=6, fig.height=6}
fitstats <- MCMCpstr(jm.sample, params = c("fit", "fit.new"), type = "chains") 
T.extreme <- fitstats$fit.new >= fitstats$fit
(p.val <- mean(T.extreme)) 
```

Let's plot predictions with credible intervals for the mean and a 95% prediction interval

```{r }
bayesCI <- MCMCpstr(jm.sample, params = "prediction",
               func = function(x) quantile(x, probs = c(0.025, 0.975)))$prediction
prds <- data.frame(y = y, sc50 = scale(elev50), low.pred=bayesCI[, 1], up.pred=bayesCI[, 2],
                   est = MCMCsummary(jm.sample, params = "prediction")$mean) 
```

```{r fig.width=6, fig.height=6}
ggplot(prds, aes(sc50, y)) + geom_point() + geom_line(aes(sc50, est)) +
  geom_line(aes(sc50, low.pred), col = "red", lty = 2) +
  geom_line(aes(sc50, up.pred), col = "red", lty = 2) +
  labs(ylab="# of Individuals", xlab="Scaled Mean Elevation") + 
  ggtitle("95% prediction interval")
```

**Comment**: comment on the appropriateness of the Poisson model in light of the goodness-of-fit test and the prediction intervals.


-------------------------------------------------------------


## Negative Binomial Model


2.  Copy the code for fitting the Poisson model using `glm` and `JAGS`. Then, modify the copied `glm` code to fit a negative binomial model with elevation and elevation$^2$ using the `glm.nb` function in the `MASS` package. 

3. Modify the JAGS code for fitting the Poisson model to allow fitting of a negative binomial model in JAGS. Calculate predicted values and a Bayesian p-value to assess goodness-of-fit of the negative binomial model. Comment on the appropriateness of the negative binomial model in light of the goodness-of-fit test and the prediction intervals.

4. Write down a set of equations describing the fitted negative binomial model.  Match any parameters in the model to the estimates obtained from JAGS. 



## Footer

R session info

```{r }
  sessionInfo()
```
 
