Fix estimation, inference and preprocessing edge cases

DESCRIPTION

@@ -14,7 +14,7 @@ Imports:
     Matrix,
     pbapply,
     ggplot2,
-    DRDID (>= 1.3.0),
+    DRDID (>= 1.4.0),
     generics,
     methods,
     tidyr,
@@ -31,7 +31,7 @@ Suggests:
     covr,
     backports,
     broom,
-    testthat (>= 3.0.0),
+    testthat (>= 3.1.8),
     withr,
     remotes,
     callr

NEWS.md

@@ -1,3 +1,42 @@
+# did (development version)
+
+  * Requires `DRDID (>= 1.4.0)` so repeated-cross-section estimation and
+    inference use the corrected upstream implementation.
+  * Correct weighted-estimation feasibility checks: observations with zero
+    normalized sampling weight no longer cause false overlap failures, and
+    regression conditioning uses the sampling-weighted control design used by
+    DRDID. Unweighted checks retain their fast path without weight scans or
+    unnecessary covariate copies.
+  * `aggte(type = "dynamic", balance_e = ...)` now retains only cohorts with
+    non-missing estimates at every observed post-treatment event time through
+    the requested horizon, including on irregular grids. Balancing precedes
+    display limits, and an error explains when no complete cohort remains.
+  * Standard errors, Wald tests, and simultaneous bands no longer discard
+    finite positive variances based on absolute outcome units. Exact-zero and
+    non-finite scales remain unavailable. `att_gt()` reuses the multiplier
+    bootstrap's critical value instead of repeating its quantiles and scaling.
+  * Fast preprocessing accepts used data columns named `args`, `tlist`, or
+    `cutoff_t` without confusing them with local variables.
+  * `conditional_did_pretest()` resolves omitted covariates and observation IDs
+    and uses the actual panel structure after preprocessing. Balanced panels
+    with `allow_unbalanced_panel = TRUE` and repeated cross sections with that
+    ignored flag no longer fail as if they were unbalanced panels.
+  * The legacy conditional pre-test now supports only its justified `"ipw"`
+    method. Regression, doubly robust, and custom methods error clearly: the
+    original outcome model need not remain correct after indicator weighting,
+    so those alternatives did not inherit ordinary outcome-model robustness.
+    `att_gt()` continues to support those estimation methods.
+  * Restore the anticipation correction from PR #276, which had merged into an
+    intermediate branch without reaching master. Anticipation-window cells use
+    the clean base period under either base-period setting and are excluded
+    from the Wald pre-test. `anticipation = 0` and post-treatment estimates are
+    unchanged.
+  * Include the small-group protection and restored tests from PR #278. Built-in
+    `dr` and `reg` repeated-cross-section cells, unbalanced panels, and balanced
+    panels with `fix_weights = "varying"` check normalized positive-weight
+    support in every required outcome regression before native fitting.
+    Undersized designs warn and return an NA cell while healthy cells continue.
+
 # did 2.5.1.902
 
   * Bug fix (`att_gt()`, balanced panel): the never-treated group is now re-checked after the panel is balanced. Balancing drops units whole, so a covariate missing in one period for every never-treated unit removed the whole group after the check had passed; every ATT(g,t) was then `NA` — silently under `faster_mode = TRUE`, and reported as "overlap condition violated" under `faster_mode = FALSE`. The last treated cohort is now used as the comparison group (as when the data contain no never-treated group), with a warning, and the estimates equal those from the balanced sample. Under `control_group = "notyettreated"` the last cohort's pre-treatment cells are no longer computed in this case (it is a comparison group only).

R/MP.R

@@ -20,6 +20,7 @@
 #'  align rows by rowname, never by position.
 #' @param n the number of unique cross-sectional units (unique values of idname)
 #' @param W the Wald statistic for pre-testing the common trends assumption
+#'  (from the pre-treatment cells, t < g - anticipation)
 #' @param Wpval the p-value of the Wald statistic for pre-testing the
 #'  common trends assumption
 #' @param aggte an aggregate treatment effects object

R/aggte.R

@@ -23,8 +23,12 @@
 #'  the sample with respect to event time.  For example, if `balance_e=2`,
 #'  `aggte` will drop groups that are not exposed to treatment for
 #'  at least three periods. (the initial period when `e=0` as well as the
-#'  next two periods when `e=1` and the `e=2`).  This ensures that
-#'  the composition of groups does not change when event time changes.
+#'  next two periods when `e=1` and the `e=2`). Each retained group must
+#'  have a non-missing ATT estimate at every observed post-treatment event
+#'  time through `balance_e`, so the composition of groups does not change
+#'  across those event times. On irregular time grids, only event times
+#'  observed among the eligible groups are required. Balancing is done before
+#'  applying `min_e` and `max_e`; an error is returned if no group remains.
 #' @param min_e The smallest event time (length of exposure, `e = t - g`,
 #'  measured in the same units as `tname`) to include.  By default,
 #'  `min_e = -Inf` so that effects at all lengths of exposure are included.

R/att_gt.R

@@ -26,7 +26,7 @@
 #'  \describe{
 #'    \item{Balanced panel}{Each 2x2 DiD comparison uses the weight from the
 #'      earlier of the two time periods involved. For post-treatment cells,
-#'      this is the base period (g-1). For pre-treatment cells with
+#'      this is the base period (g-anticipation-1). For pre-treatment cells with
 #'      \code{base_period="varying"}, this is the pre-treatment period itself.
 #'      The panel DRDID estimators are used.}
 #'    \item{Repeated cross sections and unbalanced panels}{Both periods'
@@ -44,7 +44,7 @@
 #'  \describe{
 #'    \item{\code{NULL} (default)}{For balanced panel: uses the weight from
 #'      the earlier of the two time periods in each 2x2 comparison. For
-#'      post-treatment cells, this is the base period (g-1). For
+#'      post-treatment cells, this is the base period (g-anticipation-1). For
 #'      pre-treatment cells, this depends on the \code{base_period} setting.
 #'      For RC/unbalanced panel: uses per-observation weights from both
 #'      periods.}
@@ -192,7 +192,10 @@
 #'  date (net of `anticipation`) onward are dropped, with a warning.
 #' @param anticipation The number of time periods before participating
 #'  in the treatment where units can anticipate participating in the
-#'  treatment and therefore it can affect their untreated potential outcomes
+#'  treatment and therefore it can affect their untreated potential outcomes.
+#'  Periods `g - anticipation`, ..., `g - 1` use the base period
+#'  `g - anticipation - 1` (under either `base_period`) and are not used
+#'  in the pre-test of parallel trends.
 #' @param faster_mode This option enables a faster version of `did`, optimizing
 #' computation time for large datasets by improving data management within the package.
 #' The default is set to `TRUE`. Both modes produce identical results up to
@@ -201,7 +204,7 @@
 #' @param base_period Whether to use a "varying" base period or a
 #'  "universal" base period.  Either choice results in the same
 #'  post-treatment estimates of ATT(g,t)'s.  In pre-treatment
-#'  periods, using a varying base period amounts to computing a
+#'  periods (t < g - anticipation), using a varying base period amounts to computing a
 #'  pseudo-ATT in each treatment period by comparing the change
 #'  in outcomes for a particular group relative to its comparison
 #'  group in the pre-treatment periods (i.e., in pre-treatment
@@ -596,8 +599,11 @@ att_gt <- function(yname,
   }
   se <- sqrt(Matrix::diag(V) / n)
 
-  # Zero standard error replaced by NA
-  se[se <= sqrt(.Machine$double.eps) * 10] <- NA
+  # Only unavailable or nonpositive standard errors are degenerate. An
+  # absolute tolerance would change inference when outcome units change.
+  # Magnitude alone cannot distinguish a small valid IF from nuisance-fit
+  # roundoff; numerical singularity is checked using rcond(preV) below.
+  se[!is.finite(se) | se <= 0] <- NA_real_
 
   # If a cluster variable beyond idname is requested without the bootstrap but the cluster-robust variance
   # could not be formed (e.g. cluster identifiers unavailable), fall back to i.i.d. and let the user know.
@@ -611,16 +617,14 @@ att_gt <- function(yname,
   # bootstrap variance matrix
   if (bstrap) {
     bout <- mboot(inffunc, DIDparams = dp, pl = pl, cores = cores, return_V = FALSE)
-    bres <- bout$bres
 
     if (length(zero_na_sd_entry) > 0) {
       se[-zero_na_sd_entry] <- bout$se[-zero_na_sd_entry]
     } else {
       se <- bout$se
     }
   }
-  # Zero standard error replaced by NA
-  se[se <= sqrt(.Machine$double.eps) * 10] <- NA
+  se[!is.finite(se) | se <= 0] <- NA_real_
 
 
   #-----------------------------------------------------------------------------
@@ -652,12 +656,20 @@ att_gt <- function(yname,
   }
 
   if (is.null(wald_invalid)) {
-    # select which periods are pre-treatment
-    pre <- which(group > tt)
+    # select which periods are pre-treatment (anticipation periods are not)
+    pre <- which(group - anticipation > tt)
     # number of pre-treatment cells before the variance filter below, so the
     # no-pre-periods warning can distinguish "no pre-treatment cells exist"
     # from "all pre-treatment cells were dropped for NA/zero variance"
-    n_pre_cells <- length(pre)
+    # (the universal base cell is normalized to 0 and does not count)
+    if (base_period == "universal" && length(pre) > 0L) {
+      time_grid <- if (faster_mode) dp$time_periods else dp$tlist
+      base_times <- time_grid[findInterval(group[pre] - anticipation, time_grid,
+                                           left.open = TRUE)]
+      n_pre_cells <- sum(tt[pre] != base_times)
+    } else {
+      n_pre_cells <- length(pre)
+    }
 
     # Drop group-periods that have variance equal to zero (singularity problems)
     if (length(zero_na_sd_entry) > 0) {
@@ -722,35 +734,14 @@ att_gt <- function(yname,
   # HAVE to use the bootstrap
   if (bstrap) {
     if (cband) {
-      # for uniform confidence band
-      # compute new critical value
-      # see paper for details
-      # Vectorized like mboot(): both quantiles of each bootstrap column are
-      # taken in a single type = 1 call (one sort per column instead of two),
-      # with the qnorm normalization hoisted out of the per-column loop.
-      iqr_norm <- qnorm(.75) - qnorm(.25)
-      bSigma <- apply(bres, 2, function(b) {
-        qs <- quantile(b, c(.25, .75), type = 1, na.rm = TRUE)
-        (qs[2L] - qs[1L]) / iqr_norm
-      })
-
-      bSigma[bSigma <= sqrt(.Machine$double.eps) * 10] <- NA
-
-      # sup-t confidence band: per bootstrap draw, the max over columns of
-      # |b / bSigma|, as a pmax column-fold (equivalent to the previous row-wise
-      # apply(max, na.rm = TRUE): NA entries are ignored via -Inf, and an all-NA
-      # row yields -Inf, exactly as before).
-      if (ncol(bres) == 0L) {
-        # every bootstrap dimension was degenerate (mboot dropped all columns):
-        # the previous row-wise apply(max, na.rm = TRUE) returned -Inf per row
-        bT <- rep(-Inf, nrow(bres))
-      } else {
-        scaled <- abs(bres / rep(bSigma, each = nrow(bres)))
-        scaled[is.na(scaled)] <- -Inf
-        bT <- scaled[, 1L]
-        for (j in seq_len(ncol(scaled))[-1L]) bT <- pmax(bT, scaled[, j])
+      # mboot already computed the studentized sup-t critical value from
+      # these draws. Reuse it rather than sorting and scaling a second time.
+      cval <- bout$crit.val
+      if (!is.finite(cval)) {
+        warning("Simultaneous critical value is not finite, likely because the standard errors could not be computed. Falling back to pointwise confidence intervals.")
+        cval <- qnorm(1 - alp / 2)
+        dp$cband <- FALSE
       }
-      cval <- quantile(bT, 1 - alp, type = 1, na.rm = TRUE)
       if (cval >= 7) {
         warning("Simultaneous critical value is arguably `too large' to be reliable. This usually happens when the number of observations per group is small and/or there is not much variation in outcomes.")
       }

R/compute.aggte.R

@@ -155,6 +155,12 @@ compute.aggte <- function(MP,
     stop("No post-treatment ATT(g,t) estimates are available to aggregate (all cells have t < g).")
   }
 
+  # Keep the observed cell support before na.rm removes estimates. An event
+  # time with all-NA estimates is still required for balancing; an event time
+  # absent from the original grid is not.
+  full_group <- group
+  full_t <- t
+
   if (na.rm) {
     notna <- !is.na(att)
     if (!any(notna)) {
@@ -264,7 +270,6 @@ compute.aggte <- function(MP,
   group <- orig2t_vec(originalgroup)
   glist <- orig2t_vec(originalglist)
   tlist <- unique(t)
-  maxT <- max(t)
 
   # Set the weights
   if (".w" %in% names(dta)) {
@@ -342,9 +347,6 @@ compute.aggte <- function(MP,
 
     # get standard errors from overall influence function
     simple.se <- getSE(simple.if, dp)
-    if (!is.na(simple.se)) {
-      if (simple.se <= sqrt(.Machine$double.eps) * 10) simple.se <- NA
-    }
 
 
     return(AGGTEobj(
@@ -408,7 +410,6 @@ compute.aggte <- function(MP,
 
     # recover standard errors separately by group
     selective.se.g <- unlist(BMisc::get_list_element(selective.se.inner, "se"))
-    selective.se.g[selective.se.g <= sqrt(.Machine$double.eps) * 10] <- NA
 
     # recover influence function separately by group
     selective.inf.func.g <- simplify2array(BMisc::get_list_element(selective.se.inner, "inf.func"))
@@ -465,9 +466,6 @@ compute.aggte <- function(MP,
     selective.inf.func <- as.numeric(selective.inf.func)
     # get overall standard error
     selective.se <- getSE(selective.inf.func, dp)
-    if (!is.na(selective.se)) {
-      if ((selective.se <= sqrt(.Machine$double.eps) * 10)) selective.se <- NA
-    }
 
     return(AGGTEobj(
       overall.att = selective.att,
@@ -508,12 +506,44 @@ compute.aggte <- function(MP,
 
     # if we balance the sample with resepect to event time
     if (!is.null(balance_e)) {
-      include.balanced.gt <- (t2orig(maxT) - originalgroup >= balance_e)
+      # Calendar reach alone does not guarantee that a cohort contributes at
+      # every event time: irregularly spaced periods can leave holes, and
+      # na.rm can remove individual ATT(g,t) cells. Require common observed
+      # post-treatment support through balance_e before applying min_e/max_e.
+      # Use the event times actually represented by the eligible cohorts so
+      # globally unobserved times do not rule out every cohort on sparse grids.
+      event_time <- originalt - originalgroup
+      full_event_time <- full_t - full_group
+      balance_max_t <- max(full_t)
+      eligible_full <- balance_max_t - full_group >= balance_e
+      eligible_g <- unique(full_group[eligible_full])
+      required_e <- unique(full_event_time[eligible_full &
+                                            full_event_time >= 0 & full_event_time <= balance_e])
+      n_required_e <- length(required_e)
+      complete <- logical(length(eligible_g))
+      if (n_required_e > 0L) {
+        # Count distinct (group, event-time) pairs in one pass rather than
+        # rescanning every cell for each group. Integer indices encode pairs
+        # without depending on the magnitude or spacing of the original times.
+        group_index <- match(originalgroup, eligible_g, nomatch = 0L)
+        event_index <- match(event_time, required_e, nomatch = 0L)
+        present <- group_index > 0L & event_index > 0L
+        present_g <- group_index[present]
+        pair_key <- present_g + length(eligible_g) * (as.double(event_index[present]) - 1)
+        counts <- tabulate(present_g[!duplicated(pair_key)], nbins = length(eligible_g))
+        complete <- counts == n_required_e
+      }
+      if (!any(complete)) {
+        stop("No groups have non-missing ATT(g,t) estimates at every observed ",
+             "post-treatment event time through 'balance_e'. Reduce 'balance_e' ",
+             "or inspect the missing group-time estimates.")
+      }
+      include.balanced.gt <- c(FALSE, complete)[group_index + 1L]
 
       eseq <- unique(originalt[include.balanced.gt] - originalgroup[include.balanced.gt])
       eseq <- eseq[order(eseq)]
 
-      eseq <- eseq[(eseq <= balance_e) & (eseq >= balance_e - t2orig(maxT) + t2orig(1))]
+      eseq <- eseq[(eseq <= balance_e) & (eseq >= balance_e - balance_max_t + t2orig(1))]
     }
 
     # only looks at some event times
@@ -555,7 +585,6 @@ compute.aggte <- function(MP,
     })
 
     dynamic.se.e <- unlist(BMisc::get_list_element(dynamic.se.inner, "se"))
-    dynamic.se.e[dynamic.se.e <= sqrt(.Machine$double.eps) * 10] <- NA
 
     dynamic.inf.func.e <- simplify2array(BMisc::get_list_element(dynamic.se.inner, "inf.func"))
 
@@ -597,9 +626,6 @@ compute.aggte <- function(MP,
 
     dynamic.inf.func <- as.numeric(dynamic.inf.func)
     dynamic.se <- getSE(dynamic.inf.func, dp)
-    if (!is.na(dynamic.se)) {
-      if (dynamic.se <= sqrt(.Machine$double.eps) * 10) dynamic.se <- NA
-    }
 
     return(AGGTEobj(
       overall.att = dynamic.att,
@@ -686,7 +712,6 @@ compute.aggte <- function(MP,
 
     # recover standard errors separately by time
     calendar.se.t <- unlist(BMisc::get_list_element(calendar.se.inner, "se"))
-    calendar.se.t[calendar.se.t <= sqrt(.Machine$double.eps) * 10] <- NA
     # recover influence function separately by time
     calendar.inf.func.t <- simplify2array(BMisc::get_list_element(calendar.se.inner, "inf.func"))
 
@@ -731,9 +756,6 @@ compute.aggte <- function(MP,
     calendar.inf.func <- as.numeric(calendar.inf.func)
     # get overall standard error
     calendar.se <- getSE(calendar.inf.func, dp)
-    if (!is.na(calendar.se)) {
-      if (calendar.se <= sqrt(.Machine$double.eps) * 10) calendar.se <- NA
-    }
     return(AGGTEobj(
       overall.att = calendar.att,
       overall.se = calendar.se,
@@ -873,7 +895,7 @@ getSE <- function(thisinffunc, DIDparams = NULL) {
 
   if (bstrap) {
     bout <- mboot(thisinffunc, DIDparams, return_V = FALSE)
-    return(bout$se)
+    se <- bout$se
   } else {
     # Analytical cluster-robust SE (no bootstrap): cluster sums of the aggregated influence function,
     # matching the cluster-sum aggregation of the multiplier bootstrap (CS 2021, Remark 10). Falls back to
@@ -883,8 +905,12 @@ getSE <- function(thisinffunc, DIDparams = NULL) {
     extra_cl <- if (!is.null(DIDparams)) DIDparams$clustervars[!(DIDparams$clustervars %in% c(DIDparams$idname, ""))] else character(0)
     if (length(extra_cl) > 0 && !is.null(cv) && length(cv) == length(thisinffunc)) {
       S <- rowsum(thisinffunc, cv)
-      return(sqrt(sum(S^2)) / n)
+      se <- sqrt(sum(S^2)) / n
+    } else {
+      se <- sqrt(mean((thisinffunc)^2) / n)
     }
-    return(sqrt(mean((thisinffunc)^2) / n))
   }
+  # No absolute cutoff: the same positive variance must remain available in
+  # every outcome unit. Zero and non-finite scales cannot support inference.
+  if (!is.finite(se) || se <= 0) NA_real_ else se
 }