misc

_common.R

@@ -77,3 +77,34 @@ knit_print.stanmodel <- function(x, options) {
 #   print(head(x$summary[ , stats, drop = FALSE], n))
 # }
 
+autoscale <- function(x, ...) {
+  UseMethod("autoscale")
+}
+
+autoscale.numeric <- function(x, center = TRUE, scale = TRUE) {
+  nvals <- length(unique(x))
+  if (nvals <= 1) {
+    out <- x
+  } else if (nvals == 2) {
+    out <- if (scale) {
+      (x - min(x, na.rm = TRUE)) / diff(range(x, finite = TRUE))
+    } else x
+    if (center) {
+      out <- x - mean(x)
+    }
+  } else {
+    out <- if (center) {
+      x - mean(x, na.rm = TRUE)
+    } else x
+    out <- if (scale) out / sd(out, na.rm = TRUE)
+  }
+  out
+}
+
+autoscale.matrix <- function(x, center = TRUE, scale = TRUE) {
+  apply(x, 2, autoscale, center = center, scale = scale)
+}
+
+autoscale.data_frame <- function(x, center = TRUE, scale = TRUE) {
+  mutate_if(is.numeric, autoscale.numeric, center = center, scale = scale)
+}

binomial.Rmd

@@ -175,25 +175,7 @@ library("magrittr")
 
 URL <- "https://raw.githubusercontent.com/carlislerainey/priors-for-separation/master/br-replication/data/need.csv"
 
-autoscale <- function(x, center = TRUE, scale = TRUE) {
-  nvals <- length(unique(x))
-  if (nvals <= 1) {
-    out <- x
-  } else if (nvals == 2) {
-    out <- if (scale) {
-      (x - min(x, na.rm = TRUE)) / diff(range(x, finite = TRUE))
-    } else x
-    if (center) {
-      out <- x - mean(x)
-    }
-  } else {
-    out <- if (center) {
-      x - mean(x, na.rm = TRUE)
-    } else x
-    out <- if (scale) out / sd(out, na.rm = TRUE)
-  }
-  out
-}
+
 
 f <- (oppose_expansion ~ dem_governor + obama_win + gop_leg + percent_uninsured +
       income + percent_nonwhite + percent_metro)

categorical.Rmd

@@ -1,2 +1,10 @@
 # Categorical Variables
 
+## Example: Mexico Vote Choice
+
+This example from [Zelig](http://docs.zeligproject.org/en/latest/zeligchoice-mlogit.html) uses a multinomial model 
+
+```{r}
+data("mexico", package = "Zelig")
+```
+

counts.Rmd

@@ -115,10 +115,28 @@ In **rstanarm** use `r rdoc("rstanarm", "stan_glm")` for a Poisson GLM.
 
 ## Example: Bilateral Sanctions
 
-A regression model of bilateral sanctions for the period 1939 to 1983. 
+A regression model of bilateral sanctions for the period 1939 to 1983 [@Martin1992a].
 The outcome variable is the number of countries imposing sanctions.
+
+```{r results='hide'}
+mod_poisson1 <- stan_model("stan/poisson1.stan")
+```
+
 ```{r}
 data("sanction", package = "Zelig")
+
+f <- num ~ coop + target -1L
+reg_model_data <- lm_preprocess(f, data = sanction)
+
+sanction_data <-
+  list(X = autoscale(reg_model_data$X),
+       y = reg_model_data$y)
+sanction_data$N <- nrow(sanction_data$X)
+sanction_data$K <- ncol(sanction_data$X)
+```
+
+```{r result='hide'}
+fit_sanction_pois <- sampling(mod_poisson1, data = sanction_data)
 ```
 
 
@@ -134,7 +152,7 @@ The Negative Binomial distribution has an additional parameter that allows the v
 
 The outcome is modeled as a negative binomial distribution,
 $$
-y_i \sim \dbinom(\alpha_i, \beta)
+y_i \sim \dnegbinom(\alpha_i, \beta)
 $$
 with shape $\alpha \in \R^{+}$ and inverse scale $\beta \in \R^{+}$, and $\E(y) = \alpha_i / \beta$ and $\Var(Y) = \frac{\alpha_i}{\beta^2}(\beta + 1)$.
 Then the mean can be modeled and transformed to the 
@@ -167,7 +185,34 @@ In Stan, there are multiple parameterizations of the
 
 Also, `rstanarm` supports Poisson and [negative binomial models](https://cran.r-project.org/web/packages/rstanarm/vignettes/count.html).
 
-- @BDA3 [Ch 16]
+### Example: Economic Sanctions II {ex-econ-sanctions-2}
+
+Continuing the [economic sanctions example](#ex-econ-sanctions-2) of @Martin1992a.
+
+
+```{r results='hide'}
+mod_negbin1 <- stan_model("stan/negbin1.stan")
+```
+
+
+```{r result='hide'}
+fit_sanction_nb <- sampling(mod_negbin1, data = sanction_data, control = list(adapt_delta  = 0.95))
+```
+
+```{r}
+summary(fit_sanction_nb, par = c("a", "b", "phi"))$summary
+```
+
+We can compare the 
+```{r message=FALSE}
+loo_sanction_pois <- loo(extract_log_lik(fit_sanction_pois, "log_lik"))
+loo_sanction_nb <- loo(extract_log_lik(fit_sanction_nb, "log_lik"))
+```
+
+```{r}
+loo::compare(loo_sanction_pois, loo_sanction_nb)
+```
+
 
 ### References
 
@@ -177,4 +222,3 @@ For general references on count models see
 - @McElreath2016a [Ch 10]
 - @Fox2016a [Ch. 14]
 - @BDA3 [Ch. 16]
-

includes/before_body.html

@@ -39,11 +39,10 @@
 \newcommand{\dlogistic}{\mathrm{Logistic}}
 \newcommand{\dllogistic}{\mathrm{Log-logistic}}
 \newcommand{\dpois}{\mathrm{Poisson}}
-\newcommand{\dbinom}{\mathrm{Bin}}
-\newcommand{\dbinom}{\mathrm{Bin-alt}}
+\newcommand{\dbin}{\mathrm{Bin}}
 \newcommand{\dmultinom}{\mathrm{Multinom}}
 \newcommand{\dnbinom}{\mathrm{Neg-bin}}
-\newcommand{\dnbinomalt}{\mathrm{Neg-bin-alt}}
+\newcommand{\dnbinomalt}{\mathrm{Neg-bin2}}
 \newcommand{\dbetabinom}{\mathrm{Beta-bin}}
 \newcommand{\dcauchy}{\mathrm{Cauchy}}
 \newcommand{\dhalfcauchy}{\mathrm{Cauchy}^{+}}

stan/mlogit.stan

@@ -0,0 +1,40 @@
+data {
+  int<lower = 0> N;
+  // categories
+  int<lower = 2> M;
+  int<lower = 1, upper = M> y[N];
+  // covariates
+  int<lower = 0> K;
+  vector[K] X[N];
+}
+transformed data {
+  vector<lower = 0.>[M] prior_scale_a;
+  vector<lower = 0.>[K] prior_scale_b[M];
+  prior_scale_a = rep_vector(M, 10.);
+  for (m in 1:M) {
+    prior_scale_b[m] = rep_vector(K, 2.5);
+  }
+}
+parameters {
+  vector[M] a;
+  vector[M, K] b;
+}
+transformed parameters {
+  vector[M] eta[N];
+  for (n in 1:N) {
+    eta[n] = a + X[n] * b;
+  }
+}
+model {
+  a ~ normal(0, prior_scale_a)
+  for (m in 1:M) {
+    b[m] ~ normal(0, prior_scale_b[m]);
+  }
+  for (n in 1:N) {
+    y[n] ~ ctegorical_logit(eta);
+  }
+}
+generated quantities {
+  // log_lik
+  // yrep
+}

stan/negbin1.stan

@@ -23,24 +23,27 @@ parameters {
   real<lower = 0.> reciprocal_phi;
 }
 transformed parameters {
-  vector<lower = 0., upper = 1.>[N] mu;
+  vector[N] eta;
   real<lower = 0.> phi;
   phi = 1. / reciprocal_phi;
-  mu = exp(a + X * b);
+  eta = a + X * b;
 }
 model {
   // priors
-  a ~ normal(0., 5.);
+  a ~ normal(0., 10.);
   b ~ normal(0., 2.5);
   reciprocal_phi ~ cauchy(0., 5.);
   // likelihood
-  y ~ neg_binomial_2(mu, phi);
+  y ~ neg_binomial_2_log(eta, phi);
 }
 generated quantities {
+  // expected value
+  vector[N] mu;
   // simulate data from the posterior
   vector[N] y_rep;
   // log-likelihood posterior
   vector[N] log_lik;
+  mu = exp(eta);
   for (i in 1:N) {
     y_rep[i] = neg_binomial_2_rng(mu[i], phi);
     log_lik[i] = neg_binomial_2_lpmf(y[i] | mu[i], phi);