breakfast analysis

breakfast.R

library(tidyverse)
library(rstan)
library(cmdstanr)
dat <- read.delim("data_attractiveness_juge_results.csv")

### Set subject IDs using the same factor levels
subject_levels <- unique(c(dat$ID_left,
                           dat$ID_right))

dat <- dat |>
    mutate(ID_left = factor(ID_left,
                            levels = subject_levels),
           ID_right = factor(ID_right,
                             levels = subject_levels),
           ID_left.num = as.numeric(ID_left),
           ID_right.num = as.numeric(ID_right))

###

dat <- dat |>
    mutate(ID_juge.num = as.numeric(factor(ID_juge)))

## Set up for stan
stan_data <- list(y = dat$score_attractiveness_left,
                  x = dat$diff_breakfast_type,
                  id_l = dat$ID_left.num,
                  id_r = dat$ID_right.num,
                  id_j = dat$ID_juge.num,
                  n_subjects = length(subject_levels),
                  n_judges = max(dat$ID_juge.num),
                  n_obs = nrow(dat))
                  
stan_code <- '

data {
 int<lower=0> n_subjects;
 int<lower=0> n_judges;
 int<lower=0> n_obs;
 array[n_obs] int<lower=1, upper=n_subjects> id_l;
 array[n_obs] int<lower=1, upper=n_subjects> id_r;
 array[n_obs] int<lower=1, upper=n_judges> id_j;
 array[n_obs] int<lower=0, upper=1> y;
 array[n_obs] int<lower=-1, upper=1> x;
}
parameters {
 vector [n_subjects] theta;
 vector [n_judges] kappa;
 real beta;
 real alpha;
}
model {
  theta ~ std_normal();
  kappa ~ std_normal();
  for (i in 1:n_obs){
    y[i] ~ bernoulli_logit(alpha + beta * x[i] + theta[id_l[i]] - theta[id_r[i]] + kappa[id_j[i]]);
  }
}

'

writeLines(stan_code, con = "breakfast.stan")
model <- cmdstanr::cmdstan_model(stan_file = 'breakfast.stan')

fit <- model$sample(data = stan_data,
                    parallel_chains = 4,
                    chains = 4,
                    seed = 2408)

s <- fit$summary()
s |> filter(variable == "beta")
# A tibble: 1 × 10
##   variable   mean median    sd   mad     q5    q95  rhat ess_bulk ess_tail
##   <chr>     <num>  <num> <num> <num>  <num>  <num> <num>    <num>    <num>
## 1 beta     -0.489 -0.487 0.200 0.198 -0.816 -0.166  1.00     897.    1479.