Adding exposures to Negbin sampler

Author: HWagn

vignettes/Chapter08.Rmd

@@ -739,6 +739,15 @@ Now we analyze the road safety data allowing for unobserved
 heterogeneity. We first set up the two versions of the three-block
 MH-within-Gibbs sampler.
 
+Note that the negative binomial distribution in R is specified
+as
+$$p(y|\alpha, p)={\alpha-1+y \choose \alpha-1} p^\alpha (1-p)^y$$
+or alternatively by the parameters $\alpha$ and its expected value $$\mu=\alpha (1-p)/p.$$  The expected value of 
+$p(y|\alpha, \beta)$ his given as
+$$\mu=\alpha \frac{\frac{1}{1+\alpha/e \exp(- \mathbf{x}\boldsymbol{\beta})}}
+{\frac{\alpha/e \exp(- \mathbf{x}\boldsymbol{\beta})}{1+\alpha/e \exp(- \mathbf{x}\boldsymbol{\beta})}} = e \exp( \mathbf{x}\boldsymbol{\beta}),$$
+and we will  use  $\alpha$ and $\mu$ to specify the negative binomial distribution.
+
 ```{r}
 negbin <- function(y, X, e, b0 = 0, B0 = 100, qmean, qvar, pri.alpha,
                    full.gibbs = FALSE, burnin = 1000L, M = 50000L) {
@@ -779,12 +788,12 @@ negbin <- function(y, X, e, b0 = 0, B0 = 100, qmean, qvar, pri.alpha,
      lpri_old  <- mvtnorm::dmvnorm(beta.old,  mean = b0, sigma = B0, log = TRUE)
 
      # Compute log likelihood of proposed and old value
-     lh_proposed <- dpois(y, exp(X %*% beta.proposed),  log = TRUE)
-     lh_old  <- dpois(y, exp(X %*% beta.old), log = TRUE)
-
-     maxlik <- max(lh_old,lh_proposed)
+     lh_proposed <- dpois(y, e * exp(X %*% beta.proposed), log = TRUE)
+     lh_old  <- dpois(y, e * exp(X %*% beta.old), log = TRUE)
+    
+     maxlik <- max(lh_old, lh_proposed)
      ll <- sum(lh_proposed - maxlik) - sum(lh_old - maxlik)
-
+    
      # Compute acceptance probability and accept or not
      log_acc <- min(0, ll + lpri_proposed - lpri_old + lq_old - lq_proposed)
 
@@ -848,6 +857,7 @@ We specify the prior on $\alpha$ as the Gamma distribution with shape
 2 and rate 0.5 and use both samplers to estimate the model parameters.
 
 ```{r}
+set.seed(1234)
 pri.alpha <- data.frame(shape = 2, rate = 0.5)
 
 res1 <- negbin(y, X, e, qmean = parms.proposal$mean, qvar = parms.proposal$var,