16.2_Smart_Drug_Robust.Rmd

---
title: "16.2: Smart Drug (robust estimation)"
output: html_notebook
---


## Introduction

This is the example from Section 16.2 of Kruschke on the (fictitious) effect of the "smart drug", this time using a Student t model for robust estimation. Make sure the "TwoGroupIQ.csv" file is in your project directory.


## 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

```{r}
IQ_data <- read_csv("TwoGroupIQ.csv")
IQ_data
```

For this 1-sample test, we will only look at those who took the smart drug.

```{r}
IQ1_data <- IQ_data %>%
    filter(Group == "Smart Drug")
IQ1_data
```

```{r}
ggplot(IQ1_data, aes(x = Score)) +
    geom_histogram(binwidth = 10, boundary = 50)
```

```{r}
N <- NROW(IQ1_data)
y <- IQ1_data$Score
stan_data <- list(N = N, y = y)
```


## Prior

### Simulation code

The book uses `nu_minus_one` with an exponential prior. Then we will add one to to `nu_minus_one` to constrain `nu` to take values larger than 1. For the rate `R_nu_minus_one`, we'll use `1.0/29` as this implies an "average" of 30 degrees of freedom (which is when the Student t model starts looking indistinguishable from a normal model).

```{stan, output.var = "IQ_robust_prior", cache = TRUE}
data {
    int<lower = 0> N;
    array[N] real y;
}
transformed data {
    real M;
    real<lower = 0> S;
    real<lower = 0> R_nu_minus_one;
    real<lower = 0> R_sigma;
    
    M = 100;   // mean centered at 100
    S = 50;    // mean between 0 and 200
    R_nu_minus_one = 1.0/29;  // df of 30
    R_sigma = 1.0/50;         // sd of 50
}
generated quantities {
    real<lower = 0> nu_minus_one;
    real<lower = 1> nu;
    real mu;
    real<lower = 0> sigma;
    array[N] real y_sim;
    
    nu_minus_one = exponential_rng(R_nu_minus_one);
    nu = nu_minus_one + 1;
    mu = normal_rng(M, S);
    sigma = exponential_rng(R_sigma);
    
    for(n in 1:N) {
        y_sim[n] = student_t_rng(nu, mu, sigma);
    }
}
```

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

```{r}
samples_IQ_robust_prior <- tidy_draws(fit_IQ_robust_prior)
samples_IQ_robust_prior
```

### Examine prior

```{r}
mcmc_hist(samples_IQ_robust_prior,
          pars = "nu")
```

```{r}
mcmc_hist(samples_IQ_robust_prior,
          pars = "mu")
```

```{r}
mcmc_hist(samples_IQ_robust_prior,
          pars = "sigma")
```


```{r}
mcmc_pairs(fit_IQ_robust_prior,
           pars = c("nu", "mu", "sigma"))
```

### Prior predictive distribution

```{r}
y_sim_IQ_robust_prior <- samples_IQ_robust_prior %>%
    dplyr::select(starts_with("y_sim")) %>%
    as.matrix()
```

```{r}
ppd_hist(ypred = y_sim_IQ_robust_prior[1:20,])
```

```{r}
ppd_boxplot(ypred = y_sim_IQ_robust_prior[1:10, ],
            notch = FALSE)
```

As before, these priors are likely too wide, especially now that the Student t model is capable of generating outliers.

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


## Model

In the model code below, we list `mu` and `sigma` as the parameters of interest. The `a` and `b` parameters are defined in the `transformed parameters` block. They make the model easier to read, but they are not the main objects of our inference.

```{stan, output.var = "IQ_robust", cache = TRUE}
data {
    int<lower = 0> N;
    array[N] real y;
}
transformed data {
    real M;
    real<lower = 0> S;
    real<lower = 0> R_nu_minus_one;
    real<lower = 0> R_sigma;
    
    M = 100;   // mean centered at 100
    S = 50;    // mean between 0 and 200
    R_nu_minus_one = 1.0/29;  // df of 30
    R_sigma = 1.0/50;         // sd of 50
}
parameters {
    real<lower = 0> nu_minus_one;
    real mu;
    real<lower = 0> sigma;
}
transformed parameters {
    real<lower = 1> nu;
    
    nu = nu_minus_one + 1;
}
model {
    nu_minus_one ~ exponential(R_nu_minus_one);
    mu ~ normal(M, S);
    sigma ~ exponential(R_sigma);
    y ~ student_t(nu, mu, sigma);
}
generated quantities {
    array[N] real y_sim;
    
    for (n in 1:N) {
        y_sim[n] = student_t_rng(nu, mu, sigma);
    }
}
```

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

```{r}
samples_IQ_robust <- tidy_draws(fit_IQ_robust)
samples_IQ_robust
```

## Model diagnostics

```{r}
stan_trace(fit_IQ_robust,
           pars = c("nu", "mu", "sigma"))
```

```{r}
mcmc_acf(fit_IQ_robust,
         pars = c("nu", "mu", "sigma"))
```

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

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


## Model summary

```{r}
print(fit_IQ_robust,
      pars = c("nu", "mu", "sigma"))
```


## Model visualization

```{r}
mcmc_areas(fit_IQ_robust, pars = "nu")
```

```{r}
mcmc_areas(fit_IQ_robust, pars = "mu")
```

```{r}
mcmc_areas(fit_IQ_robust, pars = "sigma")
```

```{r}
pairs(fit_IQ_robust,
      pars = c("nu", "mu", "sigma"))
```


## Posterior predictive check

```{r}
y_sim_IQ_robust <- samples_IQ_robust %>%
    dplyr::select(starts_with("y_sim")) %>%
    as.matrix()
```

```{r}
ppc_hist(y = y,
         yrep = y_sim_IQ_robust[1:19, ])
```

```{r}
ppc_boxplot(y = y,
            yrep = y_sim_IQ_robust[1:10, ],
            notch = FALSE)
```

```{r}
ppc_dens_overlay(y = y,
                 yrep = y_sim_IQ_robust[1:50, ])
```

```{r}
ppc_scatter(y = y,
            yrep = y_sim_IQ_robust[1:9, ])
```

```{r}
ppc_scatter_avg(y = y,
                yrep = y_sim_IQ_robust)
```

```{r}
ppc_intervals(y = y,
              yrep = y_sim_IQ_robust[1:19, ])
```

```{r}
ppc_stat(y = y,
         yrep = y_sim_IQ_robust)
```

```{r}
ppc_stat_2d(y = y,
            yrep = y_sim_IQ_robust)
```

There are a few simulated data sets that are outliers, including one huge one. It doesn't invalidate our inference; our estimates of the parameters of interest are stable and reasonable. If there were more outliers like this, though, we would want to reconsider our model to make it less possible to get extreme outliers like this.