Getting started with dapper: randomized response

This vignette adapts Example 1 of dapper: Data Augmentation for Private Posterior Estimation in R, in which dapper is used to analyze a privatized record-level release produced by randomized response. The example uses a subset of the UC Berkeley admissions data aggregated across departments (Bickel et al. 1975) to demonstrate privacy-aware inference for the odds ratio.

Prepare the data and privacy mechanism

We start with the four confidential cell counts in datasets::UCBAdmissions, summed over departments. For computational convenience, we take a subsample of 400 applicants. The first column encodes sex (male = 1, female = 0), and the second encodes admission status (admitted = 1, rejected = 0).

library(dapper)
library(ggplot2)

cells <- data.frame(
  sex = c(1, 1, 0, 0),
  status = c(1, 0, 1, 0)
)
counts <- c(1198, 1493, 557, 1278)
cnf_df <- cells[rep(seq_len(nrow(cells)), times = counts), ]

set.seed(1)
n <- 400
ix <- sample(seq_len(nrow(cnf_df)), n, replace = FALSE)
cnf_df <- cnf_df[ix, ]

Randomized response (Warner 1965) is applied independently to both binary attributes. Each released value matches the confidential value with probability 3/4 and is flipped with probability 1/4. This can be done using, for example, two flips of a fair coin. Applying this mechanism to both attributes gives a privacy budget of at most \(\epsilon = 2\log(3)\) for replacing one individual’s pair of attributes.

We store the release in sdp, a length-800 vector formed by stacking the sex column above the admission-status column.

ri <- as.logical(rbinom(2 * n, 1, 1/2))
ra <- rbinom(sum(ri), 1, 1/2)

sdp <- c(as.matrix(cnf_df))
sdp[ri] <- ra
prv_df <- data.frame(sex = sdp[seq_len(n)], status = sdp[n + seq_len(n)])

The following tables show the confidential and privatized counts. In an actual private-data analysis, the analyst would have access only to the privatized release and the known privacy mechanism.

admission_table <- function(x) {
  table(
    Sex = factor(x$sex, levels = c(0, 1), labels = c("Female", "Male")),
    Status = factor(x$status, levels = c(1, 0),
                    labels = c("Admitted", "Rejected"))
  )
}

knitr::kable(admission_table(cnf_df), caption = "Confidential counts")
Confidential counts
Admitted Rejected
Female 46 118
Male 109 127 
knitr::kable(admission_table(prv_df), caption = "Privatized counts")
Privatized counts
Admitted Rejected
Female 74 102
Male 104 120

Specify the four model components

Let \(\theta = (\pi_{11}, \pi_{10}, \pi_{01}, \pi_{00})\) contain the joint probabilities of the four sex/status combinations, in the order used by cells. We place a Dirichlet(1, 1, 1, 1) prior on these probabilities.

new_privacy() requires four functions, with their argument names and order as shown below.

Generate latent data: latent_f(theta)

Given the current probabilities, generate a proposed confidential dataset. Sampling from the four possible rows gives an \(n \times 2\) matrix.

latent_f <- function(theta) {
  tl <- list(c(1, 1), c(1, 0), c(0, 1), c(0, 0))
  rs <- sample(tl, n, replace = TRUE, prob = theta)
  do.call(rbind, rs)
}

Update the parameters: posterior_f(dmat, theta)

Conditional on the confidential data, the posterior is Dirichlet with parameters equal to the four cell counts plus one. Normalized gamma draws produce a draw from this posterior. The theta argument is required by the interface but is unused for this conjugate update.

posterior_f <- function(dmat, theta) {
  sex <- dmat[, 1]
  status <- dmat[, 2]
  x <- c(
    sum(sex & status),
    sum(sex & !status),
    sum(!sex & status),
    sum(!sex & !status)
  )
  t1 <- rgamma(4, shape = x + 1, rate = 1)
  t1 / sum(t1)
}

Summarize one record: statistic_f(xi, sdp, i)

For row \(i\), count the attributes that match their released values:

\[ t_i(x_i, s_{dp}) = \mathbb{1}(x_{i1} = s_{dp,i}) + \mathbb{1}(x_{i2} = s_{dp,n+i}). \]

The sampler sums these contributions across rows. This allows it to update the total efficiently when proposing a change to one latent record.

statistic_f <- function(xi, sdp, i) {
  n <- length(sdp) %/% 2
  sum(xi == sdp[c(i, n + i)])
}

Evaluate the privacy likelihood: mechanism_f(sdp, sx)

If sx is the total number of matching attributes, then length(sdp) - sx is the number that differ. The log likelihood of the randomized-response release is therefore

\[ \texttt{sx}\log(3/4) + (2n - \texttt{sx})\log(1/4). \]

mechanism_f <- function(sdp, sx) {
  sx * log(3/4) + (length(sdp) - sx) * log(1/4)
}

Run the sampler and inspect the output

Combine the components with new_privacy(), then pass the model and release to dapper_sample(). As in the paper, we run four chains with 6,000 iterations each and discard the first 1,000 iterations per chain as warmup. This leaves 20,000 draws in total. By default, the chains run sequentially.

dmod <- new_privacy(
  posterior_f = posterior_f,
  latent_f = latent_f,
  mechanism_f = mechanism_f,
  statistic_f = statistic_f,
  npar = 4,
  varnames = c("pi_11", "pi_10", "pi_01", "pi_00")
)

dp_out <- dapper_sample(
  dmod,
  sdp = sdp,
  seed = 123,
  niter = 6000,
  warmup = 1000,
  chains = 4,
  init_par = rep(0.25, 4)
)

For parallel execution, see the example in ?dapper_sample. Progress can be monitored by wrapping the sampling call in progressr::with_progress().

summary() reports posterior summaries and diagnostics using the posterior package. Inspect the R-hat values, effective sample sizes, and trace plots together before drawing conclusions. The iteration count here is chosen for a short demonstration. More iterations may be needed for precise estimates of tail probabilities or quantiles.

summary(dp_out)
#> # A tibble: 4 × 10
#>   variable  mean median     sd    mad     q5   q95  rhat ess_bulk ess_tail
#>   <chr>    <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 pi_11    0.281  0.281 0.0610 0.0625 0.182  0.382  1.02     362.     818.
#> 2 pi_10    0.336  0.335 0.0638 0.0640 0.235  0.444  1.01     431.    1191.
#> 3 pi_01    0.111  0.108 0.0548 0.0563 0.0250 0.206  1.02     282.     504.
#> 4 pi_00    0.272  0.272 0.0601 0.0616 0.172  0.372  1.02     389.     875.

plot() produces trace plots through bayesplot::mcmc_trace(). Privacy noise can make the chains mix slowly, so checking these plots is particularly useful.

plot(dp_out) +
  scale_color_manual(values = c("#00468B", "#ED0000", "#42B540", "#925E9F")) +
  theme_bw(base_size = 11)
Four panels show the sampled cell probabilities over iterations, with a different color for each chain.

Trace plots for the four cell probabilities. Colors distinguish the four chains.

Estimate the odds ratio

The odds ratio compares the odds of admission for males with those for females:

\[ \mathrm{OR} = \frac{\pi_{11}/\pi_{10}}{\pi_{01}/\pi_{00}}. \]

Values above one indicate higher aggregate odds of admission for males. Because dp_out$chain is a posterior draws_matrix, we can create this derived quantity and summarize it directly.

odds_ratio_draws <- posterior::mutate_variables(
  dp_out$chain,
  odds_ratio = (pi_11 * pi_00) / (pi_10 * pi_01)
)
odds_ratio_draws <- posterior::subset_draws(
  odds_ratio_draws, variable = "odds_ratio"
)
posterior::summarise_draws(odds_ratio_draws)
#> # A tibble: 1 × 10
#>   variable    mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
#>   <chr>      <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 odds_ratio  6.83   2.12  82.0  1.89 0.472  16.7  1.03     245.     568.

The shaded region below marks the central 80% posterior interval, while the outer curve covers the central 95%. The dashed line marks an odds ratio of one.

bayesplot::mcmc_areas(
  odds_ratio_draws,
  prob = 0.8,
  prob_outer = 0.95,
  point_est = "none"
) +
  geom_vline(xintercept = 1, linetype = "dashed", linewidth = 0.5) +
  coord_cartesian(xlim = c(0, 10)) +
  labs(x = "Odds ratio", y = NULL) +
  theme_bw(base_size = 11) +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())
Posterior density of the odds ratio, with a shaded central 80 percent interval and a dashed vertical reference line at one.

Privacy-aware posterior distribution of the odds ratio, displayed from 0 to 10.

Compare with analyses that omit the privacy mechanism

To examine the impact of the privacy mechanism, we can compare the posterior draws from dapper to an equivalent Bayesian analysis on the confidential data. Moreover, we can see the consequence of ignoring the privacy mechanism by repeating the analysis on the privatized records, but treating them as if they were the confidential data.

set.seed(1)
confidential_data <- as.matrix(cnf_df)
cps <- t(replicate(20000, posterior_f(confidential_data, NULL)))
odds_ratio_conf <- (cps[, 1] * cps[, 4]) / (cps[, 2] * cps[, 3])

set.seed(1)
noisy_data <- matrix(sdp, ncol = 2)
cps <- t(replicate(20000, posterior_f(noisy_data, NULL)))
odds_ratio_noisy <- (cps[, 1] * cps[, 4]) / (cps[, 2] * cps[, 3])

odds_ratio_private <- posterior::extract_variable(
  odds_ratio_draws, "odds_ratio"
)

The privacy-aware odds-ratio posterior has a long right tail, so we evaluate each density on a common plotting grid from 0 to 8, using all draws to estimate the density.

posterior_density <- function(draws, data, analysis) {
  estimate <- density(draws, from = 0, to = 8)
  data.frame(odds_ratio = estimate$x, density = estimate$y,
             data = data, analysis = analysis)
}

comparison <- rbind(
  posterior_density(odds_ratio_conf, "Confidential", "Privacy-aware"),
  posterior_density(odds_ratio_private, "Privatized", "Privacy-aware"),
  posterior_density(odds_ratio_conf, "Confidential", "Naïve"),
  posterior_density(odds_ratio_noisy, "Privatized", "Naïve")
)
comparison$analysis <- factor(comparison$analysis,
                              levels = c("Privacy-aware", "Naïve"))

Each panel uses the confidential-data posterior as a reference. The left compares it with the privacy-aware posterior from dapper; the right compares it with the naïve posterior that results from treating the privatized data as if they were confidential.

data_colors <- c(Confidential = "#00468B", Privatized = "#ED0000")

ggplot(comparison, aes(x = odds_ratio, y = density, fill = data, color = data)) +
  geom_area(position = "identity", alpha = 0.2, color = NA) +
  geom_line(linewidth = 0.6) +
  facet_wrap(~ analysis) +
  coord_cartesian(xlim = c(0, 8)) +
  scale_fill_manual(name = "Data", values = data_colors) +
  scale_color_manual(name = "Data", values = data_colors) +
  labs(x = "Odds ratio", y = "Density") +
  theme_bw(base_size = 11) +
  theme(legend.position = "bottom")
Two density panels compare confidential and privatized analyses. The privacy-aware posterior is wider, while the naïve posterior is concentrated closer to one.

Odds-ratio posteriors for confidential data (blue) and privatized data (red). The left panel accounts for randomized response; the right panel ignores it.

The plots demonstrate that the naïve posterior exhibits some bias towards one and is overconfident. On the other hand, the privacy-aware analysis accounts for the additional uncertainty introduced by randomized response and produces a wider posterior. This comparison illustrates why the privacy mechanism should be accounted for in a statistical analysis.

References

Bickel, Peter J., Eugene A. Hammel, and John W. O’Connell. 1975. “Sex Bias in Graduate Admissions: Data from Berkeley: Measuring Bias Is Harder Than Is Usually Assumed, and the Evidence Is Sometimes Contrary to Expectation.” Science 187 (4175): 398–404. https://doi.org/10.1126/science.187.4175.398.
Warner, Stanley L. 1965. “Randomized Response: A Survey Technique for Eliminating Evasive Answer Bias.” Journal of the American Statistical Association 60 (309): 63–69. https://doi.org/10.1080/01621459.1965.10480775.