Ideal point model in stan

idealpoint.R

# Ideal point model in stan / R
# example taken & simplified from https://medewitt.github.io/resources/stan_ideal_point.html
library(tidyverse)
library(cmdstanr)

# simulate data: 100 legislators, 150 votes
set.seed(1834)
N_legislators <- 50
N_bills <- 150

vote_probs <- c(
  rep(0.9, 21), # Liberals (majority Gov. party)
  rep(0.7, 11), # Conservatives (uneasy coalition)
  rep(0.3,  8), # Socialists (opposition)
  rep(0.25, 5), # Greens (opposition)
  rep(0.01, 3), # Religious (opposition)
  rep(0.5,  2)  # Independents (random)
)

votes_data <- tibble(
  bill_id = rep(1:N_bills, each = N_legislators),
  legislator_id = rep(1:N_legislators, N_bills),
  vote_prob = rep(vote_probs, N_bills),
  vote = rbinom(N_bills*N_legislators, 1, vote_prob),
  legislator_party = as_factor(case_when(
    legislator_id <= 20 ~ "The Classic Liberal Party",
    legislator_id > 20 & legislator_id <= 32 ~ "The Conservative Party",
    legislator_id > 32 & legislator_id <= 40 ~ "The Socialist Party",
    legislator_id > 40 & legislator_id <= 45 ~ "The Green Party",
    legislator_id > 45 & legislator_id <= 48 ~ "The Religious Party",
    TRUE ~ "Independent"
  ))
)

# transform to stan data
stan_data <- list(
  N_legislators = N_legislators,
  N_bills = N_bills,
  N_observations = nrow(votes_data),
  idx_legi = votes_data$legislator_id,
  idx_bill = votes_data$bill_id,
  y = votes_data$vote
)

# compile and fit the model!
mod <- cmdstan_model("idealpoint.stan")
fit <- mod$optimize(data = stan_data, jacobian = TRUE)

# make a plot on how well we did!
votes_data |>
  filter(bill_id == 1) |>
  mutate(alpha = res$mle("alpha")) |>
  ggplot(aes(x = legislator_id, y = alpha, color = legislator_party)) +
  geom_point() +
  geom_point(aes(y = -log(vote_prob) + mean(log(vote_prob))), color = "black", pch = 2) +
  theme_linedraw() +
  labs(
    title = "Estimation accuracy of legislator parameter",
    x = "Legislator ID",
    color = "Legislator party",
    y = "Log-odds ratio",
    subtitle = "Inaccuracies due to regularization towards the mean (prior & scaling)"
  )

ggsave("idealpoint.png", width = 9, height = 5)

idealpoint.stan

data {
  int<lower=1> N_legislators;
  int<lower=1> N_bills;
  int<lower=1> N_observations;
  array[N_observations] int<lower=1, upper=N_legislators> idx_legi; // Legislator for observation n
  array[N_observations] int<lower=1, upper=N_bills> idx_bill; // Bill for observation n
  array[N_observations] int<lower=0, upper=1> y; // Vote of observation n
}
parameters {
  vector[N_bills] gamma;
  vector[N_legislators] alpha;
  vector[N_bills] beta;
}
model {
  //priors on parameters
  alpha ~ normal(0, 1);
  beta  ~ normal(0, 10);
  gamma ~ normal(0, 10);

  // likelihood
  y ~ bernoulli_logit(gamma[idx_bill] .* alpha[idx_legi] - beta[idx_bill]);
}