09.2_Therapeutic_touch.Rmd

---
title: "9.2: Therapeutic touch"
output: html_notebook
---


## Introduction

This is the "Therapeutic touch" example in Section 9.2 of Kruschke.


## 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 "TherapeuticTouchData.csv" is in your project directory.)

```{r}
tt_data <- read_csv("TherapeuticTouchData.csv")
tt_data
```

In order to use the same Stan code as before, we'll need to aggregate the data so that it shows the number of successes for each subject (instead of the raw successes and failures).

```{r}
tt_data_sum <- tt_data %>%
    group_by(s) %>%
    summarize(n = n(), y_agg = sum(y))
tt_data_sum
```


```{r}
S <- NROW(tt_data_sum) # number of subjects
N <- tt_data_sum$n     # number of trials for each subject
y <- tt_data_sum$y_agg # number of successes for each subject
```

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


## Prior

### Simulation code

This is exactly the same Stan code as in `09.2c_Multiple_coins_from_a_singla_mint.Rmd` except that the name of the model is changed.

The one exception is at the bottom of the `generated quantities` block, where we will explicitly compute a few differences between subjects (as done in the book example.)

```{stan, output.var = "tt_prior", cache = TRUE}
data {
    int<lower = 0> S;
    array[S] int<lower = 0> N;
    array[S] int<lower = 0> y;
}
transformed data {
    real<lower = 0> A_omega;
    real<lower = 0> B_omega;
    real<lower = 0> S_kappa;
    real<lower = 0> R_kappa;

    A_omega = 2;       // hyperprior parameters for omega
    B_omega = 2;
    S_kappa =  0.01;   // hyperprior parameters for kappa
    R_kappa =  0.01;
}
generated quantities {
    real<lower = 0, upper = 1> omega;
    real<lower = 0> kappa_minus_two;
    real<lower = 2> kappa;
    array[S] real<lower = 0, upper = 1> theta;
    real<lower = 0> a;
    real<lower = 0> b;
    array[S] int<lower = 0> y_sim;
    real theta_1_14;
    real theta_1_28;
    real theta_14_28;
    
    omega = beta_rng(A_omega, B_omega);  // mode
    kappa_minus_two = gamma_rng(S_kappa, R_kappa);
    kappa = kappa_minus_two + 2;         // concentration
    a = omega * (kappa - 2) + 1;
    b = (1 - omega) * (kappa - 2) + 1;
        
    for (s in 1:S) {
        theta[s] = beta_rng(a, b);  // probability of success
    }

    y_sim = binomial_rng(N, theta); // vectorized
    
    theta_1_14 = theta[1] - theta[14];  // Some comparisons
    theta_1_28 = theta[1] - theta[28];
    theta_14_28 = theta[14] - theta[28];
}
```

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

```{r}
samples_tt_prior <- tidy_draws(fit_tt_prior)
samples_tt_prior
```

We'll examine just some subjects as in the book example. Here are some proportions calculated for subjects 1, 14, and 28.

```{r}
samples_tt_prior <- samples_tt_prior %>%
        mutate(`y_sim_prop[1]` = `y_sim[1]`/ N[1],
               `y_sim_prop[14]` = `y_sim[14]`/ N[14],
               `y_sim_prop[28]` = `y_sim[28]`/ N[28])
samples_tt_prior
```


### Examine prior

```{r}
mcmc_hist(samples_tt_prior,
          pars = "omega")
```

```{r}
mcmc_hist(samples_tt_prior,
          pars = "kappa")
```

```{r}
mcmc_hist(samples_tt_prior,
          pars = c("theta[1]", "theta[14]", "theta[28]"))
```

```{r}
mcmc_hist(samples_tt_prior,
          pars = vars(starts_with("theta_")))
```

```{r}
mcmc_pairs(fit_tt_prior,
           pars = c("omega", "kappa"))
```

```{r}
mcmc_pairs(fit_tt_prior,
           pars = vars("omega", "theta[1]", "theta[14]", "theta[28]"))
```


```{r}
mcmc_pairs(fit_tt_prior,
           pars = vars("kappa", "theta[1]", "theta[14]", "theta[28]"))
```

### Prior predictive distribution

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

You can't tell from the x-axis below, but these are PPDs for subjects 1, 14, and 28.

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


## Model

Again, the only change made below is in the `generated quantities` block to calculate a few differences between subjects.

```{stan, output.var = "tt", cache = TRUE}
data {
    int<lower = 0> S;
    array[S] int<lower = 0> N;
    array[S] int<lower = 0> y;
}
transformed data {
    real<lower = 0> A_omega;  
    real<lower = 0> B_omega;
    real<lower = 0> S_kappa;
    real<lower = 0> R_kappa;

    A_omega = 2;       // hyperprior parameters for omega
    B_omega = 2;
    S_kappa =  0.01;   // hyperprior parameters for kappa
    R_kappa =  0.01;    
}
parameters {
    real<lower = 0, upper = 1> omega;
    real<lower = 0> kappa_minus_two;
    array[S] real<lower = 0, upper = 1> theta;
}
transformed parameters {
    real<lower = 2> kappa;
    real<lower = 0> a;
    real<lower = 0> b;
    
    kappa = kappa_minus_two + 2;     // concentration
    a = omega * (kappa - 2) + 1;
    b = (1 - omega) * (kappa - 2) + 1;
}
model {
    omega ~ beta(A_omega, B_omega);  // mode
    kappa_minus_two ~ gamma(S_kappa, R_kappa);
    theta ~ beta(a, b);              // prior prob of success
    y ~ binomial(N, theta);          // likelihood
}
generated quantities {
    array[S] int<lower = 0> y_sim;
    real theta_1_14;
    real theta_1_28;
    real theta_14_28;
    
    y_sim = binomial_rng(N, theta);

    theta_1_14 = theta[1] - theta[14];  // Some comparisons
    theta_1_28 = theta[1] - theta[28];
    theta_14_28 = theta[14] - theta[28];
}
```

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

The warning above means that we have some uncertainly about how well some of the parameters have been sampled. We could follow the advice above and run longer chains (increasing the `iter` argument to the `sampling` function, set to 2000 by default). We won't do that here, but know that it's an option.

```{r}
samples_tt <- tidy_draws(fit_tt)
samples_tt
```

Again, we'll compute proportions:

```{r}
samples_tt <- samples_tt %>%
    mutate(`y_sim_prop[1]` = `y_sim[1]`/ N[1],
           `y_sim_prop[14]` = `y_sim[14]`/ N[14],
           `y_sim_prop[28]` = `y_sim[28]`/ N[28])
samples_tt
```


## Model diagnostics

```{r}
stan_trace(fit_tt,
           pars = c("omega", "kappa",
                    "theta[1]", "theta[14]", "theta[28]"))
```

```{r}
mcmc_acf(fit_tt,
         pars = vars("omega", "kappa",
                     "theta[1]", "theta[14]", "theta[28]"))
```

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

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


## Model summary

```{r}
print(fit_tt,
      pars = c("omega", "kappa", "theta"))
```


## Model visualization

```{r}
mcmc_areas(fit_tt,
           pars = c("omega", "theta[1]", "theta[14]", "theta[28]"))
```

```{r}
mcmc_areas(fit_tt,
           pars = vars(starts_with("theta_")))
```

```{r}
mcmc_hist(fit_tt,
          pars = "kappa")
```

```{r}
pairs(fit_tt,
      pars = c("omega", "theta[1]", "theta[14]", "theta[28]"))
```


## Posterior predictive check

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

```{r}
ppc_intervals(y = c(y[1]/N[1], y[14]/N[14], y[28]/N[28]),
              yrep = y_sim_tt)
```

Subject 1:

```{r}
ppc_stat(y = y[1] / N[1],
         yrep = as.matrix(y_sim_tt[ , 1]))
```

Subject 14:

```{r}
ppc_stat(y = y[14] / N[14],
         yrep = as.matrix(y_sim_tt[ , 2]))
```


Subject 28:

```{r}
ppc_stat(y = y[28] / N[28],
         yrep = as.matrix(y_sim_tt[ , 3]))
```