storable_tv.Rmd

---
title: "Regression and Other Stories: Storable"
author: "Andrew Gelman, Jennifer Hill, Aki Vehtari"
date: "`r Sys.Date()`"
output:
  github_document:
    toc: true
---
Tidyverse version by Bill Behrman.

Ordered categorical data analysis with a study from experimental
economics, on the topic of "storable votes". See Chapter 15 in
Regression and Other Stories.

-------------

```{r, message=FALSE}
# Packages
library(tidyverse)
library(rstanarm)

# Parameters
  # Data directory
dir_data <- here::here("Storable/data")
  # People to model
people <- 
  tribble(
    ~person, ~behavior,
        101, "Perfectly monotonic",
        303, "One fuzzy and one sharp cutpoint",
        409, "Monotonic with one outlier",
        405, "Only 1's and 3's",
        504, "Almost only 3's",
        112, "Erratic"
  )
  # Common code
file_common <- here::here("_common.R")

#===============================================================================

# Run common code
source(file_common)
```

# 15 Other generalized linear models

## 15.5 Ordered and unordered categorical regression

### Example of ordered categorical regression

Data

```{r, message=FALSE}
data <- 
  dir_data %>% 
  fs::dir_ls(regexp = "/\\dplayergames.csv") %>% 
  map_dfr(read_csv) %>% 
  select(person, value, vote) %>% 
  mutate(vote_ord = ordered(vote))

summary(data)
```

EDA

```{r}
data %>% 
  count(person, name = "n_rows") %>% 
  count(n_rows)
```

There are 87 people in the data, 76 have 30 observations each, 1 has 29, and 10 have 20.

Votes cast by value.

```{r, message=FALSE}
data %>% 
  ggplot(aes(value)) +
  stat_ydensity(
    aes(y = vote_ord),
    draw_quantiles = c(0.25, 0.5, 0.75),
    scale = "count"
  ) +
  geom_smooth(aes(y = vote)) +
  labs(
    title = "Votes cast by value",
    x = "Value",
    y = "Votes cast"
  )
```

Low values are associated with one vote, middle values are associated with two votes, and high values are associated with three votes.

#### Fitting the model in R

Fit model for person 401.

```{r}
set.seed(768)

fit_401 <- 
  stan_polr(
    vote_ord ~ value,
    data = data %>% filter(person == 401),
    refresh = 0,
    prior = R2(location = 0.3, what = "mean")
  )

print(fit_401, digits = 2)
```

Fit model to people in `people`.

```{r}
set.seed(768)

fit_people <- 
  data %>% 
  filter(person %in% people$person) %>% 
  nest(data = !person) %>% 
  rowwise() %>% 
  mutate(
    fit = 
      list(
        stan_polr(
          vote_ord ~ value,
          data = data,
          refresh = 0,
          prior = R2(location = 0.3, what = "mean"),
          adapt_delta = 0.9999
        )
      )
  )
```

#### Display the fitted model

Calculate cutpoints and sigma for each person.

```{r}
data_people <- 
  fit_people %>% 
  mutate(
    sims = list(as_tibble(fit)),
    c_1.5 = median(sims$`1|2` / sims$value),
    c_2.5 =  median(sims$`2|3` / sims$value),
    sigma = median(1 / sims$value)
  ) %>% 
  ungroup() %>% 
  select(!c(fit, sims)) %>% 
  mutate(person = factor(person, levels = people$person))
```

Calculate expected response curves.

```{r}
expected <- function(x, c_1.5, c_2.5, sigma) {
  p_1.5 <- plogis((x - c_1.5) / sigma)
  p_2.5 <- plogis((x - c_2.5) / sigma)
  1 * (1 - p_1.5) + 2 * (p_1.5 - p_2.5) + 3 * p_2.5
}

curves <- 
  data_people %>% 
  select(!data) %>% 
  mutate(
    x = list(seq_range(range(data$value))),
    y = pmap(list(x, c_1.5, c_2.5, sigma), expected)
  ) %>% 
  select(person, x, y) %>% 
  unnest(cols = c(x, y))
```

Cutpoint segments.

```{r}
segments <- 
  data_people %>% 
  select(person, c_1.5, c_2.5) %>% 
  pivot_longer(cols = !person, names_to = "cutpoint", values_to = "x") %>% 
  mutate(
    y = str_extract(cutpoint, pattern = "[12]") %>% as.double(),
    yend = y + 1
  ) %>% 
  filter(x >= 1, x <= 100)
```

Example individuals from storable votes study.

```{r, fig.asp=1.2}
person_label <- function(x) {
  people %>% 
    filter(person == x) %>% 
    pull(behavior)
}

data_people %>% 
  unnest(data) %>% 
  ggplot() +
  geom_segment(
    aes(x = x, xend = x, y = y, yend = yend),
    data = segments,
    color = "grey60"
  ) +
  geom_line(aes(x, y), data = curves) +
  geom_count(aes(value, vote)) +
  facet_wrap(
    facets = vars(person),
    ncol = 2,
    labeller = labeller(person = person_label)
  ) +
  scale_y_continuous(breaks = 1:3, minor_breaks = NULL) +
  scale_size(range = c(1, 3), guide = "none") +
  labs(
    title = "Example individuals from storable votes study",
    x = "Value",
    y = "Votes cast"
  )
```