Assignment: Conjoint Analysis

library(tidyverse) # ggplot, dplyr, and friends
library(haven) # Read Stata files
library(cmdstanr)
library(survey) # Panel-ish regression models
library(brms) # The best formula-based interface to Stan
library(broom) # Convert model objects to tidy data frames
library(marginaleffects)
library(scales) # Nicer labeling functions
library(broom.helpers) # Add empty reference categories to tidy model data frames
library(ggforce) # For facet_col()
library(tidybayes) # Manipulate Stan results in tidy ways


candidate <- read_stata("candidate.dta") %>%
  as_factor()

variable_lookup <- tibble::tribble(
  ~variable,    ~variable_nice,
  "atmilitary", "Military",
  "atreligion", "Religion",
  "ated",       "Education",
  "atprof",     "Profession",
  "atmale",     "Gender",
  "atinc",      "Income",
  "atrace",     "Race",
  "atage",      "Age"
) %>%
  mutate(variable_nice = fct_inorder(variable_nice))


candidate_svy_design <- survey::svydesign(
  ids = ~resID,
  weights = ~1,
  data = candidate
)

model_svy <- survey::svyglm(
  selected ~ atmilitary + atreligion + ated +
    atprof + atinc + atrace + atage + atmale,
  design = candidate_svy_design
)

# Extract all the right-hand variables
rhs_vars <- all.vars(stats::update(formula(model_svy), 0 ~ .))

# Make a data frame of all the levels of all the non-numeric things
model_variable_levels <- tibble(
  variable = rhs_vars
) %>%
  mutate(levels = map(variable, ~ {
    x <- candidate[[.x]]
    if (is.numeric(x)) {
      ""
    } else if (is.factor(x)) {
      levels(x)
    } else {
      sort(unique(x))
    }
  })) %>%
  unnest(levels) %>%
  mutate(term = paste0(variable, levels))


# BAYESIAN MODELING!
priors <- c(
  prior(normal(0, 1), class = Intercept),
  prior(normal(0, 1), class = b),
  prior(exponential(1), class = sd)
)

model_brms <- brm(
  bf(selected ~ atmilitary + atreligion + ated +
    atprof + atinc + atrace + atage + atmale +
    (1 | resID)),
  data = candidate,
  family = bernoulli(link = "logit"),
  prior = priors,
  chains = 4,
  cores = detectCores() - 1,
  iter = 2000,
  seed = 1234,
  backend = "cmdstanr", threads = threading(2), refresh = 0,
  file = "candidate_model_brms"
)

# Define the horseshoe prior for the fixed effects
hs_prior <- c(
  prior(horseshoe(df = 1), class = b), # horseshoe prior for regression coefficients
  prior(normal(0, 1), class = Intercept), # prior for the intercept
  prior(exponential(1), class = sd) # prior for random effects standard deviation
)

# Update the brm model with the horseshoe prior
model_brms_hs <- brm(
  bf(selected ~ atmilitary + atreligion + ated +
    atprof + atinc + atrace + atage + atmale +
    (1 | resID)),
  data = candidate,
  family = bernoulli(link = "logit"),
  prior = hs_prior,
  chains = 4,
  cores = detectCores() - 1,
  iter = 3000,
  seed = 1234,
  backend = "cmdstanr",
  threads = threading(2),
  refresh = 0,
  file = "candidate_model_brms_hs",
  control = list(adapt_delta = 0.99, max_treedepth = 15)
)

# Load necessary libraries
library(dplyr)
library(ggplot2)
library(broom)

# Extract coefficients
coef_brms <- as.data.frame(fixef(model_brms)) %>%
  rownames_to_column("variable") %>%
  mutate(model = "Standard")

coef_brms_hs <- as.data.frame(fixef(model_brms_hs)) %>%
  rownames_to_column("variable") %>%
  mutate(model = "Horseshoe")

# Combine coefficients into a single data frame
coef_compare <- bind_rows(coef_brms, coef_brms_hs) %>%
  select(variable, Estimate, Est.Error, model)

# Plot the comparison
ggplot(coef_compare, aes(x = variable, y = Estimate, color = model)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(aes(ymin = Estimate - Est.Error, ymax = Estimate + Est.Error),
    position = position_dodge(width = 0.5), width = 0.2
  ) +
  labs(
    title = "Coefficient Comparison: Standard vs. Horseshoe Priors",
    x = "Variables", y = "Estimate"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_blank()) # Remove x-axis text labels