5  Example: Tax Preferences

This chapter applies the structural estimator to a tax-plan conjoint experiment from Ballard-Rosa, Martin, and Scheve (2017). Respondents evaluate pairs of hypothetical income tax plans varying on six marginal tax rates (one per income bracket) plus a revenue impact indicator, and select the plan they prefer.

What makes this dataset unique among the bundled examples is that every attribute is numeric (continuous). There are no factor attributes — the six bracket rates are percentages and the revenue score is a centered integer scale. This showcases how sconjoint handles all-numeric attribute designs.

5.1 Data preparation

data(br2017, package = "sconjoint")
dim(br2017)
[1] 32000    23
head(br2017, 4)
  respondent task profile choice rate_L10 rate_10_35 rate_35_85 rate_85_175
1  177191531    1       1      1        5         25         35          25
2  177191531    1       2      0        0         35          5          15
3  177191531    2       1      1       15         25          5          25
4  177191531    2       2      0       15         15         15          15
  rate_175_375 rate_375P revenue_score resp_age resp_female resp_pid7
1           15         5             0       39           0         5
2            5        45             1       39           0         5
3           25        25             0       39           0         5
4           35         5            -1       39           0         5
  resp_educ resp_race_white resp_income resp_ineq_averse resp_work_vs_luck
1         2               0           2                0                 2
2         2               0           2                0                 2
3         2               0           2                0                 2
4         2               0           2                0                 2
  resp_taxes_harm resp_hardwork resp_high_econ_know resp_employed_ft
1               2             0                   0                1
2               2             0                   0                1
3               2             0                   0                1
4               2             0                   0                1
str(br2017)
'data.frame':   32000 obs. of  23 variables:
 $ respondent         : chr  "177191531" "177191531" "177191531" "177191531" ...
 $ task               : int  1 1 2 2 3 3 4 4 5 5 ...
 $ profile            : int  1 2 1 2 1 2 1 2 1 2 ...
 $ choice             : int  1 0 1 0 1 0 1 0 1 0 ...
 $ rate_L10           : num  5 0 15 15 25 25 15 25 5 25 ...
 $ rate_10_35         : num  25 35 25 15 15 35 35 35 25 25 ...
 $ rate_35_85         : num  35 5 5 15 15 35 35 5 5 35 ...
 $ rate_85_175        : num  25 15 25 15 35 5 15 25 25 15 ...
 $ rate_175_375       : num  15 5 25 35 35 25 25 25 15 5 ...
 $ rate_375P          : num  5 45 25 5 55 5 45 55 15 25 ...
 $ revenue_score      : num  0 1 0 -1 1 2 2 2 -1 2 ...
 $ resp_age           : num  39 39 39 39 39 39 39 39 39 39 ...
 $ resp_female        : num  0 0 0 0 0 0 0 0 0 0 ...
 $ resp_pid7          : num  5 5 5 5 5 5 5 5 5 5 ...
 $ resp_educ          : num  2 2 2 2 2 2 2 2 2 2 ...
 $ resp_race_white    : num  0 0 0 0 0 0 0 0 0 0 ...
 $ resp_income        : num  2 2 2 2 2 2 2 2 2 2 ...
 $ resp_ineq_averse   : num  0 0 0 0 0 0 0 0 0 0 ...
 $ resp_work_vs_luck  : num  2 2 2 2 2 2 2 2 2 2 ...
 $ resp_taxes_harm    : num  2 2 2 2 2 2 2 2 2 2 ...
 $ resp_hardwork      : num  0 0 0 0 0 0 0 0 0 0 ...
 $ resp_high_econ_know: num  0 0 0 0 0 0 0 0 0 0 ...
 $ resp_employed_ft   : num  1 1 1 1 1 1 1 1 1 1 ...

Each respondent contributes eight forced-choice tasks with two tax-plan profiles each (16 rows per respondent). The respondent-level moderators include demographics (age, gender, education, race, income), partisanship (resp_pid7, 1–7), and economic attitudes (inequality aversion, work-vs-luck beliefs, taxes-harm beliefs, hard-work beliefs, economic knowledge, employment status) — 12 moderators in total.

5.2 Fitting the structural model

Every attribute here is continuous, with the tax rates encoded as percentage points (0–50). The paper handles this with the varref Stage-2 prior and a low variance floor (varref_floor = 1e-3) on the raw rates, rather than the default score-based map_c5 prior. We follow the paper’s recipe here. Coefficients are in “logit utility per percentage point” units.

fit_br <- scfit(
  choice ~ rate_L10 + rate_10_35 + rate_35_85 + rate_85_175 +
           rate_175_375 + rate_375P + revenue_score |
           resp_age + resp_female + resp_pid7 + resp_educ +
           resp_race_white + resp_income + resp_ineq_averse +
           resp_work_vs_luck + resp_taxes_harm + resp_hardwork +
           resp_high_econ_know + resp_employed_ft,
  data         = br2017,
  respondent   = "respondent",
  task         = "task",
  profile      = "profile",
  K            = 5L,
  n_epochs     = 200L,
  seed         = 2024,
  stage2       = "varref",
  varref_floor = 1e-3
)
summary(fit_br)
sc_fit summary
Call: scfit(formula = choice ~ rate_L10 + rate_10_35 + rate_35_85 + 
    rate_85_175 + rate_175_375 + rate_375P + revenue_score | 
    resp_age + resp_female + resp_pid7 + resp_educ + resp_race_white + 
        resp_income + resp_ineq_averse + resp_work_vs_luck + 
        resp_taxes_harm + resp_hardwork + resp_high_econ_know + 
        resp_employed_ft, data = br2017, respondent = "respondent", 
    task = "task", profile = "profile", K = 5L, n_epochs = 200L, 
    seed = 2024, stage2 = "varref", varref_floor = 0.001)

2000 respondents | 16000 observations | K = 5 folds
hidden = 32-32-16 | epochs = 200 | seed = 2024 | device = cpu

Coefficients (DML, respondent-clustered SE):
                estimate std_error  z_value   p_value     ci_lo     ci_hi
rate_L10      -0.0426721  0.002082 -20.4981 2.240e-93 -0.046752 -0.038592
rate_10_35    -0.0320701  0.002889 -11.1018 1.230e-28 -0.037732 -0.026408
rate_35_85    -0.0188041  0.001633 -11.5157 1.099e-30 -0.022005 -0.015604
rate_85_175   -0.0008085  0.001180  -0.6855 4.930e-01 -0.003120  0.001503
rate_175_375   0.0070509  0.001047   6.7344 1.646e-11  0.004999  0.009103
rate_375P      0.0105807  0.001063   9.9504 2.513e-23  0.008497  0.012665
revenue_score  0.1362353  0.026182   5.2034 1.957e-07  0.084919  0.187551

DML/iid SE ratio (mean): 1.061

Stage 2: varref | mean(sigma_prior) = 0.001
NoteWhy stage2 = "varref" for continuous attributes?

The default map_c5 Stage-2 prior is calibrated assuming \(\mathrm{Var}(\Delta X_k) \approx 1\) (true for factor dummies under typical randomization). BR’s tax rates span 0–50 percentage points, so the score-based prior becomes loose and the MAP solver leaves per-respondent estimates with extreme tails. The paper instead uses the varref prior — \(0.5\cdot\mathrm{Var}_i(\hat\beta_{\text{ens},i,k})\), floored at varref_floor — on the raw rates. The low floor matters: the prior 0.01 default over-shrinks continuous coefficients and badly degrades recovery of individual progressivity, so continuous designs use varref_floor = 1e-3 (the paper’s setting). An alternative is normalize_deltaX = TRUE with the default map_c5 prior, which standardizes \(\Delta X\) internally; see the Advanced options chapter. For factor-dummy designs (SW, GS) neither is needed.

plot(fit_br, "loss_trace")

5.3 Population-average estimates

DML coefficients per attribute, on the logit scale, with respondent- clustered SEs. With numeric attributes the coefficient is “logit utility per unit of the attribute” (per percentage-point of the rate, or per unit of the centered revenue indicator).

plot_amce(fit_br, groups = br_groups, labels = br_labels)

5.3.1 Structural AME vs reduced-form AMCE

For continuous attributes the AME is “average effect on Pr(choice) per unit of the attribute”. Overlay the structural and LPM versions on the probability scale.

ame_struct <- sc_average(fit_br, scale = "probability")
ame_df     <- ame_struct$estimate; ame_df$source <- "Structural AME"
lpm        <- sc_baseline_lpm(fit_br)
amce_df    <- data.frame(
  dummy_name = names(stats::coef(lpm)),
  estimate   = unname(stats::coef(lpm)),
  se         = sqrt(diag(stats::vcov(lpm))),
  source     = "LPM AMCE",
  stringsAsFactors = FALSE
)
amce_df$ci_lo <- amce_df$estimate - 1.96 * amce_df$se
amce_df$ci_hi <- amce_df$estimate + 1.96 * amce_df$se
df_both <- rbind(
  ame_df[, c("dummy_name", "estimate", "se", "ci_lo", "ci_hi", "source")],
  amce_df[, c("dummy_name", "estimate", "se", "ci_lo", "ci_hi", "source")]
)
df_both$label <- ifelse(df_both$dummy_name %in% names(br_labels),
                        br_labels[df_both$dummy_name],
                        df_both$dummy_name)
df_both$label <- factor(df_both$label, levels = rev(unname(br_labels)))
df_both$group <- br_groups[df_both$dummy_name]
df_both$group <- factor(df_both$group, levels = unique(br_groups))

ggplot(df_both, aes(x = estimate, y = label, color = source)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
  geom_pointrange(aes(xmin = ci_lo, xmax = ci_hi),
                  position = position_dodge(width = 0.5),
                  size = 0.35, fatten = 2) +
  scale_color_manual(values = c("Structural AME" = "#2166AC",
                                "LPM AMCE"        = "#B2182B")) +
  facet_grid(group ~ ., scales = "free_y", space = "free_y") +
  labs(x = "Effect on Pr(choice)",
       y = NULL, color = NULL,
       title = "Structural AME vs LPM AMCE (probability scale)") +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom",
        strip.text.y    = element_text(angle = 0, hjust = 0, face = "bold"))

5.4 Individual-level preferences

Ridgeline densities of \(\hat\beta_k(\mathbf Z_i)\) per attribute across respondents. With normalize_deltaX = TRUE the MAP refinement is properly scaled — the hybrid view shows clean per-respondent distributions for every bracket.

plot(fit_br, "beta_ridgelines", groups = br_groups, labels = br_labels)

5.5 Structural quantities

5.5.1 Attribute importance

sc_importance(fit_br)
sc_quantity: importance
  estimate: data.frame with 7 rows
     attribute   share        se   ci_lo   ci_hi
      rate_L10 0.28630 0.0037662 0.27891 0.29368
    rate_10_35 0.22847 0.0038988 0.22083 0.23611
    rate_35_85 0.11375 0.0029272 0.10802 0.11949
   rate_85_175 0.04309 0.0014680 0.04021 0.04597
  rate_175_375 0.05320 0.0016437 0.04998 0.05642
     rate_375P 0.22234 0.0048923 0.21275 0.23193
 revenue_score 0.05285 0.0008282 0.05122 0.05447
plot_importance(fit_br, labels = c(
  rate_L10 = "Under $10k", rate_10_35 = "$10-35k",
  rate_35_85 = "$35-85k", rate_85_175 = "$85-175k",
  rate_175_375 = "$175-375k", rate_375P = "Over $375k",
  revenue_score = "Revenue"
))

5.5.2 Direction and intensity

sc_direction_intensity(fit_br)
sc_quantity_bivariate: direction_intensity
-- direction --
sc_quantity: direction
  estimate: data.frame with 7 rows
    dummy_name      d     se_d  ci_lo_d  ci_hi_d
      rate_L10 -0.995 0.002234 -0.99938 -0.99062
    rate_10_35 -0.948 0.007119 -0.96195 -0.93405
    rate_35_85 -0.753 0.014717 -0.78185 -0.72415
   rate_85_175 -0.007 0.022366 -0.05084  0.03684
  rate_175_375  0.433 0.020161  0.39349  0.47251
     rate_375P  0.390 0.020595  0.34963  0.43037
 revenue_score  1.000 0.000000  1.00000  1.00000
-- intensity --
sc_quantity: intensity
  estimate: data.frame with 7 rows
    dummy_name       u      se_u ci_lo_u ci_hi_u
      rate_L10 0.04154 0.0003216 0.04091 0.04218
    rate_10_35 0.03087 0.0003421 0.03020 0.03154
    rate_35_85 0.01960 0.0002885 0.01903 0.02016
   rate_85_175 0.01091 0.0001834 0.01055 0.01127
  rate_175_375 0.01203 0.0001928 0.01165 0.01240
     rate_375P 0.01943 0.0003118 0.01882 0.02004
 revenue_score 0.11425 0.0005320 0.11321 0.11529

5.5.3 Fraction preferring

plot_fraction(fit_br, groups = br_groups, labels = br_labels)

5.5.4 Heterogeneity test

sc_heterogeneity_test(fit_br, adjust = "bh")
sc_quantity: heterogeneity_test
  estimate: data.frame with 7 rows
    dummy_name  var_beta    se_var t_stat    p_value p_adjusted sig
      rate_L10 0.0002090 6.331e-06  33.01 3.441e-239 4.015e-239 ***
    rate_10_35 0.0002517 7.262e-06  34.66 1.616e-263 3.770e-263 ***
    rate_35_85 0.0002274 6.501e-06  34.98 1.972e-268 6.901e-268 ***
   rate_85_175 0.0001861 5.773e-06  32.24 2.570e-228 2.570e-228 ***
  rate_175_375 0.0001641 4.830e-06  33.97 3.152e-253 5.516e-253 ***
     rate_375P 0.0004372 1.184e-05  36.94 5.014e-299 3.510e-298 ***
 revenue_score 0.0005657 1.695e-05  33.38 1.516e-244 2.123e-244 ***
plot_hetero(fit_br, groups = br_groups, labels = br_labels)

5.5.5 Subgroup analysis by party

pid_col <- fit_br$Z[, "resp_pid7"]
plot_subgroup(
  fit_br,
  subgroup = list(
    Democrat    = pid_col <= 3,
    Independent = pid_col >= 4 & pid_col <= 4,
    Republican  = pid_col >= 5
  ),
  groups = br_groups, labels = br_labels,
  title = "Subgroup AMCE by party"
)

5.5.6 Tax schedule preferences by party

This plot is unique to the tax-plan conjoint. It extracts the six bracket coefficients \(\hat\beta_k(\mathbf Z_i)\) and shows the respondent-level preferred-rate schedule in three panels (one per party). The dark ribbon is the interquartile range across respondents within each party; the light ribbon is the 10th–90th percentile range.

## One row per respondent on the hybrid (MAP-refined) scale
beta_hat <- predict(fit_br, type = "beta")
first_row <- !duplicated(fit_br$respondent_id)
beta_resp <- beta_hat[first_row, ]

## Party assignment from moderator
pid <- fit_br$Z[first_row, "resp_pid7"]
party <- ifelse(pid <= 3, "Democrat",
         ifelse(pid >= 5, "Republican", "Independent"))

## Build schedule data: 6 brackets as x-axis, beta as y-axis
brackets <- c("Under $10k", "$10-35k", "$35-85k",
              "$85-175k", "$175-375k", "Over $375k")
rate_cols <- c("rate_L10", "rate_10_35", "rate_35_85",
               "rate_85_175", "rate_175_375", "rate_375P")

quant_data <- do.call(rbind, lapply(unique(party), function(p) {
  do.call(rbind, lapply(seq_along(brackets), function(j) {
    vals <- beta_resp[party == p, rate_cols[j]]
    data.frame(
      bracket = brackets[j],
      party   = p,
      median  = median(vals),
      q25     = quantile(vals, 0.25),
      q75     = quantile(vals, 0.75),
      q10     = quantile(vals, 0.10),
      q90     = quantile(vals, 0.90)
    )
  }))
}))
quant_data$bracket <- factor(quant_data$bracket, levels = brackets)
quant_data$party   <- factor(quant_data$party,
                              levels = c("Democrat", "Independent", "Republican"))

party_palette <- c(Democrat = "#2166AC", Independent = "gray50",
                   Republican = "#B2182B")

ggplot(quant_data, aes(x = bracket, color = party, fill = party, group = party)) +
  geom_ribbon(aes(ymin = q10, ymax = q90), alpha = 0.15, color = NA) +
  geom_ribbon(aes(ymin = q25, ymax = q75), alpha = 0.30, color = NA) +
  geom_line(aes(y = median), linewidth = 1.0) +
  geom_point(aes(y = median), size = 1.8) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  facet_wrap(~ party, ncol = 3) +
  scale_color_manual(values = party_palette, guide = "none") +
  scale_fill_manual(values  = party_palette, guide = "none") +
  labs(x = "Income bracket",
       y = expression(hat(beta)[k](Z) ~ "(utility per % point)"),
       title = "Tax schedule preferences, faceted by party") +
  theme_minimal(base_size = 11) +
  theme(axis.text.x = element_text(angle = 30, hjust = 1),
        strip.text  = element_text(face = "bold"))

Nearly all respondents exhibit progressive tax preferences (a positive slope from low to high brackets), and progressivity remains the majority pattern even among Republicans. The partisan gap is real but modest — Democrats’ mean slope is somewhat steeper than Republicans’. (These shares are computed on the bundled data; read the exact values off the rendered output rather than quoting fixed numbers.)

5.6 Validating against the homogeneous-logit AMCE

v_br <- sc_validate_amce(fit_br)
v_br
sc_validate_amce -- pooled and (optionally) subgroup comparison
Stage 2: varref
N obs: 16000, N respondents: 2000, P attributes: 7
Pooled correlation (DML theta vs homogeneous logit coef): 1

Pooled comparison (first 10 rows):
     attribute dml_theta  dml_se homog_logit_coef homog_logit_se      diff
      rate_L10 -0.042672 0.00208         -0.04093       0.000434 -0.001745
    rate_10_35 -0.032070 0.00289         -0.03108       0.000626 -0.000990
    rate_35_85 -0.018804 0.00163         -0.01803       0.000353 -0.000779
   rate_85_175 -0.000809 0.00118         -0.00117       0.000263  0.000365
  rate_175_375  0.007051 0.00105          0.00669       0.000233  0.000357
     rate_375P  0.010581 0.00106          0.00975       0.000229  0.000835
 revenue_score  0.136235 0.02618          0.12774       0.005783  0.008491
 abs_diff
 0.001745
 0.000990
 0.000779
 0.000365
 0.000357
 0.000835
 0.008491

Under correct logit specification, the DML \(\hat\theta\) on each tax bracket matches the pooled logit’s coefficient on the same bracket exactly. The pooled correlation is the cleanest summary that the structural model loses nothing the AMCE delivers.