Changes in reporting results

Author: HWagn

vignettes/Chapter08.Rmd

@@ -112,22 +112,20 @@ To compute summary statistics from the  posterior we use the following function.
 
 ```{r}
 res.mcmc <- function(x, lower = 0.025, upper = 0.975) {
-  y <- c(quantile(x, lower), mean(x), quantile(x, upper))
+  y <- cbind(quantile(x, lower), mean(x), quantile(x, upper))
   names(y) <- c(paste0(lower * 100, "%"), "Mean", paste0(upper * 100, "%"))
   y
  }    
 ```
 
-
 ```{r}
 res_beta <- apply(betas, 2, res.mcmc)
-knitr::kable(round(res_beta, 3))
+knitr::kable(round(t(res_beta), 3))
 
 (p_unemploy_base <- pnorm(res_beta[2, 1]))
 ```
-The estimated risk of unemployment for a baseline person is low and it
-is even lower for a white collar worker. It is higher for females,
-older persons and particularly for those unemployed in 1997.
+
+The estimated risk of unemployment for a baseline person is low and it is even lower for a white collar worker. It is higher for females, older persons and particularly for those unemployed in 1997.
 
 ```{r, echo = -c(1:2)}
 if (pdfplots) {
@@ -140,8 +138,7 @@ for (j in seq_len(ncol(betas))) {
 }
 ```
 
-A plot of the autocorrelation of the draws shows that although there is some
-autocorrelation, it vanishes after a few lags.
+A plot of the autocorrelation of the draws shows that although there is some autocorrelation, it vanishes after a few lags.
 
 ```{r, echo = -c(1:2)}
 if (pdfplots) {
@@ -152,12 +149,13 @@ for (j in seq_len(ncol(betas))) {
   acf(betas[, j], main = "", xlab = colnames(betas)[j], ylab = "")
 }
 ```
-The sampler is easy to implement, however there  might be problems  when the
-response variable  contains either only few or very many successes.
+The sampler is easy to implement, however there  might be problems  when the response variable  contains either only few or very many successes.
 To illustrate this issue, we use data where in $N = 500$ trials  only 1 success
 or only 1 failure is observed.
 
 ```{r}
+set.seed(1234)
+
 N <- 500
 X <- matrix(1, nrow = N)
 
@@ -168,7 +166,8 @@ y2 <- c(rep(0, N-1), 1)
 betas2 <- probit(y2, X, b0 = 0, B0 = 10000) 
 ```
 
-In both cases the autocorrelation of the draws decreases very slowly.
+In both cases the autocorrelation of the draws decreases very slowly and remains high even higher lags.
+
 ```{r, echo = -c(1:2)}
 if (pdfplots) {
   pdf("8-1_3.pdf", width = 8, height = 5)
@@ -180,13 +179,9 @@ plot(betas2, type = "l", main = "", xlab = "", ylab = "")
 acf(betas2)
 ```
 
-High autocorrelated draws in probit models not only occur if successes
-or failures are very rare, but also when a covariate (or a linear
-combination of covariates) perfectly allows to predict successes
+High autocorrelation in MCMC draws for probit models not only occur if successesor failures are very rare, but also when a covariate (or a linear combination of covariates) perfectly  allows to predict successes
 and/or failures.  Complete separation means that both successes and
-failures can be perfectly predicted by a covariate, whereas with
-quasi-complete separation only either successes or failures can be
-predicted perfectly.
+failures can be perfectly predicted by a covariate, whereas quasi-complete separation means that  either successes or failures can be predicted perfectly.
 
 # Example  8.3
 
@@ -202,14 +197,15 @@ y <- rep(c(0,1), c(ns, N - ns))
 table(x,y)
 ```
 
-We estimate the  model parameters and plot the ACF of the draws. Again the
-autocorrelations remain high even for lag 35.
+We estimate the  model parameters and plot the ACF of the draws. Again the autocorrelations remain high even for lag 35.
 
 ```{r, echo = -c(1:2)}
 if (pdfplots) {
   pdf("8-1_4.pdf", width = 8, height = 5)
 }
 par(mfrow = c(2, 2), mar = c(2.5, 1.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
+
+set.seed(1234)
 X <- cbind(rep(1, N), x)
 betas <- probit(y, X, b0 = 0, B0 = 10000)
 
@@ -220,7 +216,8 @@ acf(betas[, 2])
 ```
 # Example 8.4
 
-To illustrate quasi-seperation we use the responses as in Example 8.3., but now $x=0$ for all successes  and additionally for 100 failures.
+To illustrate quasi-seperation we use the responses as in Example 8.3., but now $x=1$ for all successes  and additionally for 100 failures. Hence  for $x=1$ always a success is observed, whereas for $x=0$ both successes and failures occur,
+
 ```{r}
 x <- rep(c(0,1), c(ns-100, N - ns+100))
 table(x, y)
@@ -234,6 +231,8 @@ if (pdfplots) {
 }
 
 par(mfrow = c(2, 2), mar = c(2.5, 1.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
+
+set.seed(1234)
 X <- cbind(rep(1, N), x)
 betas <- probit(y, X, b0 = 0, B0 = 10000)
 
@@ -243,36 +242,34 @@ plot(betas[, 2], type = "l", main = "", xlab = "", ylab = "")
 acf(betas[, 2])
 ```
 
-If we change the setting so that  x  takes values of $0$ not only for failures but also for some successes, the 
+If we change the setting so that  x  takes values of $1$ not only for failures but also for some successes, the 
 autocorrelations are low for the intercept but still high for the covariate effect.
+
 ```{r, echo = -c(1:2)}
 if (pdfplots) {
   pdf("8-1_5a.pdf", width = 8, height = 5)
 }
+
 x <- rep(c(0,1), c(ns+100, N - ns-100))
 table(x, y)
 
-par(mfrow = c(2, 2), mar = c(2.5, 1.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
+set.seed(1234)
 X <- cbind(rep(1, N), x)
 betas <- probit(y, X, b0 = 0, B0 = 10000)
 
+par(mfrow = c(2, 2), mar = c(2.5, 1.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
 plot(betas[, 1], type = "l", main = "", xlab = "", ylab = "")
 acf(betas[, 1])
 plot(betas[, 2], type = "l", main = "", xlab = "", ylab = "")
 acf(betas[, 2])
 ```
 
-
-High autocorrelations typically indicate problems with the sampler. If
-there is complete or quasi-complete separation in the data, the
-likelihood is monotone and the maximum likelihood estiamte does not
-exist.  In a Bayesian approach using a flat, improper prior on the
+High autocorrelations typically indicate problems with the sampler. If there is complete or quasi-complete separation in the data, the likelihood is monotone and the maximum likelihood estimate does not exist.  In a Bayesian approach using a flat, improper prior on the
 regression effects will result in an improper posterior
 distribution. Hence, a proper prior is required to avoid improper
 posteriors in case of separation.
 
-In the examples above we used a proper prior which is rather
-flat. With a more informative prior, the autocorrelations of the draws
+In the examples above we used a very flat but proper prior  With a more informative prior, the autocorrelations of the draws
 are lower. This can be seen in the next figure, where the simulated
 data under quasi-separation are re-analyzed with a Normal prior that
 is tighter around zero.
@@ -281,9 +278,11 @@ is tighter around zero.
 if (pdfplots) {
   pdf("8-1_6.pdf", width = 8, height = 5)
 }
-par(mfrow = c(2, 2), mar = c(2.5, 1.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
+set.seed(1234)
 betas <- probit(y, X, b0 = 0, B0 = 2.5^2)
 
+par(mfrow = c(2, 2), mar = c(2.5, 1.5, 1.5, .1), mgp = c(1.5, .5, 0), lwd = 1.5)
+
 plot(betas[, 1], type = "l", main = "", xlab = "", ylab = "")
 acf(betas[, 1])
 plot(betas[, 2], type = "l", main = "", xlab = "", ylab = "")
@@ -668,18 +667,17 @@ pri.alpha <- data.frame(shape = 2, rate = 0.5)
 res1 <- negbin(y, X, e, qmean = parms.proposal$mean, qvar = parms.proposal$var,
                pri.alpha = pri.alpha, full.gibbs = TRUE)
 
-res.negbin.full <- cbind(apply(res1$beta.post, 2, res.mcmc), 
+res.negbin.full <- rbind(t(apply(res1$beta.post, 2, res.mcmc)), 
                          res.mcmc(res1$alpha.post))
-colnames(res.negbin.full)[4] <- "alpha"
+rownames(res.negbin.full)[4] <- "alpha"
 knitr::kable(round(res.negbin.full, 3))
 
                
 res2 <- negbin(y, X, e, qmean = parms.proposal$mean, qvar = parms.proposal$var,
                pri.alpha = pri.alpha, full.gibbs = FALSE)
 
-res.negbin.partial <- cbind(apply(res2$beta.post, 2, res.mcmc), 
-                         res.mcmc(res2$alpha.post))
-colnames(res.negbin.partial)[4] <- "alpha"
+res.negbin.partial <- rbind(t(apply(res2$beta.post, 2, res.mcmc)),                          res.mcmc(res2$alpha.post))
+rownames(res.negbin.partial)[4] <- "alpha"
 knitr::kable(round(res.negbin.partial, 3))
 ```
 As expected estimation results using both samplers are rather similar.