minor changes Poisson

Author: HWagn

vignettes/Chapter08.Rmd

@@ -566,15 +566,15 @@ knitr::kable(round(res_probit.labor * pi / sqrt(3), 3))
 
 ### Example 8.8: Road safety data
 
-We fit two different Poisson regression models to the series of
+We will fit two different Poisson regression models to the series of
 monthly deadly and seriously injured children aged 6-10 in Linz
 introduced in Example 2.1:
 
 1. a small model with intercept, intervention effect and holiday dummy
    (activated in July/August);
 
 2. a larger model with intercept, intervention effect, 
-   and a seasonal dummy variables for all months except for December.
+   and a seasonal dummy variables for all months except December.
 
 The sampling performance for these two models is assessed to study how
 the acceptance rate deteriorates, when the dimension of regression effects $d$ 
@@ -607,9 +607,7 @@ gen.proposal.poisson <- function(y, X, e, b0 = 0, B0 = 100, t.max = 20) {
   betas <- matrix(NA_real_, ncol = t.max, nrow = d)
   beta.new <- matrix(c(log(mean(y)), rep(0, d - 1)), nrow = d) 
   
-  b0 <- matrix(rep(b0, length.out = d), nrow = d) 
-  B0.inv <- diag(rep(1 / B0, length.out = d), nrow = d)
-
+  B0.inv=solve(B0)
   for (t in seq_len(t.max)) {
     beta.old <- beta.new
     
@@ -638,17 +636,16 @@ gen.proposal.poisson <- function(y, X, e, b0 = 0, B0 = 100, t.max = 20) {
 
 We use a rather flat normal independence prior
 $\mathcal{N}(\mathbf{0}, 100 \mathbf{I})$ on the regression
-effects. First we determine the parameters of the proposal
+effects and determine the parameters of the proposal
 distribution.
 
 ```{r}
-parms.proposal <- gen.proposal.poisson(y, X, e, b0 = 0, B0 = 100)
+d=ncol(X)
+parms.proposal <- gen.proposal.poisson(y, X, e, b0 = rep(0,d), B0 =diag(100,d))
 parms.proposal
 ```
 
-To set up the independence Metropolis-Hastings algorithm for 
-the Poisson model, we first write a short program for the MH step
-to sample the regression effects.
+To implement the independence Metropolis-Hastings algorithm  we  write a short program for the MH step for sampling the regression effects.
 ```{r}
 sample_beta<- function(y,X,e, b0, B0, qmean, qvar, beta.old){
  
@@ -686,10 +683,11 @@ sample_beta<- function(y,X,e, b0, B0, qmean, qvar, beta.old){
 }
 ```
 
-We use this program to sample from the posterior.
+Next we combine the determination of the proposal and the MH-sampling step 
+in a  program to sample from the posterior of a Poisson regression model.
+
 ```{r}
-poisson <- function(y, X, e, b0 = 0, B0 = 100, qmean, qvar,
-                    burnin = 1000L, M = 10000L) {
+poisson <- function(y, X, e, b0 = 0, B0 = 100, burnin = 1000L, M = 10000L) {
   d <- ncol(X)
 
   b0 <- rep(b0, length.out = d)
@@ -699,6 +697,10 @@ poisson <- function(y, X, e, b0 = 0, B0 = 100, qmean, qvar,
   colnames(beta.post) <- colnames(X)
   acc <- numeric(length = M)
   
+  parms.proposal<- gen.proposal.poisson(y, X, e, b0 , B0)
+  qmean <- parms.proposal$mean
+  qvar<-parms.proposal$var
+  
   beta <- as.vector(mvtnorm::rmvnorm(1, mean = qmean, sigma = qvar))
   
   for (m in seq_len(burnin + M)) {
@@ -719,8 +721,7 @@ We perform MCMC and report the results.
 
 ```{r}
 set.seed(1234)
-res1 <- poisson(y, X, e, b0 = 0, B0 = 100,
-                qmean = parms.proposal$mean, qvar = parms.proposal$var)
+res1 <- poisson(y, X, e, b0 = 0, B0 = 100)
 
 res.poisson1 <- cbind(t(round(apply(res1$beta.post, 2, res.mcmc), 3)),
                       "exp(beta)" = round(exp(colMeans(res1$beta.post)), 5))
@@ -750,22 +751,14 @@ seas <- rbind(diag(1, 11), rep(0, 11))
 seas.dummies <- matrix(rep(t(seas), 16), ncol = 11, byrow = TRUE)
 colnames(seas.dummies) <- c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul",
                             "Aug", "Sep", "Oct", "Nov")
-X.large <- cbind(X[, -3],
-                 seas.dummies)
+X.large <- cbind(X[, -3], seas.dummies)
 ```
 
-We set the prior parameters and compute the parameters of the proposal
-distribution.
-
-```{r}
-parms.proposal2 <- gen.proposal.poisson(y, X.large, e, b0 = 0, B0 = 100)
-```
-Next we fit the model.
+We set the prior parameters and fit the model.
 
 ```{r}
 set.seed(1234)
-res2 <- poisson(y, X.large, e, b0 = 0, B0 = 100,
-                qmean = parms.proposal2$mean, qvar = parms.proposal2$var)
+res2 <- poisson(y, X.large, e, b0 = 0, B0 = 100)
 
 res.poisson2 <- cbind(t(round(apply(res2$beta.post, 2, res.mcmc), 3)),
                       "exp(beta)" = round(exp(colMeans(res2$beta.post)), 5))
@@ -797,7 +790,7 @@ be estimated.
 ### Example 8.9: Road safety data
 
 Now we re-analyze the road safety data allowing for unobserved
-heterogeneity. We first set up the two versions of the three-block
+heterogeneity. We will first set up the two versions of the three-block
 MH-within-Gibbs sampler.