v2.3.1.902: bug fixes, new broom methods, and documentation improvements

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: did
 Title: Treatment Effects with Multiple Periods and Groups
-Version: 2.3.1.902
+Version: 2.3.1.903
 Authors@R: c(person("Brantly", "Callaway", email = "brantly.callaway@uga.edu", role = c("aut", "cre")), person("Pedro H. C.", "Sant'Anna", email="pedro.santanna@emory.edu", role = c("aut")))
 URL: https://bcallaway11.github.io/did/, https://github.com/bcallaway11/did/
 Description: The standard Difference-in-Differences (DID) setup involves two periods and two groups -- a treated group and untreated group.  Many applications of DID methods involve more than two periods and have individuals that are treated at different points in time.  This package contains tools for computing average treatment effect parameters in Difference in Differences setups with more than two periods and with variation in treatment timing using the methods developed in Callaway and Sant'Anna (2021) <doi:10.1016/j.jeconom.2020.12.001>.  The main parameters are group-time average treatment effects which are the average treatment effect for a particular group at a a particular time.  These can be aggregated into a fewer number of treatment effect parameters, and the package deals with the cases where there is selective treatment timing, dynamic treatment effects, calendar time effects, or combinations of these.  There are also functions for testing the Difference in Differences assumption, and plotting group-time average treatment effects.
@@ -30,4 +30,5 @@ Suggests:
     here,
     knitr,
     covr,
-    backports
+    backports,
+    broom

NAMESPACE

@@ -4,6 +4,8 @@ S3method(ggdid,AGGTEobj)
 S3method(ggdid,MP)
 S3method(glance,AGGTEobj)
 S3method(glance,MP)
+S3method(nobs,AGGTEobj)
+S3method(nobs,MP)
 S3method(print,AGGTEobj)
 S3method(print,MP)
 S3method(summary,AGGTEobj)
@@ -53,4 +55,5 @@ importFrom(generics,glance)
 importFrom(generics,tidy)
 importFrom(methods,as)
 importFrom(methods,is)
+importFrom(stats,nobs)
 importFrom(tidyr,gather)

R/att_gt.R

@@ -378,46 +378,76 @@ att_gt <- function(yname,
   # compute Wald pre-test
   #-----------------------------------------------------------------------------
 
-  # select which periods are pre-treatment
-  pre <- which(group > tt)
-
-  # Drop group-periods that have variance equal to zero (singularity problems)
-  if (length(zero_na_sd_entry) > 0) {
-    pre <- pre[!(pre %in% zero_na_sd_entry)]
-  }
-  # pseudo-atts in pre-treatment periods
-  preatt <- as.matrix(att[pre])
-
-  # covariance matrix of pre-treatment atts
-  preV <- as.matrix(V[pre, pre])
-
-  # check if there are actually any pre-treatment periods
-  if (length(preV) == 0) {
-    msg <- paste0(
-      "No pre-treatment periods available for the Wald pre-test of parallel trends. ",
-      "This can happen when all groups are first treated early in the panel ",
-      "(e.g., in the second time period) so that no pre-treatment ATT(g,t) estimates exist."
+  # Determine whether the analytical V is a valid basis for the Wald pre-test.
+  # V = t(inffunc) %*% inffunc / n is valid for balanced panels, unbalanced
+  # panels (influence functions are aggregated to the unit level), and repeated
+  # cross-sections (CLT applies per independent observation).
+  #
+  # It is invalid when an extra cluster variable (beyond idname) is specified:
+  # the analytical V ignores between-cluster correlation, making the Wald stat
+  # anti-conservative.  In that case we skip the test and message the user to
+  # rely on bootstrap confidence intervals instead.
+
+  extra_clustervars <- clustervars[!(clustervars %in% c(idname, ""))]
+
+  wald_invalid <- NULL
+  if (length(extra_clustervars) > 0) {
+    wald_invalid <- paste0(
+      "The Wald pre-test is not reported when clustering beyond the unit level ",
+      "(clustervars = '", paste(extra_clustervars, collapse = "', '"), "') ",
+      "because the analytical variance matrix does not account for between-cluster ",
+      "correlation. Use the bootstrap confidence intervals to assess pre-trends."
     )
-    if (anticipation > 0) {
-      msg <- paste0(msg, " Note: anticipation=", anticipation, " further reduces the number of available pre-treatment periods.")
-    }
-    warning(msg)
-    W <- NULL
-    Wpval <- NULL
-  } else if (sum(is.na(preV))) {
-    warning("Not returning pre-test Wald statistic due to NA pre-treatment values")
-    W <- NULL
-    Wpval <- NULL
-  } else if (rcond(preV) <= .Machine$double.eps) {
-    # singular covariance matrix for pre-treatment periods
-    warning("Not returning pre-test Wald statistic due to singular covariance matrix")
+  }
+
+  if (!is.null(wald_invalid)) {
+    message(wald_invalid)
     W <- NULL
     Wpval <- NULL
-  } else {
-    # everything is working...
-    W <- n * t(preatt) %*% solve(preV) %*% preatt
-    q <- length(pre) # number of restrictions
-    Wpval <- round(1 - pchisq(W, q), 5)
+  }
+
+  if (is.null(wald_invalid)) {
+    # select which periods are pre-treatment
+    pre <- which(group > tt)
+
+    # Drop group-periods that have variance equal to zero (singularity problems)
+    if (length(zero_na_sd_entry) > 0) {
+      pre <- pre[!(pre %in% zero_na_sd_entry)]
+    }
+    # pseudo-atts in pre-treatment periods
+    preatt <- as.matrix(att[pre])
+
+    # covariance matrix of pre-treatment atts
+    preV <- as.matrix(V[pre, pre])
+
+    # check if there are actually any pre-treatment periods
+    if (length(preV) == 0) {
+      msg <- paste0(
+        "No pre-treatment periods available for the Wald pre-test of parallel trends. ",
+        "This can happen when all groups are first treated early in the panel ",
+        "(e.g., in the second time period) so that no pre-treatment ATT(g,t) estimates exist."
+      )
+      if (anticipation > 0) {
+        msg <- paste0(msg, " Note: anticipation=", anticipation, " further reduces the number of available pre-treatment periods.")
+      }
+      warning(msg)
+      W <- NULL
+      Wpval <- NULL
+    } else if (sum(is.na(preV))) {
+      warning("Not returning pre-test Wald statistic due to NA pre-treatment values")
+      W <- NULL
+      Wpval <- NULL
+    } else if (rcond(preV) <= .Machine$double.eps) {
+      # singular covariance matrix for pre-treatment periods
+      warning("Not returning pre-test Wald statistic due to singular covariance matrix")
+      W <- NULL
+      Wpval <- NULL
+    } else {
+      # everything is working...
+      W <- n * t(preatt) %*% solve(preV) %*% preatt
+      q <- length(pre) # number of restrictions
+      Wpval <- round(1 - pchisq(W, q), 5)
+    }
   }
 
 

R/pre_process_did.R

@@ -202,22 +202,12 @@ pre_process_did <- function(yname,
     #  make sure id is numeric
     if (! (is.numeric(data[, idname])) ) stop("The id variable '", idname, "' must be numeric. Please convert it.")
 
-    ## # checks below are useful, but removing due to being slow
-    ## # might ought to figure out a way to do these faster later
-    ## # these checks are also closely related to making sure
-    ## # that we have a well-balanced panel, so it might make
-    ## # sense to move them over to the BMisc package
-
-    ## # Check if idname is unique by tname
-    ## n_id_year = all( table(data[, idname], data[, tname]) <= 1)
-    ## if (! n_id_year) stop("The value of idname must be the unique (by tname)")
-
-    ## # make sure gname doesn't change across periods for particular individuals
-    ## if (!all(sapply( split(data, data[,idname]), function(df) {
-    ##   length(unique(df[,gname]))==1
-    ## }))) {
-    ##   stop("The value of gname must be the same across all periods for each particular individual.")
-    ## }
+    # Check that gname is time-invariant within each unit (treatment irreversibility).
+    # Uses unique() + anyDuplicated() for O(n) performance.
+    id_g_unique <- unique(data[, c(idname, gname)])
+    if (anyDuplicated(id_g_unique[[idname]]) != 0L) {
+      stop("The value of gname (treatment variable) must be the same across all periods for each particular unit. The treatment must be irreversible.")
+    }
   }
 
 

R/pre_process_did2.R

@@ -19,6 +19,14 @@ validate_args <- function(args, data){
   name_message <- "__ARG__ must be a character scalar and a name of a column from the dataset."
   dreamerr::check_set_arg(args$tname, args$gname, args$yname, "match", .choices = data_names, .message = name_message, .up = 1)
 
+  # Strip idname from clustervars before validation: users may naturally pass
+  # clustervars = c(idname, extra_var), and idname clustering is implicit/redundant.
+  # This must happen before the dreamerr check which enforces scalar length.
+  if (!is.null(args$clustervars) && !is.null(args$idname) && (args$idname %in% args$clustervars)) {
+    args$clustervars <- setdiff(args$clustervars, args$idname)
+    if (length(args$clustervars) == 0L) args$clustervars <- NULL
+  }
+
   # Flag for clustervars and weightsname
   checkvar_message <- "__ARG__ must be NULL or a character scalar that is a name of a column from the dataset."
   dreamerr::check_set_arg(args$weightsname, args$clustervars, "NULL | match", .choices = data_names, .message = checkvar_message, .up = 1)
@@ -60,13 +68,9 @@ validate_args <- function(args, data){
   dreamerr::check_set_arg(args$base_period, "match", .choices = c("universal", "varying"), .message = base_period_message, .up = 1)
 
   # Flags for cluster variable
+  # Note: idname was already stripped from clustervars above (before dreamerr check)
   if (!is.null(args$clustervars)) {
-    # dropping idname from cluster
-    if ((!is.null(args$idname)) && (args$idname %in% args$clustervars)){
-      args$clustervars <- setdiff(args$clustervars, args$idname)
-    }
-
-    # check if user is providing more than 2 cluster variables (different than idname)
+    # check if user is providing more than 1 cluster variable (beyond idname)
     if (length(args$clustervars) > 1) {
       stop("You can only provide 1 cluster variable additionally to the one provided in idname. Please check your arguments.")
     }
@@ -584,6 +588,14 @@ pre_process_did2 <- function(yname,
   # run error checking on arguments
   validate_args(args, data)
 
+  # Strip idname from clustervars in the caller's copy of args.
+  # validate_args works on a local copy, so we must repeat this here.
+  # Clustering at the unit level (idname) is already done implicitly.
+  if (!is.null(args$clustervars) && !is.null(args$idname) && (args$idname %in% args$clustervars)) {
+    args$clustervars <- setdiff(args$clustervars, args$idname)
+    if (length(args$clustervars) == 0L) args$clustervars <- NULL
+  }
+
   # put in blank xformla if no covariates
   if (is.null(args$xformla)) {
     args$xformla <- ~1

R/tidy.R

@@ -6,11 +6,60 @@ generics::tidy
 #' @export
 generics::glance
 
-#' tidy results from MP objects
+#' Number of unique cross-sectional units in an MP object
+#'
+#' @importFrom stats nobs
+#' @param object a model of class MP produced by the [att_gt()] function
+#' @param ... additional arguments (ignored)
+#' @return Integer. The number of unique cross-sectional units in the data.
+#' @export
+nobs.MP <- function(object, ...) {
+  as.integer(object$n)
+}
+
+#' Number of unique cross-sectional units in an AGGTEobj object
+#'
+#' @importFrom stats nobs
+#' @param object a model of class AGGTEobj produced by the [aggte()] function
+#' @param ... additional arguments (ignored)
+#' @return Integer. The number of unique cross-sectional units in the data.
+#' @export
+nobs.AGGTEobj <- function(object, ...) {
+  if (object$DIDparams$faster_mode) {
+    as.integer(object$DIDparams$id_count)
+  } else {
+    as.integer(object$DIDparams$n)
+  }
+}
+
+#' Tidy an MP object into a data frame
+#'
+#' Returns a tidy data frame of group-time average treatment effect estimates
+#' from an [att_gt()] result.
 #'
 #' @importFrom generics tidy
 #' @param x a model of class MP produced by the [att_gt()] function
 #' @param ... Additional arguments to tidying method.
+#'
+#' @return A data frame with one row per ATT(g,t) estimate and columns:
+#' \describe{
+#'   \item{term}{ATT(g,t) label}
+#'   \item{group}{the treatment cohort g}
+#'   \item{time}{the time period t}
+#'   \item{estimate}{the ATT(g,t) point estimate}
+#'   \item{std.error}{standard error}
+#'   \item{statistic}{t-statistic (\code{estimate / std.error})}
+#'   \item{p.value}{two-sided pointwise p-value
+#'     (\code{2 * (1 - pnorm(abs(statistic)))}).
+#'     Marginal per-estimate; does \strong{not} account for multiple testing
+#'     across ATT(g,t) cells.}
+#'   \item{conf.low, conf.high}{simultaneous confidence band limits, using the
+#'     bootstrap uniform critical value when \code{bstrap=TRUE} and
+#'     \code{cband=TRUE}, otherwise pointwise}
+#'   \item{point.conf.low, point.conf.high}{pointwise confidence interval limits
+#'     using \code{qnorm(1 - alp/2)}, always}
+#' }
+#'
 #' @export
 tidy.MP <- function(x, ...) {
   out <- data.frame(
@@ -19,6 +68,8 @@ tidy.MP <- function(x, ...) {
     time      = x$t,
     estimate  = x$att,
     std.error = x$se,
+    statistic = x$att / x$se,
+    p.value   = 2 * (1 - stats::pnorm(abs(x$att / x$se))),
     conf.low  = x$att - x$c * x$se,
     conf.high = x$att + x$c * x$se,
     point.conf.low  = x$att - stats::qnorm(1 - x$alp/2) * x$se,
@@ -42,11 +93,44 @@ glance.MP <- function(x, ...) {
   out
 }
 
-#' tidy results from AGGTEobj objects
+#' Tidy an AGGTEobj into a data frame
+#'
+#' Returns a tidy data frame of aggregated treatment effect estimates from an
+#' [aggte()] result.
 #'
 #' @importFrom generics tidy
 #' @param x a model of class AGGTEobj produced by the [aggte()] function
 #' @param ... Additional arguments to tidying method.
+#'
+#' @return A data frame whose columns depend on \code{type}:
+#' \describe{
+#'   \item{type}{the aggregation type: \code{"simple"}, \code{"dynamic"},
+#'     \code{"group"}, or \code{"calendar"}}
+#'   \item{term}{label for each estimate}
+#'   \item{estimate}{point estimate}
+#'   \item{std.error}{standard error}
+#'   \item{statistic}{t-statistic (\code{estimate / std.error})}
+#'   \item{p.value}{two-sided pointwise p-value
+#'     (\code{2 * (1 - pnorm(abs(statistic)))}).
+#'     Marginal per-estimate; does \strong{not} account for multiple testing
+#'     across event times or groups.}
+#'   \item{conf.low, conf.high}{simultaneous confidence band limits.  When
+#'     \code{bstrap=TRUE} and \code{cband=TRUE} these use the bootstrap uniform
+#'     critical value (\code{crit.val.egt}); otherwise they equal the pointwise
+#'     intervals.  For \code{type="simple"} and the overall average row of
+#'     \code{type="group"}, a single scalar is returned so simultaneous and
+#'     pointwise coincide.}
+#'   \item{point.conf.low, point.conf.high}{pointwise confidence interval limits
+#'     always using \code{qnorm(1 - alp/2)}.}
+#' }
+#'
+#' @details
+#' The key distinction between \code{conf.low}/\code{conf.high} and
+#' \code{point.conf.low}/\code{point.conf.high} is that the former accounts for
+#' multiple testing across all estimates (simultaneous coverage), while the
+#' latter provides marginal (per-estimate) coverage only.  Use the simultaneous
+#' bands when you want to make joint inferences across all event times or groups.
+#'
 #' @export
 tidy.AGGTEobj<- function(x, ...) {
   if(x$type == "dynamic"){
@@ -56,6 +140,8 @@ tidy.AGGTEobj<- function(x, ...) {
       event.time= x$egt,
       estimate  = x$att.egt,
       std.error = x$se.egt,
+      statistic = x$att.egt / x$se.egt,
+      p.value   = 2 * (1 - stats::pnorm(abs(x$att.egt / x$se.egt))),
       conf.low  = x$att.egt - x$crit.val.egt * x$se.egt,
       conf.high = x$att.egt + x$crit.val.egt * x$se.egt,
       point.conf.low  = x$att.egt - stats::qnorm(1 - x$DIDparams$alp/2) * x$se.egt,
@@ -68,6 +154,8 @@ tidy.AGGTEobj<- function(x, ...) {
       group    = c('Average', x$egt),
       estimate  = c(x$overall.att, x$att.egt),
       std.error = c(x$overall.se, x$se.egt),
+      statistic = c(x$overall.att, x$att.egt) / c(x$overall.se, x$se.egt),
+      p.value   = 2 * (1 - stats::pnorm(abs(c(x$overall.att, x$att.egt) / c(x$overall.se, x$se.egt)))),
       conf.low  = c(x$overall.att - stats::qnorm(1 - x$DIDparams$alp/2) * x$overall.se, x$att.egt - x$crit.val.egt * x$se.egt),
       conf.high = c(x$overall.att + stats::qnorm(1 - x$DIDparams$alp/2) * x$overall.se, x$att.egt + x$crit.val.egt * x$se.egt),
       point.conf.low  = c(x$overall.att - stats::qnorm(1 - x$DIDparams$alp/2) * x$overall.se, x$att.egt - stats::qnorm(1 - x$DIDparams$alp/2) * x$se.egt),
@@ -81,18 +169,22 @@ tidy.AGGTEobj<- function(x, ...) {
       term = paste0('ATT(', x$egt, ")"),
       estimate  = x$att.egt,
       std.error = x$se.egt,
+      statistic = x$att.egt / x$se.egt,
+      p.value   = 2 * (1 - stats::pnorm(abs(x$att.egt / x$se.egt))),
       conf.low  = x$att.egt - x$crit.val.egt * x$se.egt,
       conf.high = x$att.egt + x$crit.val.egt * x$se.egt,
       point.conf.low  = x$att.egt - stats::qnorm(1 - x$DIDparams$alp/2) * x$se.egt,
       point.conf.high = x$att.egt + stats::qnorm(1 - x$DIDparams$alp/2) * x$se.egt)
   }
-  
+
   if(x$type == "simple"){
     out <- data.frame(
       type      = x$type,
       term      = 'ATT(simple average)',
       estimate  = x$overall.att,
       std.error = x$overall.se,
+      statistic = x$overall.att / x$overall.se,
+      p.value   = 2 * (1 - stats::pnorm(abs(x$overall.att / x$overall.se))),
       conf.low  = x$overall.att - stats::qnorm(1 - x$DIDparams$alp/2) * x$overall.se,
       conf.high = x$overall.att + stats::qnorm(1 - x$DIDparams$alp/2) * x$overall.se,
       point.conf.low  = x$overall.att - stats::qnorm(1 - x$DIDparams$alp/2) * x$overall.se,