Fairness Analyses

Testing Qwen 3 (22B-235B)

Setup

Load Dependencies

# https://github.com/jmgirard/rocker-bayes/releases/tag/v6.0

library(tidyverse)
library(brms)
library(rstan)
library(easystats)
library(tidybayes)
library(emmeans)
library(gt)
library(posterior)

silently <- function(x) {
  suppressWarnings(suppressMessages(x))
}

my_priors <- c(
  set_prior("normal(0, 3)", class = "b"),
  set_prior("student_t(3, 0, 2.5)", class = "Intercept"),
  set_prior("student_t(3, 0, 2.5)", class = "sd"),
  set_prior("lkj_corr_cholesky(2)", class = "L"),
  set_prior("normal(0, 2)", class = "Intercept", dpar = "phi"),
  set_prior("normal(0, 2)", class = "b", dpar = "phi")
)

Prep Data

# Sessions used as fewshot examples (excluded from evaluation)
source("load_data.R")
excluded_ids <- load_exemplars()

# Predictions (public, ../data) joined to human ratings (NDA, ../nda)
ablation_raw <- load_prediction_sheets(
  conditions = "full",
  models = "Qwen 3 (22B-235B)"
)

# Collect item score labels and predictions (one sheet per item, three seeds each)
dat_items <-
  map(
    .x = 1:10,
    .f = \(i) {
      ablation_raw[[i + 1]] |>
        filter(session %in% excluded_ids[[i + 1]] == FALSE) |>
        transmute(
          session,
          patient,
          item = sprintf("%02d", i),
          label = ground_truth,
          seed_1 = rating_0,
          seed_2 = rating_1,
          seed_3 = rating_2
        ) |>
        pivot_longer(
          starts_with("seed_"),
          names_to = "seed",
          names_prefix = "seed_",
          values_to = "pred"
        )
    }
  ) |>
  bind_rows() |>
  transmute(
    session,
    patient,
    item = factor(item),
    seed = factor(seed),
    pred,
    label
  )

# Indirect total scores: first non-missing prediction per item, summed per session
dat_sums <-
  dat_items |>
  summarize(
    .by = c(session, patient, item),
    pred = first(pred, na_rm = TRUE),
    label = first(label)
  ) |>
  summarize(
    .by = c(session, patient),
    pred = sum(pred),
    label = sum(label)
  )

# Participant characteristics from NDA (ndar_subject01)
demo <-
  load_subjects() |>
  transmute(
    patient,
    diagnosis = fct_collapse(diagnosis,
      Depression = c("MDD", "MDD w/ psychosis"),
      Bipolar = c("BP1 (manic)", "BP1 (hypomanic)", "BP1 (depressed)",
                  "BP1 (mixed)", "BP2 (depressed)", "BP2 (hypomanic)", "BP NOS"),
      Psychosis = c("SZ", "SZA (bp)", "SZA (dep)", "Psychosis NOS"),
      other_level = "Substance"
    ),
    diagnosis = fct_relevel(diagnosis, after = 0,
      "Depression", "Bipolar", "Psychosis", "Substance"
    ),
    age,
    sex = factor(sex, levels = c("Female", "Male")),
    race = fct_collapse(race,
      White = "White",
      Black = "Black or African American",
      other_level = "Other"
    ),
    race = fct_relevel(race, after = 0,
      "White", "Black", "Other"
    ),
    education = factor(education, ordered = TRUE),
    education0 = as.integer(education) - 1
  ) |>
  # Standardize person-level predictors within the analytic sample only
  filter(patient %in% unique(dat_sums$patient)) |>
  mutate(agez = standardize(age))

# Interviewer identity (RA1-RA4) per session, from the public sessions file
ready_sums <-
  dat_sums |>
  left_join(select(load_sessions(), session, interviewer), by = "session") |>
  left_join(demo, by = "patient") |>
  mutate(
    interviewer = factor(
      interviewer,
      levels = c("RA1", "RA2", "RA3", "RA4"),
      labels = c("Interviewer 1", "Interviewer 2", "Interviewer 3", "Interviewer 4")
    ),
    patient = factor(patient),
    labelz = standardize(label)
  )

Create Functions

Postprocess Model Results

# Extract bias, sensitivity, and precision deviations from an audit model.
# Categorical predictors get per-level relative deviations (vs. grand mean);
# numeric predictors get a single slope-style deviation per metric.
audit_fairness <- function(model, group_var, contrast_label = group_var) {

  if (is.numeric(model$data[[group_var]])) {
    # ==========================================
    # CONTINUOUS LOGIC (Overall Trend)
    # ==========================================
    draws <- as_draws_df(model)

    bias_col <- grep(paste0("^b_", group_var, "$"), names(draws), value = TRUE)
    sens_col <- grep(paste0("^b_labelz:", group_var, "$|^b_", group_var, ":labelz$"), names(draws), value = TRUE)
    prec_col <- grep(paste0("^b_phi_", group_var, "$"), names(draws), value = TRUE)

    process_param <- function(col_name, metric_name) {
      if (length(col_name) == 0) return(NULL)
      vals <- draws[[col_name[1]]]
      tibble(
        contrast = contrast_label,
        Metric = metric_name,
        Deviation = mean(vals),
        CI_Low = quantile(vals, 0.025),
        CI_High = quantile(vals, 0.975),
        pd = as.numeric(bayestestR::p_direction(vals)$pd),
        p_val = 2 * (1 - pd)
      )
    }

    bind_rows(
      process_param(bias_col, "Bias (Sum Score)"),
      process_param(sens_col, "Sensitivity (Slope)"),
      process_param(prec_col, "Precision (Log Phi)")
    ) |>
      mutate(contrast = factor(contrast))

  } else {
    # ==========================================
    # CATEGORICAL LOGIC (Per-Level)
    # ==========================================
    orig_levels <- levels(model$data[[group_var]])

    # Per-level relative deviations from the per-draw grand mean
    # (emmeans/tidybayes draws arrive grouped; .by requires ungrouped data)
    process_posterior <- function(post_draws, type_name) {
      post_draws |>
        ungroup() |>
        mutate(contrast = as.character(contrast)) |>
        mutate(grand_mean = mean(.value), .by = .draw) |>
        mutate(deviation = .value - grand_mean) |>
        summarize(
          .by = contrast,
          Metric = type_name,
          Deviation = mean(deviation),
          CI_Low = quantile(deviation, 0.025),
          CI_High = quantile(deviation, 0.975),
          pd = as.numeric(bayestestR::p_direction(deviation)$pd),
          p_val = 2 * (1 - pd)
        )
    }

    # Bias: expected sum score per level at average severity
    grid <- tibble(!!group_var := orig_levels, labelz = 0)
    draws_bias <- model |>
      epred_draws(newdata = grid, re_formula = NA) |>
      ungroup() |>
      transmute(.draw, contrast = .data[[group_var]], .value = .epred)

    # Sensitivity: severity slope per level
    draws_sens <- silently(emtrends(model, specs = group_var, var = "labelz")) |>
      gather_emmeans_draws() |>
      ungroup() |>
      rename(contrast = all_of(group_var))

    # Precision: log phi per level
    draws_prec <- silently(emmeans(model, specs = group_var, dpar = "phi")) |>
      gather_emmeans_draws() |>
      ungroup() |>
      rename(contrast = all_of(group_var))

    bind_rows(
      process_posterior(draws_bias, "Bias (Sum Score)"),
      process_posterior(draws_sens, "Sensitivity (Slope)"),
      process_posterior(draws_prec, "Precision (Log Phi)")
    ) |>
      mutate(contrast = factor(contrast, levels = orig_levels)) |>
      arrange(Metric, contrast)
  }
}

Create Results Summary Table

format_deviation_table <- function(deviations_df) {
  deviations_df |>
    mutate(
      Metric_Short = case_when(
        grepl("Bias", Metric) ~ "Bias",
        grepl("Sensitivity", Metric) ~ "Sens",
        grepl("Precision", Metric) ~ "Prec",
        TRUE ~ "Other"
      ),
      CI_Label = sprintf("[%.2f, %.2f]", CI_Low, CI_High)
    ) |>
    select(contrast, Metric_Short, Deviation, p_val, CI_Label) |>
    pivot_wider(
      names_from = Metric_Short,
      values_from = c(Deviation, p_val, CI_Label),
      names_glue = "{Metric_Short}_{.value}"
    ) |>
    select(
      contrast,
      starts_with("Bias"),
      starts_with("Sens"),
      starts_with("Prec")
    ) |>
    gt() |>
    tab_spanner(label = "Bias (Intercept)", columns = starts_with("Bias")) |>
    tab_spanner(label = "Sensitivity (Slope)", columns = starts_with("Sens")) |>
    tab_spanner(label = "Precision (Noise)", columns = starts_with("Prec")) |>
    fmt_number(columns = ends_with("Deviation"), decimals = 2) |>
    fmt_number(columns = ends_with("p_val"), decimals = 3) |>
    text_transform(
      locations = cells_body(columns = ends_with("p_val")),
      fn = function(x) {
        vals <- suppressWarnings(as.numeric(x))
        ifelse(is.na(vals), x,
               ifelse(vals < 0.001, "< .001", sprintf("%.3f", vals))
        )
      }
    ) |>
    cols_label(
      contrast = "Group",
      Bias_Deviation = "Dev.", Bias_p_val = "p-val", Bias_CI_Label = "95% CI",
      Sens_Deviation = "Dev.", Sens_p_val = "p-val", Sens_CI_Label = "95% CI",
      Prec_Deviation = "Dev.", Prec_p_val = "p-val", Prec_CI_Label = "95% CI"
    ) |>
    cols_align(align = "left", columns = contrast) |>
    tab_header(
      title = "Fairness Audit Results",
      subtitle = "Systematic Deviations by Metric"
    ) |>
    tab_style(
      style = cell_text(weight = "bold"),
      locations = list(
        cells_body(columns = "Bias_p_val", rows = Bias_p_val < 0.05),
        cells_body(columns = "Sens_p_val", rows = Sens_p_val < 0.05),
        cells_body(columns = "Prec_p_val", rows = Prec_p_val < 0.05)
      )
    ) |>
    tab_options(table.font.size = 14)
}

# Helper: format an Est + 95% CI summary as a gt table
gt_est_table <- function(df, group_col, group_label, est_label, title, subtitle = NULL) {
  df |>
    mutate(CI_Label = sprintf("[%.2f, %.2f]", CI_Low, CI_High)) |>
    select(all_of(group_col), Est, CI_Label, any_of("p_val")) |>
    gt() |>
    fmt_number(columns = Est, decimals = 2) |>
    cols_label(.list = setNames(list(group_label, est_label, "95% CI"),
                                c(group_col, "Est", "CI_Label"))) |>
    cols_align(align = "left", columns = all_of(group_col)) |>
    tab_header(title = title, subtitle = subtitle) |>
    tab_options(table.font.size = 14)
}

Run Analyses

All audit models share the same structure; only the group variable differs. Fits are cached to .rds files, so settings changes require deleting the corresponding file to take effect.

fit_audit <- function(group_term, file, seed = NA) {
  brm(
    bf(
      as.formula(paste0(
        "pred | trials(60) ~ labelz * ", group_term, " + (1 + labelz | patient)"
      )),
      as.formula(paste0("phi ~ ", group_term))
    ),
    data = ready_sums,
    family = beta_binomial(link = "logit", link_phi = "log"),
    prior = my_priors,
    init = 0.1,
    warmup = 3000,
    iter = 5000,
    cores = 4,
    chains = 4,
    refresh = 0,
    seed = seed,
    file = file.path(fits_dir, file),
    file_refit = "on_change",
    control = list(adapt_delta = 0.99, max_treedepth = 12),
    backend = "cmdstanr"
  )
}

fit_sex  <- fit_audit("sex", file = "sex")
fit_age  <- fit_audit("agez", file = "age")
fit_race <- fit_audit("race", file = "race")
fit_edu  <- fit_audit("education0", file = "edu", seed = 1)

Warning: Rows containing NAs were excluded from the model.

fit_diag <- fit_audit("diagnosis", file = "diag")
fit_intv <- fit_audit("interviewer", file = "intv")

Warning: Rows containing NAs were excluded from the model.

Convergence Diagnostics

check_convergence <- function(fit, name) {
  n_divergent <- nuts_params(fit) |>
    filter(Parameter == "divergent__") |>
    pull(Value) |>
    sum()

  summarise_draws(as_draws_df(fit), "rhat", "ess_bulk", "ess_tail") |>
    filter(!variable %in% c("lprior", "lp__"), !str_starts(variable, "r_")) |>
    mutate(
      class = if_else(
        str_starts(variable, "sd_") | str_starts(variable, "cor_"),
        "Random", "Fixed"
      ),
      model = name
    ) |>
    summarize(
      .by = c(model, class),
      max_rhat = max(rhat, na.rm = TRUE),
      min_ess_bulk = min(ess_bulk, na.rm = TRUE),
      min_ess_tail = min(ess_tail, na.rm = TRUE),
      n_params = n(),
      n_divergent = n_divergent
    )
}

bind_rows(
  check_convergence(fit_sex, "Sex"),
  check_convergence(fit_age, "Age"),
  check_convergence(fit_race, "Race"),
  check_convergence(fit_edu, "Education"),
  check_convergence(fit_diag, "Diagnosis"),
  check_convergence(fit_intv, "Interviewer")
) |>
  arrange(class, model) |>
  gt(groupname_col = "class") |>
  fmt_number(columns = max_rhat, decimals = 3) |>
  fmt_number(columns = starts_with("min_ess"), decimals = 0) |>
  cols_label(
    model = "Model",
    max_rhat = "Max R-hat",
    min_ess_bulk = "Min Bulk ESS",
    min_ess_tail = "Min Tail ESS",
    n_params = "Parameters",
    n_divergent = "Divergences"
  ) |>
  cols_align(align = "left", columns = model) |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = list(
      cells_body(columns = max_rhat, rows = max_rhat >= 1.01),
      cells_body(columns = min_ess_bulk, rows = min_ess_bulk < 400),
      cells_body(columns = min_ess_tail, rows = min_ess_tail < 400),
      cells_body(columns = n_divergent, rows = n_divergent > 0)
    )
  ) |>
  tab_header(
    title = "Convergence Diagnostics by Model",
    subtitle = "Worst-case values across parameters (patient-level effects excluded)"
  ) |>
  tab_source_note("Bold marks values to inspect: R-hat >= 1.01, ESS < 400, or any divergent transitions.") |>
  tab_options(table.font.size = 14)

Convergence Diagnostics by Model
Worst-case values across parameters (patient-level effects excluded)
Model	Max R-hat	Min Bulk ESS	Min Tail ESS	Parameters	Divergences
Fixed
Age	1.002	1,823	3,351	8	0
Diagnosis	1.002	2,012	3,374	14	0
Education	1.002	2,137	3,269	8	0
Interviewer	1.002	1,714	2,575	14	0
Race	1.002	1,773	3,293	11	0
Sex	1.003	2,491	3,957	8	0
Random
Age	1.003	710	754	3	0
Diagnosis	1.006	650	752	3	0
Education	1.009	704	1,092	3	0
Interviewer	1.003	685	898	3	0
Race	1.006	702	928	3	0
Sex	1.009	650	910	3	0
Bold marks values to inspect: R-hat >= 1.01, ESS < 400, or any divergent transitions.

Collect and Format Results

res_deviations <-
  bind_rows(
    audit_fairness(fit_sex, group_var = "sex"),
    audit_fairness(fit_age, group_var = "agez", contrast_label = "Age (per SD)"),
    audit_fairness(fit_race, group_var = "race"),
    audit_fairness(fit_edu, group_var = "education0", contrast_label = "Education (per level)"),
    audit_fairness(fit_diag, group_var = "diagnosis"),
    audit_fairness(fit_intv, group_var = "interviewer")
  )

# The "Substance" factor level is retained in the fitted models (renaming it
# would invalidate the cached fits); clinical record review showed the group
# is heterogeneous, so it is reported as "Other". Displayed here as "Other
# Diagnosis" because the racial group also has an "Other" level and duplicate
# contrast names break the pivot in format_deviation_table()
res_deviations <- res_deviations |>
  mutate(contrast = fct_recode(contrast, "Other Diagnosis" = "Substance"))

format_deviation_table(res_deviations)

Fairness Audit Results
Systematic Deviations by Metric
Group	Bias (Intercept)			Sensitivity (Slope)			Precision (Noise)
	Dev.	p-val	95% CI	Dev.	p-val	95% CI	Dev.	p-val	95% CI
Female	0.30	0.260	[-0.22, 0.82]	−0.04	0.087	[-0.09, 0.01]	−0.33	0.103	[-0.80, 0.07]
Male	−0.30	0.260	[-0.82, 0.22]	0.04	0.087	[-0.01, 0.09]	0.33	0.103	[-0.07, 0.80]
Age (per SD)	−0.04	0.067	[-0.08, 0.00]	0.02	0.485	[-0.03, 0.07]	−0.31	0.169	[-0.82, 0.12]
White	0.30	0.472	[-0.53, 1.09]	−0.01	0.891	[-0.08, 0.07]	0.14	0.604	[-0.56, 0.74]
Black	0.20	0.733	[-0.95, 1.36]	0.04	0.499	[-0.07, 0.15]	0.01	0.929	[-0.81, 1.13]
Other	−0.50	0.374	[-1.59, 0.62]	−0.03	0.537	[-0.13, 0.07]	−0.15	0.667	[-0.95, 0.75]
Education (per level)	−0.07	0.001	[-0.11, -0.03]	0.00	0.925	[-0.05, 0.04]	−0.42	0.030	[-0.90, -0.04]
Depression	−0.14	0.871	[-1.69, 1.34]	−0.17	0.010	[-0.31, -0.04]	0.22	0.518	[-0.54, 0.96]
Bipolar	−1.46	0.052	[-3.00, 0.01]	−0.11	0.101	[-0.25, 0.02]	−0.11	0.752	[-0.90, 0.67]
Psychosis	−1.62	0.033	[-3.19, -0.15]	−0.08	0.215	[-0.22, 0.05]	0.51	0.264	[-0.38, 1.72]
Other Diagnosis	3.21	0.099	[-0.65, 7.28]	0.37	0.034	[0.03, 0.72]	−0.62	0.273	[-1.66, 0.93]
Interviewer 1	−0.43	0.328	[-1.29, 0.45]	0.01	0.803	[-0.07, 0.10]	−0.37	0.480	[-1.49, 0.72]
Interviewer 2	0.02	0.967	[-0.84, 0.89]	−0.07	0.095	[-0.16, 0.01]	−1.44	0.001	[-2.70, -0.48]
Interviewer 3	0.91	0.186	[-0.44, 2.32]	0.05	0.510	[-0.10, 0.19]	0.60	0.607	[-1.11, 2.94]
Interviewer 4	−0.50	0.457	[-1.81, 0.84]	0.01	0.847	[-0.12, 0.15]	1.21	0.214	[-0.59, 3.46]

Expected Predictions (Raw Scale)

Audit Group Sizes

ready_sums |>
  distinct(patient, diagnosis) |>
  count(diagnosis, name = "Patients") |>
  mutate(diagnosis = fct_recode(diagnosis, Other = "Substance")) |>
  gt() |>
  cols_label(diagnosis = "Audit Group") |>
  cols_align(align = "left", columns = diagnosis) |>
  tab_header(title = "Patients per Diagnostic Audit Group") |>
  tab_options(table.font.size = 14)

Patients per Diagnostic Audit Group
Audit Group	Patients
Depression	96
Bipolar	85
Psychosis	83
Other	13

ready_sums |>
  count(interviewer, name = "Sessions") |>
  gt() |>
  cols_label(interviewer = "Interviewer") |>
  cols_align(align = "left", columns = interviewer) |>
  tab_header(
    title = "Sessions per Interviewer",
    subtitle = "NA = interviewer identity unresolved (excluded from interviewer model)"
  ) |>
  tab_options(table.font.size = 14)

Sessions per Interviewer
NA = interviewer identity unresolved (excluded from interviewer model)
Interviewer	Sessions
Interviewer 1	206
Interviewer 2	243
Interviewer 3	39
Interviewer 4	40
NA	13

Sensitivity by Diagnostic Group

grid_diag <- expand_grid(
  diagnosis = levels(ready_sums$diagnosis),
  labelz = c(-1, 0, 1)
)

diag_display <- c("Depression", "Bipolar", "Psychosis", "Other")

slope_draws <- fit_diag |>
  epred_draws(newdata = grid_diag, re_formula = NA) |>
  ungroup() |>
  select(diagnosis, labelz, .draw, .epred) |>
  mutate(diagnosis = recode(diagnosis, Substance = "Other")) |>
  pivot_wider(names_from = labelz, values_from = .epred, names_prefix = "SD_") |>
  mutate(slope_points = SD_1 - SD_0)

slope_draws |>
  mutate(diagnosis = factor(diagnosis, levels = diag_display)) |>
  summarize(
    .by = diagnosis,
    Est = mean(slope_points),
    CI_Low = quantile(slope_points, 0.025),
    CI_High = quantile(slope_points, 0.975)
  ) |>
  arrange(diagnosis) |>
  gt_est_table(
    group_col = "diagnosis", group_label = "Group",
    est_label = "Points per SD",
    title = "Expected Sensitivity by Diagnostic Group",
    subtitle = "Raw MADRS points per SD of true severity"
  )

Expected Sensitivity by Diagnostic Group
Raw MADRS points per SD of true severity
Group	Points per SD	95% CI
Depression	10.65	[9.45, 11.89]
Bipolar	11.30	[9.98, 12.66]
Psychosis	11.67	[10.28, 13.11]
Other	18.44	[11.89, 23.33]

contrast_draws <- slope_draws |>
  select(diagnosis, .draw, slope_points) |>
  pivot_wider(names_from = diagnosis, values_from = slope_points) |>
  mutate(diff = Other - Depression)

contrast_draws |>
  summarize(
    Contrast = "Other - Depression",
    Est = mean(diff),
    CI_Label = sprintf("[%.2f, %.2f]", quantile(diff, 0.025), quantile(diff, 0.975)),
    p_val = 2 * (1 - as.numeric(bayestestR::p_direction(diff)$pd))
  ) |>
  gt() |>
  fmt_number(columns = Est, decimals = 2) |>
  text_transform(
    locations = cells_body(columns = p_val),
    fn = \(x) {
      vals <- suppressWarnings(as.numeric(x))
      ifelse(vals < 0.001, "< .001", sprintf("%.3f", vals))
    }
  ) |>
  cols_label(Est = "Est.", CI_Label = "95% CI", p_val = "p-val") |>
  tab_style(
    style = cell_text(weight = "bold"),
    locations = cells_body(columns = p_val, rows = p_val < 0.05)
  ) |>
  tab_header(
    title = "Sensitivity Contrast in Raw Points",
    subtitle = "Difference in expected points per SD of true severity"
  ) |>
  tab_options(table.font.size = 14)

Sensitivity Contrast in Raw Points
Difference in expected points per SD of true severity
Contrast	Est.	95% CI	p-val
Other - Depression	7.78	[1.08, 12.88]	0.025

Bias and Precision by Education

edu_labels <- levels(ready_sums$education)

grid_edu <- expand_grid(education0 = 0:4, labelz = 0)

fit_edu |>
  epred_draws(newdata = grid_edu, re_formula = NA) |>
  ungroup() |>
  summarize(
    .by = education0,
    Est = mean(.epred),
    CI_Low = quantile(.epred, 0.025),
    CI_High = quantile(.epred, 0.975)
  ) |>
  arrange(education0) |>
  mutate(Education = edu_labels[education0 + 1]) |>
  gt_est_table(
    group_col = "Education", group_label = "Education Level",
    est_label = "Expected Score",
    title = "Expected MADRS Prediction by Education",
    subtitle = "At average severity (labelz = 0)"
  )

Expected MADRS Prediction by Education
At average severity (labelz = 0)
Education Level	Expected Score	95% CI
Less than High School	19.10	[17.86, 20.36]
High School/GED	18.21	[17.39, 19.03]
Part College or 2-year degree	17.34	[16.75, 17.92]
4-year College degree	16.50	[15.77, 17.23]
Part or completed Graduate degree	15.69	[14.65, 16.76]

grid_ext <- tibble(education0 = c(0, 4), labelz = 0)

mu_draws <- posterior_epred(fit_edu, newdata = grid_ext, re_formula = NA) / 60
phi_draws <- posterior_linpred(fit_edu, newdata = grid_ext, re_formula = NA,
                               dpar = "phi", transform = TRUE)
# Beta-binomial SD in raw MADRS points (n = 60 trials)
sd_draws <- sqrt(60 * mu_draws * (1 - mu_draws) * (1 + 59 / (phi_draws + 1)))

tibble(
  education0 = rep(c(0, 4), each = nrow(sd_draws)),
  sd = c(sd_draws[, 1], sd_draws[, 2])
) |>
  summarize(
    .by = education0,
    Est = mean(sd),
    CI_Low = quantile(sd, 0.025),
    CI_High = quantile(sd, 0.975)
  ) |>
  arrange(education0) |>
  mutate(Education = edu_labels[education0 + 1]) |>
  gt_est_table(
    group_col = "Education", group_label = "Education Level",
    est_label = "Prediction SD",
    title = "Model-Implied Prediction Noise by Education",
    subtitle = "SD of predictions in raw MADRS points"
  )

Model-Implied Prediction Noise by Education
SD of predictions in raw MADRS points
Education Level	Prediction SD	95% CI
Less than High School	4.01	[3.67, 4.56]
Part or completed Graduate degree	4.90	[4.18, 5.72]

Precision by Interviewer

intv_levels <- levels(ready_sums$interviewer)
grid_intv <- tibble(interviewer = intv_levels, labelz = 0)

mu_intv <- posterior_epred(fit_intv, newdata = grid_intv, re_formula = NA) / 60
phi_intv <- posterior_linpred(fit_intv, newdata = grid_intv, re_formula = NA,
                              dpar = "phi", transform = TRUE)
sd_intv <- sqrt(60 * mu_intv * (1 - mu_intv) * (1 + 59 / (phi_intv + 1)))

tibble(
  interviewer = rep(intv_levels, each = nrow(sd_intv)),
  sd = as.numeric(sd_intv)
) |>
  summarize(
    .by = interviewer,
    Est = mean(sd),
    CI_Low = quantile(sd, 0.025),
    CI_High = quantile(sd, 0.975)
  ) |>
  gt_est_table(
    group_col = "interviewer", group_label = "Interviewer",
    est_label = "Prediction SD",
    title = "Model-Implied Prediction Noise by Interviewer",
    subtitle = "SD of predictions in raw MADRS points"
  )

Model-Implied Prediction Noise by Interviewer
SD of predictions in raw MADRS points
Interviewer	Prediction SD	95% CI
Interviewer 1	4.06	[3.60, 4.63]
Interviewer 2	4.87	[4.35, 5.43]
Interviewer 3	3.94	[3.53, 4.99]
Interviewer 4	3.68	[3.43, 4.34]