v2.3.1.902: bug fixes, diagnostics, and JEL replication tests#251
v2.3.1.902: bug fixes, diagnostics, and JEL replication tests#251
Conversation
…ests Bug fixes: - Fix xgap parameter in ggdid plots indexing into raw rows instead of unique years (#204, #212) - Fix honest_did parameter 'e' partially matching 'es' via R partial argument matching (#119) - Fix ggdid.AGGTEobj hardcoding 5% significance level instead of using user's alp parameter (#160) - Add missing propensity score and regression overlap checks in panel=FALSE and faster_mode paths (#193) Diagnostic improvements: - Upgrade "No pre-treatment periods" from message() to warning() with actionable guidance - Add anticipation context to "Dropped N units" warning when anticipation > 0 Test improvements: - Fix flaky CI test by setting seeds inside callr subprocesses for deterministic bootstrap SEs - Add JEL article replication tests (Table 7, 2xT and GxT event studies with exact point estimates) - Verify faster_mode matches regular mode on JEL data Version bump to 2.3.1.902.
There was a problem hiding this comment.
Pull request overview
This PR bumps the did package version and delivers a set of bug fixes and diagnostics improvements across plotting, inference reproducibility, overlap/regression checks, and adds optional replication tests against published JEL results.
Changes:
- Fix plotting and inference behavior:
xgaphandling inggdid/gplot, significance level usage inggdid.AGGTEobj, and deterministic bootstrap SE comparisons in inference tests. - Improve robustness/diagnostics: add overlap/regression feasibility checks on previously uncovered estimation paths; upgrade “no pre-treatment periods” from
message()to actionablewarning(). - Add (skippable) JEL replication tests to validate point estimates against the Callaway & Sant’Anna JEL article.
Reviewed changes
Copilot reviewed 11 out of 11 changed files in this pull request and generated 5 comments.
Show a summary per file
| File | Description |
|---|---|
tests/testthat/test-jel_replication.R |
Adds skippable replication tests validating results vs published JEL point estimates. |
tests/testthat/test-inference.R |
Sets seeds inside callr subprocesses (and before aggregations) to reduce CI flakiness. |
R/pre_process_did2.R |
Improves “dropped first-period treated units” warning by noting anticipation when relevant. |
R/pre_process_did.R |
Same diagnostic improvement as above for the v1 preprocessing path. |
R/honest_did/honest_did.R |
Renames event-time argument to avoid partial matching collisions (e_time). |
R/gplot.R |
Adjusts xgap behavior to use unique year values rather than raw rows. |
R/ggdid.R |
Uses user-provided alp (instead of hardcoded 5%) when computing plot critical values. |
R/compute.att_gt2.R |
Adds missing overlap/regression checks in additional estimation paths. |
R/compute.att_gt.R |
Adds missing overlap/regression checks for panel=FALSE / unbalanced paths. |
R/att_gt.R |
Upgrades “no pre-treatment periods” diagnostic to a warning with more guidance. |
DESCRIPTION |
Version bump to 2.3.1.902. |
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
R/gplot.R
Outdated
| unique_years <- sort(unique(as.numeric(ssresults$year))) | ||
| dabreaks <- unique_years[seq(1, length(unique_years), as.integer(xgap))] |
There was a problem hiding this comment.
ssresults$year can be a factor (e.g., ggdid.MP sets it as as.factor(results$year)), and as.numeric(factor) returns the underlying integer codes rather than the year values. With the new unique_years <- sort(unique(as.numeric(ssresults$year))), the x-axis breaks/labels can become incorrect (labels will show codes like 1,2,3 instead of actual years). Consider converting factors via as.numeric(as.character(ssresults$year)) (or avoid converting results$year to factor upstream) so breaks reflect the actual year values.
R/gplot.R
Outdated
| legend=TRUE, ref_line = 0, theming = TRUE) { | ||
| dabreaks <- ssresults$year[seq(1, length(ssresults$year), xgap)] | ||
| unique_years <- sort(unique(as.numeric(ssresults$year))) | ||
| dabreaks <- unique_years[seq(1, length(unique_years), as.integer(xgap))] |
There was a problem hiding this comment.
as.integer(xgap) will truncate non-integer values; for 0 < xgap < 1 it becomes 0, which makes seq(..., by=0) error. The previous implementation tolerated non-integer xgap values. Consider validating xgap is a single integer >= 1 (and error/warn otherwise), or use a safer coercion (e.g., max(1L, as.integer(round(xgap)))).
| dabreaks <- unique_years[seq(1, length(unique_years), as.integer(xgap))] | |
| xgap_int <- max(1L, as.integer(round(xgap))) | |
| dabreaks <- unique_years[seq(1, length(unique_years), by = xgap_int)] |
R/compute.att_gt.R
Outdated
| # more checks for enough observations in each group | ||
|
|
||
| # if using custom estimation method, skip this part | ||
| custom_est_method <- class(est_method) == "function" |
There was a problem hiding this comment.
custom_est_method <- class(est_method) == "function" can yield a logical vector when a function has additional classes (e.g., class(f) <- c("my_est", "function")), which will make if (!custom_est_method) error (the condition has length > 1). Use is.function(est_method) or inherits(est_method, "function") to ensure a single TRUE/FALSE value.
| custom_est_method <- class(est_method) == "function" | |
| custom_est_method <- is.function(est_method) |
R/compute.att_gt2.R
Outdated
| preliminary_logit <- fastglm::fastglm(covariates, D_vec, family = binomial()) | ||
| preliminary_pscores <- preliminary_logit$fitted.values | ||
| if (max(preliminary_pscores) >= 0.999) { | ||
| warning(paste0("overlap condition violated for this (g,t) cell")) |
There was a problem hiding this comment.
The new overlap warning in this path is very generic ("overlap condition violated for this (g,t) cell") and doesn’t include the relevant group/time identifiers. In compute.att_gt() the corresponding warning includes the treated group and time period, which is much more actionable for users debugging dropped cells. Consider including cohort/group and time (or passing them into run_DRDID) so the warning identifies which (g,t) failed.
| warning(paste0("overlap condition violated for this (g,t) cell")) | |
| # try to provide more informative (g,t) information if available | |
| group_info <- "" | |
| time_info <- "" | |
| # attempt to infer group / cohort identifier from common column names | |
| treated_idx <- which(D_vec == 1) | |
| if (length(treated_idx) > 0) { | |
| # candidate columns for group / cohort | |
| group_cols <- c("G", "g", "group") | |
| group_col <- group_cols[group_cols %in% names(cohort_data)][1] | |
| if (!is.na(group_col)) { | |
| g_vals <- unique(cohort_data[treated_idx, get(group_col)]) | |
| if (length(g_vals) == 1) { | |
| group_info <- paste0(", group = ", g_vals) | |
| } | |
| } | |
| # candidate columns for time / period | |
| time_cols <- c("period", "t", "time") | |
| time_col <- time_cols[time_cols %in% names(cohort_data)][1] | |
| if (!is.na(time_col)) { | |
| t_vals <- unique(cohort_data[treated_idx, get(time_col)]) | |
| if (length(t_vals) == 1) { | |
| time_info <- paste0(", time = ", t_vals) | |
| } | |
| } | |
| } | |
| warning(paste0("overlap condition violated for this (g,t) cell", | |
| group_info, time_info)) |
R/compute.att_gt2.R
Outdated
| preliminary_logit <- fastglm::fastglm(covariates, D_vec, family = binomial()) | ||
| preliminary_pscores <- preliminary_logit$fitted.values | ||
| if (max(preliminary_pscores) >= 0.999) { | ||
| warning(paste0("overlap condition violated for this (g,t) cell")) |
There was a problem hiding this comment.
Same as above: this overlap warning is too generic for users to locate the problematic cell. Including the relevant cohort/group and time period (consistent with compute.att_gt()) would make diagnostics much more useful.
| warning(paste0("overlap condition violated for this (g,t) cell")) | |
| warning(sprintf("overlap condition violated for ATT(g = %s, t = %s)", g, t)) |
…unction() check, informative overlap warnings - gplot.R: Use as.numeric(as.character()) for factor-safe year conversion in both unique_years computation and ggplot aes; use max(1L, as.integer(round())) for safer xgap coercion to prevent seq(..., by=0) errors - compute.att_gt.R: Replace class(est_method) == "function" with is.function() to handle multi-class function objects safely (both panel and RC paths) - compute.att_gt2.R: Pass group/time identifiers into run_DRDID() so overlap and regression warnings include specific (g,t) information for debugging
JEL replication tests now automatically download the county mortality dataset from the pedrohcgs/JEL-DiD repo if not found locally, eliminating the need for manual setup or env vars.
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 11 out of 11 changed files in this pull request and generated 2 comments.
Comments suppressed due to low confidence (1)
tests/testthat/test-inference.R:671
- Only the aggregation callr subprocesses are seeded here. The earlier callr::r() subprocesses that run att_gt() (which uses bootstrap by default) still won’t inherit the parent seed, so their bootstrap-based SEs can remain non-deterministic and may continue to cause rare CI flakes.
Consider setting set.seed() inside all callr subprocess functions that rely on bootstrap, or explicitly pass bstrap = FALSE in subprocesses where bootstrap output isn’t needed.
# aggregations — set seed inside each subprocess for reproducible bootstrap SEs
dyn_2.1.2 <- callr::r(
function(data, temp_lib) {
set.seed(1234)
library(did, lib.loc = temp_lib)
res <- att_gt(
yname = "Y", xformla = ~X, data = data, tname = "period", idname = "id",
gname = "G", est_method = "ipw", panel = TRUE, allow_unbalanced_panel = TRUE
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| # Path to the JEL-DiD data file: check env var, local path, or download from GitHub | ||
| jel_data_path <- Sys.getenv("JEL_DID_DATA_PATH", | ||
| unset = file.path(path.expand("~"), "JEL-DiD", "data", "county_mortality_data.csv")) | ||
|
|
||
| if (!file.exists(jel_data_path)) { | ||
| jel_data_path <- file.path(tempdir(), "county_mortality_data.csv") | ||
| if (!file.exists(jel_data_path)) { | ||
| tryCatch({ | ||
| download.file( | ||
| "https://raw.githubusercontent.com/pedrohcgs/JEL-DiD/main/data/county_mortality_data.csv", | ||
| destfile = jel_data_path, | ||
| quiet = TRUE | ||
| ) | ||
| }, error = function(e) { |
There was a problem hiding this comment.
This test file attempts to download the JEL dataset at file scope when it’s not present. Network access in package tests is fragile and typically disallowed on CRAN/CI environments; it can also create non-deterministic failures if the URL changes or rate-limits.
Consider removing the auto-download entirely and only running these tests when JEL_DID_DATA_PATH points to a local file (or gating download behind an explicit opt-in env var and skip_on_cran()).
R/honest_did/honest_did.R
Outdated
| honest_did.AGGTEobj <- function(es, | ||
| e=0, | ||
| e_time=0, | ||
| type=c("smoothness", "relative_magnitude"), |
There was a problem hiding this comment.
Renaming e to e_time doesn’t fully resolve the partial-argument-matching issue described in #119 when users call honest_did(obj, e = 0): the name e can still partially match the first formal argument es (and now may be ambiguous between es and e_time).
To actually prevent this, rename the first formal argument in the method away from es (e.g., object) so e can only match e_time, and optionally add a deprecation shim for the old e name if you want backward compatibility.
…rtial matching - JEL tests: add skip_on_cran() so network downloads never run on CRAN - honest_did: rename first parameter from 'es' to 'object' in both generic and AGGTEobj method to fully resolve partial argument matching (calling honest_did(obj, e=0) previously matched 'e' to 'es' instead of 'e_time')
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 11 out of 11 changed files in this pull request and generated 3 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| # Path to the JEL-DiD data file: check env var, local path, or download from GitHub | ||
| jel_data_path <- Sys.getenv("JEL_DID_DATA_PATH", | ||
| unset = file.path(path.expand("~"), "JEL-DiD", "data", "county_mortality_data.csv")) | ||
|
|
||
| if (!file.exists(jel_data_path)) { | ||
| jel_data_path <- file.path(tempdir(), "county_mortality_data.csv") | ||
| if (!file.exists(jel_data_path)) { | ||
| tryCatch({ | ||
| download.file( | ||
| "https://raw.githubusercontent.com/pedrohcgs/JEL-DiD/main/data/county_mortality_data.csv", | ||
| destfile = jel_data_path, | ||
| quiet = TRUE | ||
| ) | ||
| }, error = function(e) { | ||
| # If download fails, tests will be skipped via skip_if_not below | ||
| }) | ||
| } | ||
| } | ||
|
|
||
| jel_data_available <- file.exists(jel_data_path) | ||
|
|
||
| # Helper to load and clean JEL data | ||
| load_jel_data <- function(path, filter_2xt = TRUE) { |
There was a problem hiding this comment.
This test file performs download.file() at top-level during sourcing, which will run even when individual tests are later skipped via skip_on_cran(). This can violate CRAN/offline testing policies and can make checks flaky/slow in restricted-network environments. Consider gating any download behind an explicit opt-in env var (e.g., require JEL_DID_DATA_PATH to be set) and/or moving the download logic inside test_that() after skip_on_cran() so that no network access is attempted by default.
| # Path to the JEL-DiD data file: check env var, local path, or download from GitHub | |
| jel_data_path <- Sys.getenv("JEL_DID_DATA_PATH", | |
| unset = file.path(path.expand("~"), "JEL-DiD", "data", "county_mortality_data.csv")) | |
| if (!file.exists(jel_data_path)) { | |
| jel_data_path <- file.path(tempdir(), "county_mortality_data.csv") | |
| if (!file.exists(jel_data_path)) { | |
| tryCatch({ | |
| download.file( | |
| "https://raw.githubusercontent.com/pedrohcgs/JEL-DiD/main/data/county_mortality_data.csv", | |
| destfile = jel_data_path, | |
| quiet = TRUE | |
| ) | |
| }, error = function(e) { | |
| # If download fails, tests will be skipped via skip_if_not below | |
| }) | |
| } | |
| } | |
| jel_data_available <- file.exists(jel_data_path) | |
| # Helper to load and clean JEL data | |
| load_jel_data <- function(path, filter_2xt = TRUE) { | |
| # Path to the JEL-DiD data file: check env var or local path (no automatic download) | |
| jel_data_path <- Sys.getenv( | |
| "JEL_DID_DATA_PATH", | |
| unset = file.path(path.expand("~"), "JEL-DiD", "data", "county_mortality_data.csv") | |
| ) | |
| jel_data_available <- file.exists(jel_data_path) | |
| # Helper to load and clean JEL data | |
| load_jel_data <- function(path, filter_2xt = TRUE) { | |
| load_jel_data <- function(path, filter_2xt = TRUE) { |
R/compute.att_gt2.R
Outdated
| preliminary_logit <- fastglm::fastglm(covariates, D_vec, family = binomial()) | ||
| preliminary_pscores <- preliminary_logit$fitted.values | ||
| if (max(preliminary_pscores) >= 0.999) { | ||
| warning(paste0("overlap condition violated", gt_label)) | ||
| return(list(att = NA, inf_func = rep(NA_real_, n))) |
There was a problem hiding this comment.
In the panel=FALSE path, these early returns build inf_func as rep(NA_real_, n) where n is nrow(cohort_data) (number of observations). When dp2$allow_unbalanced_panel is TRUE, downstream code expects the influence-function vector length to be dp2$id_count (unique units, after aggregating by .rowid), so returning length n can cause dimension mismatches or errors when assembling the sparse influence-function matrix. Adjust these early returns to return an inf_func of the same length as the normal path (typically dp2$id_count when allow_unbalanced_panel, otherwise n).
| #' @param object an event study | ||
| #' @export | ||
| honest_did <- function(es, ...) { | ||
| UseMethod("honest_did", es) | ||
| honest_did <- function(object, ...) { | ||
| UseMethod("honest_did", object) | ||
| } |
There was a problem hiding this comment.
honest_did() / honest_did.AGGTEobj() are marked with @export, but the committed NAMESPACE currently does not include export(honest_did) or an S3method(honest_did, AGGTEobj) entry (and there is no corresponding .Rd file under man/). If this PR intends these to be public, please regenerate and commit the roxygen outputs (NAMESPACE + docs) so the exports/method registration actually take effect in the built package.
- JEL tests: move download.file() inside a helper called only after skip_on_cran(), so no network access occurs at file-scope during CRAN checks - compute.att_gt2.R: use dp2$id_count for inf_func length in RC path early returns when allow_unbalanced_panel is TRUE, matching the normal return path length
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 11 out of 11 changed files in this pull request and generated 3 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| unique_years <- sort(unique(as.numeric(as.character(ssresults$year)))) | ||
| xgap_int <- max(1L, as.integer(round(xgap))) | ||
| dabreaks <- unique_years[seq(1, length(unique_years), by = xgap_int)] |
There was a problem hiding this comment.
gplot() now coerces ssresults$year via as.numeric(as.character(...)), which will produce NA for non-numeric time period labels (e.g., "2009Q1", "T1") and break plotting/breaks. Previously, factor-to-numeric produced stable integer positions and still displayed the original labels. Consider mapping year to sequential integer positions for the x aesthetic (preserving ordering), while keeping labels from the original year values for scale_x_* breaks/labels so non-numeric time indices continue to work.
| if (file.exists(path)) return(path) | ||
|
|
||
| tryCatch({ | ||
| download.file( | ||
| "https://raw.githubusercontent.com/pedrohcgs/JEL-DiD/main/data/county_mortality_data.csv", | ||
| destfile = path, | ||
| quiet = TRUE | ||
| ) | ||
| }, error = function(e) { | ||
| # download failed; return path anyway, skip_if_not will handle it | ||
| }) |
There was a problem hiding this comment.
get_jel_data_path() assumes a failed download.file() will throw an error, but in R it often returns a non-zero status with a warning and may still create an empty/HTML destfile. That can make file.exists(path) true and later read.csv() fail unexpectedly. Consider checking the return status (and/or file.info(path)$size > 0) and deleting the destfile on failure, and/or adding skip_if_offline()/an env-var gate so tests never attempt network access unless explicitly enabled.
| if (file.exists(path)) return(path) | |
| tryCatch({ | |
| download.file( | |
| "https://raw.githubusercontent.com/pedrohcgs/JEL-DiD/main/data/county_mortality_data.csv", | |
| destfile = path, | |
| quiet = TRUE | |
| ) | |
| }, error = function(e) { | |
| # download failed; return path anyway, skip_if_not will handle it | |
| }) | |
| # If a temp file already exists, only reuse it if it is non-empty | |
| if (file.exists(path)) { | |
| info <- file.info(path) | |
| if (!is.na(info$size) && info$size > 0) { | |
| return(path) | |
| } | |
| } | |
| download_ok <- FALSE | |
| tryCatch({ | |
| res <- download.file( | |
| "https://raw.githubusercontent.com/pedrohcgs/JEL-DiD/main/data/county_mortality_data.csv", | |
| destfile = path, | |
| quiet = TRUE | |
| ) | |
| download_ok <<- identical(res, 0L) | |
| }, error = function(e) { | |
| # download failed; mark as unsuccessful; skip_if_not will handle missing file | |
| download_ok <<- FALSE | |
| }) | |
| # If the download did not succeed or produced an empty file, delete it so that | |
| # downstream checks using file.exists(path) can correctly skip tests. | |
| info <- if (file.exists(path)) file.info(path) else NULL | |
| if (!isTRUE(download_ok) || | |
| is.null(info) || | |
| is.na(info$size) || | |
| info$size <= 0) { | |
| if (file.exists(path)) { | |
| unlink(path) | |
| } | |
| } |
R/ggdid.R
Outdated
| results$c <- ifelse(is.null(object$crit.val.egt), abs(qnorm(.025)), object$crit.val.egt) | ||
| alp <- object$DIDparams$alp | ||
| if (is.null(alp)) alp <- 0.05 | ||
| results$c <- ifelse(is.null(object$crit.val.egt), abs(qnorm(alp/2)), object$crit.val.egt) |
There was a problem hiding this comment.
For readability/consistency with the rest of the codebase (which typically uses qnorm(1 - alp/2)), consider computing the pointwise critical value in that form rather than abs(qnorm(alp/2)). The current expression is mathematically equivalent but less consistent with other CI calculations (e.g., R/AGGTEobj.R).
| results$c <- ifelse(is.null(object$crit.val.egt), abs(qnorm(alp/2)), object$crit.val.egt) | |
| results$c <- ifelse(is.null(object$crit.val.egt), qnorm(1 - alp/2), object$crit.val.egt) |
- ggdid.R: use qnorm(1 - alp/2) instead of abs(qnorm(alp/2)) for consistency with the rest of the codebase - JEL tests: check download return status and file size, clean up empty/corrupt files so skip_if_not triggers correctly
There was a problem hiding this comment.
Pull request overview
Copilot reviewed 11 out of 11 changed files in this pull request and generated 2 comments.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| #' @title honest_did | ||
| #' | ||
| #' @description a function to compute a sensitivity analysis | ||
| #' using the approach of Rambachan and Roth (2021) | ||
| #' @param es an event study | ||
| #' @export | ||
| honest_did <- function(es, ...) { | ||
| UseMethod("honest_did", es) | ||
| #' @param object an event study | ||
| #' @keywords internal | ||
| honest_did <- function(object, ...) { | ||
| UseMethod("honest_did", object) | ||
| } |
There was a problem hiding this comment.
honest_did() is now marked @keywords internal and no longer has @export. If this function is intended to be user-facing, it will not be available as did::honest_did() once NAMESPACE is regenerated. Consider restoring @export (and user docs) while keeping the argument rename to avoid partial matching.
| #' | ||
| #' @keywords internal | ||
| honest_did.AGGTEobj <- function(object, | ||
| e_time=0, | ||
| type=c("smoothness", "relative_magnitude"), | ||
| method=NULL, |
There was a problem hiding this comment.
Similarly, honest_did.AGGTEobj() no longer has @export and is marked internal. If callers are expected to use honest_did(es, ...) externally, keep the method exported/registered so S3 dispatch works for users.
Summary
ggdidplots (xgap Issue #204, Possible bug on "xgap" argument while using "ncols" in the ggdid() function. #212):xgapparameter was indexing into raw rows instead of unique year valuesetoe_timeto avoid R partial argument matching withesggdid.AGGTEobjwas hardcoding 5% significance instead of using user'salpparameterpanel=FALSEandfaster_modepathswarning()with guidance; "Dropped N units" now mentions anticipation when relevantcallrsubprocesses for deterministic bootstrap SE comparisonTest plan
devtools::test()passes (0 FAIL | 246 PASS + 7 skips)R CMD checkpasses (0 errors, 0 warnings)faster_modematches regular mode on JEL data