Trolley Data

Data Set

Rate the moral appropriateness of the action in the story. Response: 1 (Forbidden) to 7 (Obligatory)

Setup

library(tidyverse)
library(tidybayes)
library(bayesplot)
library(cmdstanr)
cmdstanr::register_knitr_engine()

library(dagitty)
library(ggdag)

df <- read.csv("new_trolley.csv") |>
  mutate(
    response = as_factor(response),
    id = as_factor(id),
    case = as_factor(case),
    story = as_factor(story),
    gender = ifelse(male, "female", "male") |> as_factor(),
    age_std = scale(age) |> drop(),
    edu = as_factor(edu) |>
      fct_relevel(
        "Elementary School",
        "Middle School",
        "Some High School",
        "High School Graduate",
        "Some College",
        "Bachelor's Degree",
        "Master's Degree",
        "Graduate Degree"
      )
  )

Model 1: A, I, C

dagitty(
  'dag {
"A,I,C" [exposure,pos="-0.204,-0.022"]
Response [outcome,pos="-0.083,-0.022"]
"A,I,C" -> Response
}'
) |>
  tidy_dagitty() |>
  ggdag_status(text_size = 2) +
  theme_dag_gray(legend.position = "bottom")

\[ \begin{align} \mathit{R}_i &\sim \mathrm{OrderedLogit}(\phi_i, \alpha) \\ \phi_i &= \beta_{A}A_i + \beta_{I}I_i + \beta_{C}C_i \\ \beta &\sim \mathrm{Normal}(0, 0.5) \\ \alpha_j &\sim \mathrm{Normal}(0, 1) \end{align} \]

new_trolley.stan

data {
  int<lower=1> n;
  int<lower=1> n_response;
  
  array[n] int<lower=1, upper=n_response> response;
  vector<lower=0, upper=1>[n] contact;
  vector<lower=0, upper=1>[n] intention;
  vector<lower=0, upper=1>[n] action;
}
parameters {
  ordered[n_response - 1] alpha;
  real bA;
  real bI;
  real bC;
}
model {
  vector[n] phi = bA * action + bI * intention + bC * contact;
  
  alpha ~ normal(0, 1);
  bA ~ normal(0, 0.5);
  bI ~ normal(0, 0.5);
  bC ~ normal(0, 0.5);
  
  response ~ ordered_logistic(phi, alpha);
}

Fitting and diagnostics

model <- cmdstan_model("new_trolley.stan")
fit <- compose_data(df) |> model$sample(chains = 4, parallel_chains = 4)

mcmc_trace(fit$draws(), regex_pars = "alpha")

mcmc_trace(fit$draws(), regex_pars = "b")

mcmc_pairs(fit$draws(), regex_pars = "alpha")

mcmc_pairs(fit$draws(), regex_pars = "b")

Interpretation

posterior <- fit$draws() |>
  recover_types(df) |>
  spread_draws(alpha[n], bA, bI, bC) |>
  pivot_wider(names_from = n, names_prefix = "alpha_", values_from = alpha)

prediction_grid <- expand_grid(
  action = c(0, 1),
  intention = c(0, 1),
  contact = c(0, 1)
)

predictions <- posterior |>
  crossing(prediction_grid) |>
  mutate(
    phi = bA * action + bI * intention + bC * contact,
    p1 = plogis(alpha_1 - phi),
    p2 = plogis(alpha_2 - phi) - plogis(alpha_1 - phi),
    p3 = plogis(alpha_3 - phi) - plogis(alpha_2 - phi),
    p4 = plogis(alpha_4 - phi) - plogis(alpha_3 - phi),
    p5 = plogis(alpha_5 - phi) - plogis(alpha_4 - phi),
    p6 = plogis(alpha_6 - phi) - plogis(alpha_5 - phi),
    p7 = 1 - plogis(alpha_6 - phi)
  ) |>
  select(-c(phi)) |>
  pivot_longer(
    starts_with("p"),
    names_to = "cat",
    values_to = "prob"
  ) |>
  mutate(category = readr::parse_number(cat))

prediction_summary <- predictions |>
  group_by(action, intention, contact, category) |>
  summarise(
    mean_prob = mean(prob),
    lower = quantile(prob, 0.1),
    upper = quantile(prob, 0.9),
    .groups = "drop"
  ) |>
  mutate(
    combo = paste0(action, ",", intention, ",", contact) |>
      as_factor() |>
      fct_relevel(
        "0,0,0",
        "1,0,0",
        "0,1,0",
        "0,0,1",
        "1,1,0",
        "1,0,1",
        "0,1,1",
        "1,1,1"
      )
  )

prediction_summary |>
  filter(action + intention + contact <= 1) |>
  ggplot(
    aes(x = factor(category), y = mean_prob, ymin = lower, ymax = upper),
  ) +
  geom_col(width = 0.4, fill = "lightblue", color = "royalblue") +
  geom_errorbar(width = 0.2, alpha = 0.5) +
  facet_wrap(~combo, nrow = 2) +
  labs(x = "Action, Intention, Contact", y = NULL)

prediction_summary |>
  filter(action + intention + contact > 1) |>
  ggplot(
    aes(x = factor(category), y = mean_prob, ymin = lower, ymax = upper),
  ) +
  geom_col(width = 0.4, fill = "lightblue", color = "royalblue") +
  geom_errorbar(width = 0.2, alpha = 0.5) +
  facet_wrap(~combo, nrow = 2) +
  labs(x = "Action, Intention, Contact", y = NULL)

Model 2: + Gender

new_trolley2.stan

data {
  int<lower=1> n;
  
  int<lower=1> n_response;
  int<lower=1> n_gender;
  
  array[n] int<lower=1, upper=n_response> response;
  array[n] int<lower=0, upper=1> contact;
  array[n] int<lower=0, upper=1> intention;
  array[n] int<lower=0, upper=1> action;
  
  array[n] int<lower=1, upper=n_gender> gender;
}
parameters {
  vector[n_gender] bA;
  vector[n_gender] bI;
  vector[n_gender] bC;
  
  ordered[n_response - 1] alpha;
}
model {
  // linear predictor
  vector[n] phi;
  for (i in 1 : n) 
    phi[i] = bA[gender[i]] * action[i] + bI[gender[i]] * intention[i]
             + bC[gender[i]] * contact[i];
  
  // priors
  bA[n_gender] ~ normal(0, 0.5);
  bI[n_gender] ~ normal(0, 0.5);
  bC[n_gender] ~ normal(0, 0.5);
  
  alpha ~ normal(0, 1);
  
  // likelihood function
  response ~ ordered_logistic(phi, alpha);
}

model2 <- cmdstan_model("new_trolley2.stan")
fit2 <- compose_data(df) |>
  model2$sample(chains = 4, parallel_chains = 4)

fit2$summary()

# A tibble: 13 × 10
   variable       mean     median     sd    mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>      <dbl>  <dbl>  <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__     -18486.    -18486.    2.46   2.35   -1.85e+4 -1.85e+4 1.00     1792.
 2 bA[1]        -0.893     -0.893 0.0512 0.0522 -9.78e-1 -8.09e-1 1.000    4304.
 3 bA[2]        -0.535     -0.535 0.0484 0.0477 -6.13e-1 -4.56e-1 1.00     3790.
 4 bI[1]        -0.900     -0.901 0.0487 0.0491 -9.80e-1 -8.18e-1 1.000    5510.
 5 bI[2]        -0.560     -0.560 0.0438 0.0448 -6.32e-1 -4.88e-1 1.00     5446.
 6 bC[1]        -1.08      -1.08  0.0678 0.0691 -1.19e+0 -9.70e-1 1.00     4678.
 7 bC[2]        -0.843     -0.843 0.0631 0.0619 -9.49e-1 -7.39e-1 1.00     5637.
 8 alpha[1]     -2.84      -2.84  0.0471 0.0478 -2.92e+0 -2.76e+0 1.00     3070.
 9 alpha[2]     -2.16      -2.16  0.0426 0.0424 -2.23e+0 -2.09e+0 1.00     3170.
10 alpha[3]     -1.57      -1.57  0.0397 0.0401 -1.64e+0 -1.51e+0 1.00     3240.
11 alpha[4]     -0.541     -0.541 0.0368 0.0359 -6.02e-1 -4.79e-1 1.00     3494.
12 alpha[5]      0.136      0.136 0.0372 0.0372  7.51e-2  1.98e-1 1.00     3733.
13 alpha[6]      1.05       1.05  0.0406 0.0408  9.84e-1  1.12e+0 1.000    4220.
# ℹ 1 more variable: ess_tail <dbl>

Model 3: + Age and Education

new_trolley3.stan

data {
  int<lower=1> n;
  
  int<lower=1> n_response;
  int<lower=1> n_gender;
  int<lower=1> n_edu;
  
  int<lower=1> n_story;
  int<lower=1> n_id;
  
  array[n] int<lower=1, upper=n_response> response;
  array[n] int<lower=0, upper=1> contact;
  array[n] int<lower=0, upper=1> intention;
  array[n] int<lower=0, upper=1> action;
  
  array[n] int<lower=1, upper=n_gender> gender;
  array[n] int age;
  vector[n] age_std;
  array[n] int edu;
  
  array[n] int<lower=1, upper=n_story> story;
  array[n] int<lower=1, upper=n_id> id;
}
parameters {
  vector[n_gender] bA;
  vector[n_gender] bI;
  vector[n_gender] bC;
  
  real bAge;
  real bEdu;
  
  ordered[n_response - 1] alpha;
  simplex[n_edu] delta;
}
model {
  // linear predictor
  vector[n] phi;
  for (i in 1 : n) 
    phi[i] = bA[gender[i]] * action[i] + bI[gender[i]] * intention[i]
             + bC[gender[i]] * contact[i] + bAge * age_std[i]
             + bEdu * sum(delta[1 : edu[i]]);
  
  // priors
  bA[n_gender] ~ normal(0, 0.5);
  bI[n_gender] ~ normal(0, 0.5);
  bC[n_gender] ~ normal(0, 0.5);
  
  bAge ~ normal(0, 0.5);
  bEdu ~ normal(0, 0.5);
  delta ~ dirichlet(rep_vector(2, n_edu));
  
  alpha ~ normal(0, 1);
  
  // likelihood function
  response ~ ordered_logistic(phi, alpha);
}

model3 <- cmdstan_model("new_trolley3.stan")
fit3 <- compose_data(df) |>
  model3$sample(chains = 4, parallel_chains = 4)

options(tibble.print_max = 50)
fit3$summary()

# A tibble: 23 × 10
   variable        mean    median     sd    mad       q5      q95  rhat ess_bulk
   <chr>          <dbl>     <dbl>  <dbl>  <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__     -18514.      -1.85e+4 3.37   3.34   -1.85e+4 -1.85e+4 1.00     1761.
 2 bA[1]        -0.890   -8.90e-1 0.0519 0.0518 -9.76e-1 -8.04e-1 1.00     2945.
 3 bA[2]        -0.543   -5.42e-1 0.0491 0.0490 -6.25e-1 -4.62e-1 1.00     3180.
 4 bI[1]        -0.901   -9.02e-1 0.0494 0.0496 -9.80e-1 -8.19e-1 1.000    3755.
 5 bI[2]        -0.564   -5.63e-1 0.0443 0.0454 -6.37e-1 -4.93e-1 1.00     3394.
 6 bC[1]        -1.08    -1.08e+0 0.0693 0.0686 -1.19e+0 -9.65e-1 1.00     3499.
 7 bC[2]        -0.849   -8.49e-1 0.0618 0.0613 -9.50e-1 -7.47e-1 0.999    2817.
 8 bAge         -0.0842  -8.46e-2 0.0216 0.0220 -1.19e-1 -4.85e-2 1.00     2245.
 9 bEdu          0.197    2.14e-1 0.146  0.121  -7.76e-2  3.98e-1 1.00     1193.
10 alpha[1]     -2.70    -2.70e+0 0.118  0.101  -2.91e+0 -2.53e+0 1.00     1152.
11 alpha[2]     -2.02    -2.01e+0 0.117  0.0997 -2.22e+0 -1.84e+0 1.00     1141.
12 alpha[3]     -1.43    -1.42e+0 0.116  0.0976 -1.64e+0 -1.26e+0 1.00     1133.
13 alpha[4]     -0.397   -3.88e-1 0.117  0.0983 -6.02e-1 -2.22e-1 1.00     1137.
14 alpha[5]      0.281    2.89e-1 0.117  0.0976  7.75e-2  4.54e-1 1.00     1156.
15 alpha[6]      1.19     1.20e+0 0.118  0.0973  9.88e-1  1.37e+0 1.00     1206.
16 delta[1]      0.140    1.23e-1 0.0881 0.0821  2.81e-2  3.11e-1 1.00     3843.
17 delta[2]      0.108    9.49e-2 0.0711 0.0663  1.98e-2  2.44e-1 1.000    3976.
18 delta[3]      0.120    1.05e-1 0.0787 0.0735  2.10e-2  2.73e-1 1.00     4651.
19 delta[4]      0.0922   7.75e-2 0.0630 0.0538  1.82e-2  2.15e-1 1.000    3447.
20 delta[5]      0.0744   5.80e-2 0.0610 0.0434  1.26e-2  1.98e-1 1.00     2131.
21 delta[6]      0.284    2.90e-1 0.133  0.137   5.09e-2  4.97e-1 1.00     1526.
22 delta[7]      0.0819   6.89e-2 0.0593 0.0504  1.40e-2  2.00e-1 1.00     2796.
23 delta[8]      0.0984   8.47e-2 0.0657 0.0605  1.81e-2  2.25e-1 1.00     4032.
# ℹ 1 more variable: ess_tail <dbl>

options(tibble.print_max = NULL)