09.5_Batting_average.Rmd

---
title: "9.5: Batting average"
output: html_notebook
---


## Introduction

This is the "Batting average" example in Section 9.5 of Kruschke (but with some changes).


## Preliminaries

Load necessary packages:

```{r}
library(tidyverse)
library(rstan)
library(tidybayes)
library(bayesplot)
```

Set Stan to save compiled code.

```{r}
rstan_options(auto_write = TRUE)
```

Set Stan to use parallel processing where possible.

```{r}
options(mc.cores = parallel::detectCores())
```


## Data

First we import the data. (Make sure the file "BattingAverage.csv" is in your project directory.)

```{r}
ba_data <- read_csv("BattingAverage.csv")
ba_data
```

We have to tell Stan the number of categories `C`: there are nine possible positions. We'll pass the position number (`PriPosNumber`) to Stan as well.

```{r}
C <- 9
S <- NROW(ba_data) # number of subjects
N <- ba_data$AtBats # number of trials for each subject
y <- ba_data$Hits # number of successes for each subject
cats <- ba_data$PriPosNumber
```

```{r}
stan_data <- list(C = C, S = S, N = N, y = y, cats = cats)
```


## Prior

### Simulation code

We incorporate player position to allow different values of $\omega$ and $\kappa$ for each position. There will still be a parameter called $\omega_{0}$ that will serve as the mode of the beta "hyper-hyperprior" whereas $\omega_{c}$ ($1 \leq c \leq 9$) will be the mode of the hyperprior for each category of player (informing the $\theta$ for each player). The same will be true for $\kappa_{0}$ versus the collection of $\kappa_{c}$.

The gamma hyperpriors on $\kappa$ turn out to cause divergent transitions. We'll use an exponential distribution (with rate parameter 1) instead, from which Stan seems to be able to sample more easily.

We will also use vectors instead of arrays. This allows us to multiply things out without using `for` loops. For example, the lines

```
a = omega .* (kappa - 2) + 1;
b = (1 - omega) .* (kappa - 2) + 1;
```

use the operator `.*` instead of just `*`. That's because `a` and `b` are actually vectors of length 9, as are `omega` and `kappa`. So each equation actually represents 9 different equations:

```
a[1] = omega[1] * (kappa[1] - 2) + 1;
a[2] = omega[2] * (kappa[2] - 2) + 1;
a[3] = omega[3] * (kappa[3] - 2) + 1;
...
```

The sampling statement for `theta` also looks a little weird:

```
theta[s] = beta_rng(a[cats[s]], b[cats[s]]);
```

This is sampling 948 values, but `a[cats[s]]` and `b[cats[s]]` refer to only the nine possible values of `a` or `b` depending on the category (position) of the player for whom we're sampling `theta`.


```{stan, output.var = "ba_prior", cache = TRUE}
data {
    int<lower = 0> C;
    int<lower = 0> S;
    array[S] int<lower = 0> N;
    array[S] int<lower = 0> y;
    array[S] int<lower = 0> cats;
}
transformed data {
    real<lower = 0> A_omega_0;
    real<lower = 0> B_omega_0;
    real<lower = 0> R_kappa_0;

    A_omega_0 = 2;       // hyperprior parameters for omega
    B_omega_0 = 2;
    R_kappa_0 = 1;       // hyperprior parameter for kappa
}
generated quantities {
    real<lower = 0, upper = 1> omega_0;
    vector<lower = 0, upper = 1>[C] omega;
    vector<lower = 0>[C] kappa_minus_two;
    real<lower = 0> kappa_minus_two_0;
    real<lower = 2> kappa_0;
    vector<lower = 2>[C] kappa;
    real<lower = 0> a_0;
    real<lower = 0> b_0;
    vector<lower = 0>[C] a;
    vector<lower = 0>[C] b;
    vector<lower = 0, upper = 1>[S] theta;
    real omega_7_4;
    real omega_4_1;
    real theta_75_156;
    real theta_159_844;
    real theta_494_754;
    real theta_573_428;
    array[S] int<lower = 0> y_sim;

    omega_0 = beta_rng(A_omega_0, B_omega_0);
    kappa_minus_two_0 = exponential_rng(R_kappa_0);
    kappa_0 = kappa_minus_two_0 + 2;
    

    a_0 = omega_0 * (kappa_0 - 2) + 1;
    b_0 = (1 - omega_0) * (kappa_0 - 2) + 1;
    for(c in 1:C) {
        omega[c] = beta_rng(a_0, b_0);
    }
    
    for(c in 1:C) {
        kappa_minus_two[c] = exponential_rng(R_kappa_0);
        kappa[c] = kappa_minus_two[c] + 2;
    }
    
    
    a = omega .* (kappa - 2) + 1;       // Componentwise mult.
    b = (1 - omega) .* (kappa - 2) + 1; // of vectors

    for(s in 1:S) {
        theta[s] = beta_rng(a[cats[s]], b[cats[s]]);
    }
    
    omega_7_4 = omega[7] - omega[4];
    omega_4_1 = omega[4] - omega[1];
    theta_75_156 = theta[75] - theta[156];
    theta_159_844 = theta[159] - theta[844];
    theta_494_754 = theta[494] - theta[754];
    theta_573_428 = theta[573] - theta[428];

    y_sim = binomial_rng(N, theta); // vectorized
}
```

```{r, cache = TRUE}
fit_ba_prior <- sampling(ba_prior,
                         data = stan_data,
                         chains = 1,
                         algorithm = "Fixed_param",
                         seed = 11111,
                         refresh = 0)
```

```{r}
samples_ba_prior <- tidy_draws(fit_ba_prior)
samples_ba_prior
```

### Examine prior

```{r}
mcmc_hist(samples_ba_prior,
          pars = "omega_0")
```

```{r}
mcmc_hist(samples_ba_prior,
          pars = vars(starts_with("omega[")))
```

```{r}
mcmc_hist(samples_ba_prior,
          pars = "kappa_0")
```

```{r}
mcmc_hist(samples_ba_prior,
          pars = vars(starts_with("kappa[")))
```

```{r}
mcmc_pairs(fit_ba_prior,
           pars = c("omega_0", "kappa_0"))
```

```{r}
mcmc_hist(samples_ba_prior,
          pars = c("theta[75]", "theta[156]",
                   "theta[159]", "theta[844]",
                   "theta[494]", "theta[754]",
                   "theta[573]", "theta[428]"))
```

```{r}
mcmc_hist(samples_ba_prior,
          pars = c("omega_7_4", "omega_4_1",
                   "theta_75_156", "theta_159_844",
                   "theta_494_754", "theta_573_428"))
```

```{r}
mcmc_pairs(fit_ba_prior,
           pars = c("omega[1]", "omega[4]", "omega[7]"))
```

### Prior predictive distribution

Use proportions for select players:

```{r}
samples_ba_prior <- samples_ba_prior %>%
    mutate(`y_sim_prop[75]` =  `y_sim[75]`/N[75],
           `y_sim_prop[156]` = `y_sim[156]`/N[156],
           `y_sim_prop[159]` = `y_sim[159]`/N[159],
           `y_sim_prop[844]` = `y_sim[844]`/N[844],
           `y_sim_prop[494]` = `y_sim[494]`/N[494],
           `y_sim_prop[754]` = `y_sim[754]`/N[754],
           `y_sim_prop[573]` = `y_sim[573]`/N[573],
           `y_sim_prop[428]` = `y_sim[428]`/N[428])
```

```{r}
y_sim_ba_prior <- samples_ba_prior %>%
    dplyr::select(starts_with("y_sim_prop")) %>%
    as.matrix()
```

We can look at the data for just the eight players from the book:

```{r}
ppd_intervals(ypred = y_sim_ba_prior)
```


## Model

```{stan, output.var = "ba", cache = TRUE}
data {
    int<lower = 0> C;
    int<lower = 0> S;
    array[S] int<lower = 0> N;
    array[S] int<lower = 0> y;
    array[S] int<lower = 0> cats;
}
transformed data {
    real<lower = 0> A_omega_0;
    real<lower = 0> B_omega_0;
    real<lower = 0> R_kappa_0;

    A_omega_0 = 2;       // hyperprior parameters for omega
    B_omega_0 = 2;
    R_kappa_0 = 1;       // hyperprior parameter for kappa
}
parameters {
    real<lower = 0, upper = 1> omega_0;
    vector<lower = 0, upper = 1>[C] omega;
    real<lower = 0> kappa_minus_two_0;
    vector<lower = 0>[C] kappa_minus_two;
    vector<lower = 0, upper = 1>[S] theta;
}
transformed parameters {
    real<lower = 2> kappa_0;
    vector<lower = 2>[C] kappa;
    real<lower = 0> a_0;
    real<lower = 0> b_0;
    vector<lower = 0>[C] a;
    vector<lower = 0>[C] b;
    
    kappa_0 = kappa_minus_two_0 + 2;     
    kappa = kappa_minus_two + 2;     
    
    a_0 = omega_0 * (kappa_0 - 2) + 1;
    b_0 = (1 - omega_0) * (kappa_0 - 2) + 1;
    a = omega .* (kappa - 2) + 1;       // Componentwise mult.
    b = (1 - omega) .* (kappa - 2) + 1; // of vectors
}
model {
    omega_0 ~ beta(A_omega_0, B_omega_0);
    kappa_minus_two_0 ~ exponential(R_kappa_0);
    omega ~ beta(a_0, b_0);
    kappa_minus_two ~ exponential(R_kappa_0);
    theta ~ beta(a[cats], b[cats]);              
    y ~ binomial(N, theta);          // likelihood
}
generated quantities {
    real omega_7_4;
    real omega_4_1;
    real theta_75_156;
    real theta_159_844;
    real theta_494_754;
    real theta_573_428;
    array[S] int<lower = 0> y_sim;
    
    omega_7_4 = omega[7] - omega[4];
    omega_4_1 = omega[4] - omega[1];
    theta_75_156 = theta[75] - theta[156];
    theta_159_844 = theta[159] - theta[844];
    theta_494_754 = theta[494] - theta[754];
    theta_573_428 = theta[573] - theta[428];
    
    y_sim = binomial_rng(N, theta);
}
```

```{r, cache = TRUE}
fit_ba <- sampling(ba,
                   data = stan_data,
                   seed = 11111,
                   refresh = 0)
```

```{r}
samples_ba <- tidy_draws(fit_ba)
samples_ba
```

Again, we'll compute proportions:

```{r}
samples_ba <- samples_ba %>%
    mutate(`y_sim_prop[75]` =  `y_sim[75]`/N[75],
           `y_sim_prop[156]` = `y_sim[156]`/N[156],
           `y_sim_prop[159]` = `y_sim[159]`/N[159],
           `y_sim_prop[844]` = `y_sim[844]`/N[844],
           `y_sim_prop[494]` = `y_sim[494]`/N[494],
           `y_sim_prop[754]` = `y_sim[754]`/N[754],
           `y_sim_prop[573]` = `y_sim[573]`/N[573],
           `y_sim_prop[428]` = `y_sim[428]`/N[428])
```


## Model diagnostics

```{r}
stan_trace(fit_ba,
           pars = c("omega_0", "kappa_0"))
```

```{r}
stan_trace(fit_ba,
           pars = c("omega"))
```

```{r}
stan_trace(fit_ba,
           pars = c("theta[75]", "theta[156]",
                    "theta[159]", "theta[844]",
                    "theta[494]",  "theta[754]",
                    "theta[573]", "theta[428]"))
```

```{r}
mcmc_acf(fit_ba,
         pars = "omega_0", "kappa_0")
```


```{r}
mcmc_acf(fit_ba,
         pars = c("omega[1]", "omega[2]", "omega[3]"))
```

```{r}
mcmc_acf(fit_ba,
         pars = c("omega[4]", "omega[5]", "omega[6]"))
```

```{r}
mcmc_acf(fit_ba,
         pars = c("omega[7]", "omega[8]", "omega[9]"))
```

```{r}
mcmc_acf(fit_ba,
         pars = c("kappa[1]", "kappa[2]", "kappa[3]"))
```

```{r}
mcmc_acf(fit_ba,
         pars = c("kappa[4]", "kappa[5]", "kappa[6]"))
```

```{r}
mcmc_acf(fit_ba,
         pars = c("kappa[7]", "kappa[8]", "kappa[9]"))
```

```{r}
mcmc_rhat(rhat(fit_ba))
```

```{r}
mcmc_neff(neff_ratio(fit_ba))
```


## Model summary

```{r}
print(fit_ba,
      pars = c("omega_0", "kappa_0", "omega", "kappa"))
```


## Model visualization

```{r}
mcmc_areas(fit_ba,
           pars = "omega_0")
```

```{r}
mcmc_areas(fit_ba,
           pars = vars(starts_with("omega[")))
```

```{r}
mcmc_areas(fit_ba,
           pars = "kappa_0")
```

```{r}
mcmc_areas(fit_ba,
           pars = vars(starts_with("kappa[")))
```

```{r}
mcmc_areas(fit_ba,
           pars = c("theta[75]", "theta[156]",
                   "theta[159]", "theta[844]",
                   "theta[494]", "theta[754]",
                   "theta[573]", "theta[428]"))
```

```{r}
pairs(fit_ba,
      pars = c("omega_0", "kappa_0"))
```

Compare different positions. Note that these don't match the book graphs because the book has an error here. The numbers indexing `omega` do not match the position descriptions in the book. In the actual data, Position 1 is Pitcher, Position 4 is 2nd Base, and Position 7 is Left Field.

```{r}
mcmc_areas(fit_ba,
           pars = "omega_7_4", "omega_4_1")
```

Comparisons of individual players:

```{r}
mcmc_areas(fit_ba,
           pars = c("theta_75_156", "theta_159_844",
           "theta_494_754", "theta_573_428"))
```


## Posterior predictive check

```{r}
y_sim_ba <- samples_ba %>%
    dplyr::select(starts_with("y_sim_prop")) %>%
    as.matrix()
```

```{r}
ppc_intervals(y = c(y[75]/N[75], y[156]/N[156],
                    y[159]/N[159], y[844]/N[844],
                    y[494]/N[494], y[754]/N[754],
                    y[573]/N[573], y[428]/N[428]),
              yrep = y_sim_ba)
```

There is some shrinkage visible. It would be more prominent for people with low values of `AtBats`. (Less data means the prior is more powerful.) For those with high `AtBats`, the data speaks loudly, so the posterior data draws are centered pretty much right at the data value.