DP Regression and Privacy-Aware Confidence Intervals

library(DPrivStats)
set.seed(7)

Closed-form DP linear regression

dp_lm() fits OLS on clipped data and releases coefficients with the analytic Gaussian mechanism calibrated to the L2 sensitivity \(\Delta_2 = C \cdot D / \lambda_{\min}(X^\top X)\) (see ?dp_lm_sensitivity).

data(example_microdata)
fit <- dp_lm(income ~ education + age + hours, example_microdata,
             epsilon = 2.0, delta = 1e-6,
             bounds = list(y = c(0, 500000)))
fit
#> 
#> Differentially Private Linear Regression
#> epsilon = 2.000, delta = 1.0e-06
#> 
#> (Intercept)   education         age       hours 
#>     6589456    -3460229    -2008413    -1192446

Privacy-aware confidence intervals

Standard intervals fail under DP because the point estimate is biased by privacy noise and the standard error is itself privatized. dp_confint() adds the known privacy-noise variance to the sampling variance:

ci <- dp_confint(fit)
ci
#>             estimate     lower    upper
#> (Intercept)  6589456 -49344972 62523883
#> education   -3460229  -9636144  2715686
#> age         -2008413  -7715307  3698480
#> hours       -1192446  -6941864  4556972

Bootstrap alternatives:

dp_confint(fit, method = "parametric_bootstrap", B = 200)
#>             estimate     lower    upper
#> (Intercept)  6589456 -43119727 64925930
#> education   -3460229  -9996642  1888128
#> age         -2008413  -7612340  3207671
#> hours       -1192446  -6797279  4294436

DP-SGD for GLMs

For logistic regression, dp_glm() implements DP-SGD with per-sample gradient clipping:

d <- transform(example_microdata,
               high_income = as.numeric(income > 60000))
glm_fit <- dp_glm(high_income ~ education + age, d, binomial(),
                  epsilon = 2.0, delta = 1e-6,
                  max_grad_norm = 1, n_iter = 300, lr = 0.05,
                  bounds = list(y = c(0, 1)))
glm_fit
#> 
#> Differentially Private GLM (DP-SGD)
#> family = binomial, iterations = 300
#> epsilon = 2.000, delta = 1.0e-06
#> 
#> (Intercept)   education         age 
#>      0.0197     -0.0563     -0.1577

Coverage validation

validate_coverage() checks that the analytical intervals attain nominal coverage across Monte Carlo replicates:

gen <- function(n) {
  data.frame(x = rnorm(n, 2, 1))
}
set.seed(1)
vc <- validate_coverage(y ~ x, c(`(Intercept)` = 1, x = 2),
                        sigma = 1, data_gen = gen, n = 400, n_sims = 50,
                        epsilon = 3, delta = 1e-6,
                        y_bounds = c(-20, 20))
vc$coverage_rate
#> [1] 0.98