global bounds; incorporate into local CIs

Author: CoryMcCartan

DESCRIPTION

@@ -8,11 +8,11 @@ Authors@R: c(
            comment = c(ORCID = "0000-0002-5687-2647"))
     )
 Description: Efficient and user-friendly routines for modern ecological
-    inference. Implements the methods described in McCartan & Kuriwaki (2025)
-    <arxiv:2507.XXXXX>, which generalize ecological regression as introduced by
+    inference. Implements the methods described in McCartan & Kuriwaki (2025+)
+    <arxiv:2509.20194>, which generalize ecological regression as introduced by
     Goodman (1953) <doi:10.2307/2088121>. Includes routines for preprocessing,
     synthetic data generation, double/debiased machine learning (DML)
-    estimation, and sensitivity analysis.
+    estimation, partial identification bounds, and sensitivity analysis.
 Depends:
     R (>= 3.5.0)
 Imports:
@@ -41,7 +41,7 @@ RoxygenNote: 7.3.3
 LazyData: true
 SystemRequirements: C++17
 Config/testthat/edition: 3
-URL: https://github.com/CoryMcCartan/seine, https://corymccartan.com/seine/
+URL: https://corymccartan.com/seine/, https://github.com/CoryMcCartan/seine
 BugReports: https://github.com/CoryMcCartan/seine/issues
 VignetteBuilder: knitr
 Config/build/compilation-database: true

NEWS.md

@@ -1,5 +1,8 @@
 # seine 0.1.0
 
 * Double/debiased ML for ecological inference
+* Sensitivity analysis for EI
+* Partial identification bounds
+* Local (unit-level) estimates and confidence intervals
 * Preprocessing helper functions
 * Synthetic EI data

R/ei_bounds.R

@@ -18,28 +18,34 @@
 #'   containing the total number of observations in each aggregate unit. For
 #'   example, the column containing the total number of voters. Required for
 #'   computing weights unless `x` is an [ei_spec()] object.
+#' @param global If `TRUE`, aggregate the bounds across units to produce bounds
+#'   on the global estimands.
 #' @param ... Not currently used, but required for extensibility.
 #'
 #' @returns A data frame with bounds. The `.row` column in the output
-#'   corresponds to the observation index in the input. The `wt` column contains
-#'   the product of the predictor variable and total for each observation.
-#'   Taking a weighted average of the bounds against this column will produce
-#'   global bounds. The `min` and `max` columns contain the minimum and maximum
-#'   values for each local estimand. It has class `ei_bounds`.
+#'   corresponds to the observation index in the input. The `min` and `max`
+#'   columns contain the minimum and maximum values for each local estimand.
+#'   The `wt` column contains the product of the predictor variable and total
+#'   for each observation. Taking a weighted average of the bounds against this
+#'   column will produce global bounds. It has class `ei_bounds`.
 #'
 #' @examples
 #' data(elec_1968)
 #'
 #' spec = ei_spec(elec_1968, vap_white:vap_other, pres_dem_hum:pres_abs,
 #'                total = pres_total, covariates = c(state, pop_urban, farm))
 #'
-#' bounds = ei_bounds(spec, bounds = c(0, 1))
-#' print(bounds)
+#' ei_bounds(spec, bounds = c(0, 1))
+#' ei_bounds(spec, bounds = c(0, 1), global = TRUE)
 #'
-#' # aggregate min/max
+#' # Infer bounds
+#' ei_bounds(pres_ind_wal ~ vap_white, data = elec_1968, total = pres_total, bounds = NULL)
+#'
+#' # manually aggregate min/max
 #' # easier with dplyr:
 #' # summarize(across(min:max, ~ weighted.mean(.x, wt)), .by=c(predictor, outcome))
-#' do.call(rbind, lapply(split(bounds, ~ predictor + outcome), function(b) {
+#' grp_units = split(ei_bounds(spec, bounds = c(0, 1)), ~ predictor + outcome)
+#' do.call(rbind, lapply(grp_units, function(b) {
 #'     tibble::tibble(
 #'         predictor = b$predictor[1],
 #'         outcome = b$outcome[1],
@@ -48,16 +54,14 @@
 #'     )
 #' }))
 #'
-#' # Infer bounds
-#' ei_bounds(pres_ind_wal ~ vap_white, data = elec_1968, total = pres_total, bounds = NULL)
 #' @export
-ei_bounds <- function(x, ..., total, contrast = NULL, bounds = c(0, 1)) {
+ei_bounds <- function(x, ..., total, contrast = NULL, bounds = c(0, 1), global = FALSE) {
     UseMethod("ei_bounds")
 }
 
 #' @export
 #' @rdname ei_bounds
-ei_bounds.ei_spec <- function(x, total, contrast = NULL, bounds = c(0, 1), ...) {
+ei_bounds.ei_spec <- function(x, total, contrast = NULL, bounds = c(0, 1), global = FALSE, ...) {
     spec = x
     validate_ei_spec(spec)
 
@@ -75,12 +79,12 @@ ei_bounds.ei_spec <- function(x, total, contrast = NULL, bounds = c(0, 1), ...)
         total = as.numeric(eval_tidy(enquo(total), spec))
     }
 
-    ei_bounds_impl(x_mat, y_mat, total, contrast, bounds)
+    ei_bounds_bridge(x_mat, y_mat, total, contrast, bounds, global)
 }
 
 #' @export
 #' @rdname ei_bounds
-ei_bounds.formula <- function(formula, data, total, contrast = NULL, bounds = c(0, 1), ...) {
+ei_bounds.formula <- function(formula, data, total, contrast = NULL, bounds = c(0, 1), global = FALSE, ...) {
     forms = ei_forms(formula)
     form_preds = terms(rlang::new_formula(lhs = NULL, rhs = forms$predictors))
     form_out = terms(rlang::new_formula(forms$outcome, rhs = NULL))
@@ -103,13 +107,13 @@ ei_bounds.formula <- function(formula, data, total, contrast = NULL, bounds = c(
 
     total = as.numeric(eval_tidy(enquo(total), data))
 
-    ei_bounds_impl(x, y, total, contrast, processed$blueprint$bounds)
+    ei_bounds_bridge(x, y, total, contrast, processed$blueprint$bounds, global)
 }
 
 
 #' @export
 #' @rdname ei_bounds
-ei_bounds.data.frame <- function(x, y, total, contrast = NULL, bounds = c(0, 1), ...) {
+ei_bounds.data.frame <- function(x, y, total, contrast = NULL, bounds = c(0, 1), global = FALSE, ...) {
     x_mat = as.matrix(x)
     check_preds(x_mat, call = rlang::new_call(rlang::sym("ei_bounds")))
     y_mat = as.matrix(y)
@@ -122,13 +126,13 @@ ei_bounds.data.frame <- function(x, y, total, contrast = NULL, bounds = c(0, 1),
 
     bounds = check_bounds(bounds, y_mat)
 
-    ei_bounds_impl(x_mat, y_mat, total, contrast, bounds)
+    ei_bounds_bridge(x_mat, y_mat, total, contrast, bounds, global)
 }
 
 #' @export
 #' @rdname ei_bounds
-ei_bounds.matrix <- function(x, y, total, contrast = NULL, bounds = c(0, 1), ...) {
-    ei_bounds.data.frame(x, y, total, contrast, bounds, ...)
+ei_bounds.matrix <- function(x, y, total, contrast = NULL, bounds = c(0, 1), global = FALSE, ...) {
+    ei_bounds.data.frame(x, y, total, contrast, bounds, global, ...)
 }
 
 #' @export
@@ -142,37 +146,55 @@ ei_bounds.default <- function(x, ...) {
 
 # Implementation --------------------------------------------------------------
 
-ei_bounds_impl <- function(x, y, total, contrast, bounds) {
+ei_bounds_bridge <- function(x, y, total, contrast, bounds, global = FALSE) {
     n = nrow(x)
     n_x = ncol(x)
     n_y = ncol(y)
 
-    if (!is.null(contrast)) {
-        cli_abort("{.arg contrast} is not yet implemented for {.fn ei_bounds}.")
-    }
-
     if (identical(bounds, c(-Inf, Inf))) {
         cli_abort("At least one bound must be provided for {.fn ei_bounds}.", call=parent.frame())
     }
     if (any(is.na(x))) cli_abort("Missing values found in predictors.", call=parent.frame())
     if (any(is.na(y))) cli_abort("Missing values found in outcome.", call=parent.frame())
 
-    result = R_bounds_lp(x, y, bounds)
+    result = ei_bounds_impl(x, y, total, contrast, bounds)
 
     x_nm = colnames(x)
     y_nm = colnames(y)
 
-    tibble::new_tibble(
-        list(
-            .row = rep(seq_len(n), n_x * n_y),
-            predictor = rep(rep(x_nm, each = n), n_y),
-            outcome = rep(y_nm, each = n * n_x),
-            wt = if (is.null(contrast)) rep(c(x * total), n_y) else NULL,
-            min = c(result$min),
-            max = c(result$max)
-        ),
-        class = "ei_bounds"
-    )
+    if (isTRUE(global)) {
+        wt = total * (matrix(1, 1, n_y) %x% x)
+        wt = scale_cols(wt, 1 / colSums(wt))
+        tibble::new_tibble(
+            list(
+                predictor = rep(x_nm, n_y),
+                outcome = rep(y_nm, each = n_x),
+                min = colSums(result$min * wt),
+                max = colSums(result$max * wt)
+            ),
+            class = "ei_bounds"
+        )
+    } else {
+        tibble::new_tibble(
+            list(
+                .row = rep(seq_len(n), n_x * n_y),
+                predictor = rep(rep(x_nm, each = n), n_y),
+                outcome = rep(y_nm, each = n * n_x),
+                wt = if (is.null(contrast)) rep(c(x * total), n_y) else NULL,
+                min = c(result$min),
+                max = c(result$max)
+            ),
+            class = "ei_bounds"
+        )
+    }
+}
+
+ei_bounds_impl <- function(x, y, total, contrast, bounds) {
+    if (!is.null(contrast)) {
+        cli_abort("{.arg contrast} is not yet implemented for {.fn ei_bounds}.")
+    }
+
+    R_bounds_lp(x, y, as.double(bounds))
 }
 
 #' @describeIn ei_bounds Format bounds as an array with dimensions

R/ei_local.R

@@ -60,7 +60,7 @@
 #'   corresponds to the observation index in the input. The `wt` column contains
 #'   the product of the predictor variable and total for each observation.
 #'   Taking a weighted average of the estimate against this column will produce
-#'   a global estimate. It has class `ei_est_local`, supporting several methods.
+#'   a global estimate. It has class `ei_est_local`.
 #'
 #' @inherit ei_est references
 #'
@@ -118,6 +118,7 @@ ei_est_local = function(
     if (is.null(sum_one) && all(bounds == c(0, 1))) {
         sum_one = isTRUE(all.equal(rowSums(y), rep(1, nrow(y))))
     }
+    has_bounds = !identical(bounds, c(-Inf, Inf))
 
     if (missing(total)) {
         if (inherits(data, "ei_spec")) {
@@ -145,11 +146,7 @@ ei_est_local = function(
         eta_proj = eta_proj %*% contr$m
     }
 
-    sds = if (!isFALSE(conf_level)) {
-        local_sds(rl$x, b_cov, rl$vcov_u, regr$sigma2, contr$m, !is.null(contrast))
-    } else {
-        NULL
-    }
+    sds = local_sds(rl$x, b_cov, rl$vcov_u, regr$sigma2, contr$m, !is.null(contrast))
 
     ests = list(
         .row = rep(seq_len(n), length(x_nm)),
@@ -169,20 +166,25 @@ ei_est_local = function(
         class = "ei_est_local"
     )
 
+    if (has_bounds && is.null(contrast)) {
+        bb = ei_bounds_bridge(rl$x, y, total, contrast, bounds)
+        fac = if (isTRUE(unimodal)) sqrt(1/12) else 0.5
+        ests$std.error = pmin(ests$std.error, fac * (bb$max - bb$min))
+    }
+
     if (!isFALSE(conf_level)) {
         fac = if (isTRUE(unimodal)) 4 / 9 else 1
         chebyshev = sqrt(fac / (1 - conf_level))
         ests$conf.low = ests$estimate - chebyshev * ests$std.error
         ests$conf.high = ests$estimate + chebyshev * ests$std.error
 
-        if (is.null(contrast)) {
-            ests$conf.low[ests$conf.low < bounds[1]] = bounds[1]
-            ests$conf.high[ests$conf.high < bounds[1]] = bounds[1]
-            ests$conf.low[ests$conf.low > bounds[2]] = bounds[2]
-            ests$conf.high[ests$conf.high > bounds[2]] = bounds[2]
+        if (has_bounds && is.null(contrast)) {
+            ests$conf.low = pmax(ests$conf.low, bb$min)
+            ests$conf.low = pmin(ests$conf.low, bb$max)
+            ests$conf.high = pmax(ests$conf.high, bb$min)
+            ests$conf.high = pmin(ests$conf.high, bb$max)
         }
     }
-    ests$std.error = NULL
 
     ests
 }

explore/local.R

@@ -15,6 +15,13 @@ spec = ei_spec(elec_1968, c(vap_white, vap_black, vap_other), c(pres_dem_hum, pr
 m = ei_ridge(spec, bounds=F, sum_one=F)
 rr = ei_riesz(spec, penalty = m$penalty)
 
+left_join(ei_est(m, rr, spec), ei_bounds(spec, bounds=0:1, global=TRUE)) |>
+    mutate(
+        within = min <= estimate & estimate <= max,
+        se_max = (max - min) / sqrt(12),
+        se_within = std.error <= se_max
+    )
+
 (colMeans(est_check_regr(m, spec, nrow(spec), NULL, 4, T)$vcov) |>
     matrix(3, 3) * m$sigma2[1]) |>
     diag() |>
@@ -40,8 +47,8 @@ ei_est_local(m, spec, b_cov = b_cov2, bounds=c(0, 1), sum_one = TRUE) |>
 
 
 ei_est_local(m, spec, b_cov = b_cov2, bounds = c(0, 1), sum_one = TRUE, conf_level = 0.95, regr_var = T) |>
-    # dplyr::filter(.row == 592) |>
-    dplyr::filter(.row == 1) |>
+    dplyr::filter(.row == 592) |>
+    # dplyr::filter(.row == 1) |>
     ggplot(aes(estimate, paste0(outcome, " | ", predictor))) +
     geom_pointrange(aes(xmin = conf.low, xmax = conf.high), size=0.25)