Skip to content

Commit 8b10dcf

Browse files
committed
ensembleHTE 0.3.0: decouple SE clustering from fold splitting via se_cluster_id
- Add se_cluster_id to ensemble_hte() and ensemble_pred() so cluster-robust SEs can be computed at a different level than the fold-splitting unit. - If only one of individual_id / se_cluster_id is given, it is used for both. - Print an early message reporting the split level and clustering level. - Store $se_cluster_id and name labels; print shows both levels separately. - Downstream analyses read se_cluster_id with fallback to individual_id. - Fix: cluster the inner cross-validated ensemble split by the panel id too.
1 parent a35204f commit 8b10dcf

12 files changed

Lines changed: 393 additions & 106 deletions

DESCRIPTION

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
Package: ensembleHTE
22
Type: Package
33
Title: Ensemble Methods for Learning Features of Heterogeneous Treatment Effects
4-
Version: 0.2.1
4+
Version: 0.3.0
55
Authors@R: c(
66
person("Bruno", "Fava", email = "brunovnfava@gmail.com", role = c("aut", "cre"),
77
comment = c(ORCID = "0009-0002-8218-516X"))

NEWS.md

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,31 @@
1+
# ensembleHTE 0.3.0
2+
3+
## New features
4+
5+
* `ensemble_hte()` and `ensemble_pred()` gain a `se_cluster_id` argument that
6+
decouples the level of cluster-robust standard errors from the level used to
7+
form cross-fitting folds. `individual_id` now controls only how folds are
8+
built (all of a unit's rows stay in one fold), while `se_cluster_id` controls
9+
the clustering of standard errors in the downstream analyses (`blp()`,
10+
`gates()`, `clan()`, `gavs()`, ...). If only one of the two is supplied, it is
11+
used for both roles. When a fit starts, an informative message reports the
12+
fold-splitting level and the SE-clustering level. Fit objects now expose
13+
`$se_cluster_id` (and the `$individual_id_name` / `$se_cluster_id_name`
14+
labels), and `print()` shows the two levels separately. A typical use is
15+
predicting outcomes for unobserved units within observed clusters (split by
16+
individual, cluster SEs by village).
17+
18+
## Bug fixes
19+
20+
* Panel data (`individual_id`) is now respected in the cross-validated
21+
ensemble step, not just the outer cross-fitting split. Previously, when
22+
`ensemble_strategy = "cv"`, the inner ensemble-weight CV split ignored the
23+
cluster identifier, so observations from the same individual could land in
24+
different ensemble folds and appear on both the fitting and prediction sides
25+
of that inner CV. The ensemble split (and the baseline ensemble that reuses
26+
it) in `ensemble_hte()` and `ensemble_pred()` now cluster by
27+
`individual_id`, matching the outer split.
28+
129
# ensembleHTE 0.2.1
230

331
## Bug fixes

R/analysis_hte.R

Lines changed: 9 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -665,8 +665,8 @@ gates <- function(ensemble_fit, n_groups = 3, outcome = NULL, treatment = NULL,
665665
check_strata <- if (!is.null(strata_vec)) strata_vec[check_ref] else NULL
666666
.check_small_cells(check_folds, n_groups, restrict_by = check_strata, func_name = "GATES")
667667

668-
# Extract individual_id for cluster-robust SEs (if panel data)
669-
cluster_id <- ensemble_fit$individual_id
668+
# Extract the SE-clustering identifier (if panel/clustered data)
669+
cluster_id <- .fit_cluster_id(ensemble_fit)
670670

671671
gates_by_rep <- lapply(1:M, function(m) {
672672
eff_controls <- controls
@@ -1127,8 +1127,9 @@ blp <- function(ensemble_fit, outcome = NULL, treatment = NULL,
11271127
# Extract components from ensemble_fit
11281128
M <- ensemble_fit$M
11291129

1130-
# Extract individual_id for cluster-robust SEs (if panel data)
1131-
cluster_id <- if (!is.null(ensemble_fit$individual_id)) ensemble_fit$individual_id[use_idx] else NULL
1130+
# Extract the SE-clustering identifier (if panel/clustered data)
1131+
.cid <- .fit_cluster_id(ensemble_fit)
1132+
cluster_id <- if (!is.null(.cid)) .cid[use_idx] else NULL
11321133

11331134
# Resolve baseline_as_control
11341135
has_stored_baseline <- inherits(ensemble_fit, "ensemble_hte_fit") &&
@@ -1906,8 +1907,8 @@ clan <- function(ensemble_fit, variables = NULL, n_groups = 3, na_rm = FALSE, sc
19061907
check_folds <- lapply(splits, function(s) s[check_ref])
19071908
.check_small_cells(check_folds, n_groups, func_name = "CLAN")
19081909

1909-
# Extract individual_id for cluster-robust SEs (if panel data)
1910-
cluster_id <- ensemble_fit$individual_id
1910+
# Extract the SE-clustering identifier (if panel/clustered data)
1911+
cluster_id <- .fit_cluster_id(ensemble_fit)
19111912

19121913
# Check for NAs in analysis variables and warn if na_rm = FALSE
19131914
if (!na_rm) {
@@ -2385,8 +2386,8 @@ gavs <- function(ensemble_fit, n_groups = 3, outcome = NULL, subset = NULL,
23852386
check_strata <- if (!is.null(strata_vec)) strata_vec[check_ref] else NULL
23862387
.check_small_cells(check_folds, n_groups, restrict_by = check_strata, func_name = "GAVS")
23872388

2388-
# Extract individual_id for cluster-robust SEs (if panel data)
2389-
cluster_id <- ensemble_fit$individual_id
2389+
# Extract the SE-clustering identifier (if panel/clustered data)
2390+
cluster_id <- .fit_cluster_id(ensemble_fit)
23902391

23912392
gavs_by_rep <- lapply(1:M, function(m) {
23922393
.gavs_single(

R/analysis_pred.R

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -230,8 +230,9 @@ blp_pred <- function(ensemble_fit, outcome = NULL, subset = NULL) {
230230
# Extract components from ensemble_fit
231231
M <- ensemble_fit$M
232232

233-
# Extract individual_id for cluster-robust SEs (if panel data)
234-
cluster_id <- if (!is.null(ensemble_fit$individual_id)) ensemble_fit$individual_id[use_idx] else NULL
233+
# Extract the SE-clustering identifier (if panel/clustered data)
234+
.cid <- .fit_cluster_id(ensemble_fit)
235+
cluster_id <- if (!is.null(.cid)) .cid[use_idx] else NULL
235236

236237
# Compute BLP for each repetition
237238
blp_by_rep <- lapply(1:M, function(m) {

R/compare_restricted.R

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -580,8 +580,8 @@ gavs_restricted <- function(ensemble_fit, restrict_by, n_groups = 3, outcome = N
580580
.check_small_cells(check_folds, n_groups, restrict_by = NULL, func_name = "GAVS Comparison (unrestricted)")
581581
.check_small_cells(check_folds, n_groups, restrict_by = check_strata_rest, func_name = "GAVS Comparison (restricted)")
582582

583-
# Extract individual_id for cluster-robust SEs (if panel data)
584-
cluster_id <- ensemble_fit$individual_id
583+
# Extract the SE-clustering identifier (if panel/clustered data)
584+
cluster_id <- .fit_cluster_id(ensemble_fit)
585585

586586
results_by_rep <- lapply(1:M, function(m) {
587587
gavs_unrest <- .gavs_single(
@@ -1333,8 +1333,8 @@ gates_restricted <- function(ensemble_fit, restrict_by, n_groups = 3, outcome =
13331333
.check_small_cells(check_folds, n_groups, restrict_by = NULL, func_name = "GATES Comparison (unrestricted)")
13341334
.check_small_cells(check_folds, n_groups, restrict_by = check_strata_rest, func_name = "GATES Comparison (restricted)")
13351335

1336-
# Extract individual_id for cluster-robust SEs (if panel data)
1337-
cluster_id <- ensemble_fit$individual_id
1336+
# Extract the SE-clustering identifier (if panel/clustered data)
1337+
cluster_id <- .fit_cluster_id(ensemble_fit)
13381338

13391339
# Compute GATES for each repetition (both unrestricted and restricted)
13401340
results_by_rep <- lapply(1:M, function(m) {

R/ensemble_hte.R

Lines changed: 88 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -229,15 +229,32 @@
229229
#' See \strong{Ensemble Strategy} section for details.
230230
#' @param individual_id Required when the dataset is a panel (e.g., individuals
231231
#' observed over multiple time periods). Specifies the column that identifies
232-
#' individuals so that (1) all observations for the same individual are placed
233-
#' in the same cross-fitting fold, and (2) cluster-robust standard errors are
234-
#' used in all downstream analyses.
232+
#' the unit used to build cross-fitting folds, so that all observations for
233+
#' the same unit are placed in the same fold. By default this identifier is
234+
#' also used to cluster the standard errors in all downstream analyses; supply
235+
#' \code{se_cluster_id} to cluster at a different level (see below).
235236
#'
236237
#' Example: for a panel of students observed across semesters, set
237238
#' \code{individual_id = student_id}.
238239
#'
239240
#' Can be an unquoted column name, a quoted string (\code{"student_id"}),
240241
#' or a vector of identifiers.
242+
#' @param se_cluster_id Optional identifier for the level at which cluster-robust
243+
#' standard errors are computed in the downstream analyses (\code{\link{blp}},
244+
#' \code{\link{gates}}, \code{\link{clan}}, ...). Decouples SE clustering from
245+
#' fold splitting: \code{individual_id} controls how folds are formed, while
246+
#' \code{se_cluster_id} controls SE clustering. If only one of
247+
#' \code{individual_id} / \code{se_cluster_id} is supplied, that identifier is
248+
#' used for both roles. A message reporting the split level and the clustering
249+
#' level is printed when the fit starts.
250+
#'
251+
#' A typical use is predicting outcomes for unobserved units within observed
252+
#' clusters (e.g., new villagers in sampled villages): split at the individual
253+
#' level (\code{individual_id = id}) but cluster SEs at the village level
254+
#' (\code{se_cluster_id = village}).
255+
#'
256+
#' Accepts an unquoted column name, a quoted string, or a vector of
257+
#' identifiers. Defaults to \code{individual_id}.
241258
#' @param n_cores Integer. Number of cores for parallel processing of repetitions.
242259
#' Default is 1 (sequential). Set to higher values to parallelize the M repetitions.
243260
#' Uses the \code{future} framework, so users can also set up their own parallel
@@ -297,7 +314,9 @@
297314
#' \item{task_type}{Task type (regr or classif)}
298315
#' \item{scale_covariates}{Whether covariates were scaled}
299316
#' \item{tune, tune_params}{Tuning settings}
300-
#' \item{individual_id}{Vector of individual identifiers (if panel data)}
317+
#' \item{individual_id}{Vector of fold-splitting identifiers (if panel data)}
318+
#' \item{se_cluster_id}{Vector of SE-clustering identifiers (if supplied or
319+
#' inherited from \code{individual_id})}
301320
#' \item{n_cores}{Number of cores used for parallel processing}
302321
#' }
303322
#'
@@ -432,6 +451,7 @@ ensemble_hte <- function(formula = NULL, treatment = NULL, data = NULL,
432451
train_idx = NULL,
433452
ensemble_strategy = c("cv", "average"),
434453
individual_id = NULL,
454+
se_cluster_id = NULL,
435455
n_cores = 1,
436456
store_baseline = c("ensemble", "none", "all")) {
437457

@@ -763,42 +783,57 @@ ensemble_hte <- function(formula = NULL, treatment = NULL, data = NULL,
763783
stop("Y must not have NA values for training observations (where train_idx = TRUE)")
764784
}
765785

766-
# Handle individual_id for panel data
767-
individual_id_vec <- NULL
786+
# Handle panel / cluster identifiers.
787+
# `individual_id` governs fold splitting (all its observations stay together
788+
# in one cross-fitting fold); `se_cluster_id` governs cluster-robust standard
789+
# errors in the downstream analyses. If only one is supplied, it is used for
790+
# BOTH roles.
791+
id_ind <- NULL
768792
if (!is.null(individual_id)) {
769-
# Try to resolve individual_id as column name (like treatment)
770-
id_expr <- substitute(individual_id)
771-
id_resolved <- tryCatch(
772-
parse_column_name(id_expr, parent.frame(), "individual_id", data),
773-
error = function(e) NULL
774-
)
775-
776-
if (!is.null(id_resolved) && id_resolved %in% names(data)) {
777-
individual_id_vec <- data[[id_resolved]]
778-
} else if (is.character(individual_id) && length(individual_id) == 1 && individual_id %in% names(data)) {
779-
individual_id_vec <- data[[individual_id]]
780-
} else if (length(individual_id) == n) {
781-
individual_id_vec <- individual_id
782-
} else {
783-
stop("individual_id must be a column name in the data or a vector of length n (", n, ")")
793+
id_ind <- .resolve_panel_id(individual_id, substitute(individual_id),
794+
parent.frame(), data, n, "individual_id")
795+
}
796+
id_se <- NULL
797+
if (!is.null(se_cluster_id)) {
798+
id_se <- .resolve_panel_id(se_cluster_id, substitute(se_cluster_id),
799+
parent.frame(), data, n, "se_cluster_id")
800+
}
801+
802+
split_id_vec <- NULL; split_id_name <- NULL
803+
cluster_id_vec <- NULL; cluster_id_name <- NULL
804+
if (!is.null(id_ind) && !is.null(id_se)) {
805+
split_id_vec <- id_ind$vec; split_id_name <- id_ind$name
806+
cluster_id_vec <- id_se$vec; cluster_id_name <- id_se$name
807+
} else if (!is.null(id_ind)) {
808+
split_id_vec <- id_ind$vec; split_id_name <- id_ind$name
809+
cluster_id_vec <- id_ind$vec; cluster_id_name <- id_ind$name
810+
} else if (!is.null(id_se)) {
811+
split_id_vec <- id_se$vec; split_id_name <- id_se$name
812+
cluster_id_vec <- id_se$vec; cluster_id_name <- id_se$name
813+
}
814+
815+
if (!is.null(split_id_vec)) {
816+
if (anyNA(split_id_vec)) {
817+
stop("The fold-splitting identifier ('", split_id_name, "') contains NA ",
818+
"values. All observations must have an identifier.")
784819
}
785-
786-
if (anyNA(individual_id_vec)) {
787-
stop("individual_id contains NA values. All observations must have an individual identifier.")
820+
if (anyNA(cluster_id_vec)) {
821+
stop("The SE-clustering identifier ('", cluster_id_name, "') contains NA ",
822+
"values. All observations must have an identifier.")
788823
}
789-
790-
n_individuals <- length(unique(individual_id_vec))
791-
if (n_individuals == n) {
792-
warning("Every observation has a unique individual_id (", n_individuals,
793-
" unique IDs for ", n, " observations). ",
794-
"This is equivalent to no panel structure. ",
795-
"If your data is not panel data, you can omit individual_id.")
824+
n_split_groups <- length(unique(split_id_vec))
825+
if (n_split_groups == n && !(!is.null(id_ind) && !is.null(id_se))) {
826+
warning("Every observation has a unique '", split_id_name, "' value (",
827+
n_split_groups, " for ", n, " observations). ",
828+
"This is equivalent to no panel structure for fold splitting.")
796829
}
797-
if (n_individuals < K) {
798-
stop("Number of unique individuals (", n_individuals,
830+
if (n_split_groups < K) {
831+
stop("Number of unique '", split_id_name, "' groups (", n_split_groups,
799832
") must be at least K (", K, ")")
800833
}
801-
834+
message("ensembleHTE: splitting cross-fitting folds by '", split_id_name,
835+
"' (", n_split_groups, " groups); clustering standard errors by '",
836+
cluster_id_name, "' (", length(unique(cluster_id_vec)), " clusters).")
802837
}
803838

804839
# Propensity score handling
@@ -874,7 +909,7 @@ ensemble_hte <- function(formula = NULL, treatment = NULL, data = NULL,
874909
# Split the sample - stratify by D and train_idx
875910
stratify_var <- interaction(D, train_idx, drop = TRUE)
876911
splits <- create_folds(n, M = M, K = K, stratify_var = stratify_var,
877-
cluster_id = individual_id_vec)
912+
cluster_id = split_id_vec)
878913

879914
# Train learners and compute ensemble predictions
880915
ite_cols <- paste0("ite_", algorithms)
@@ -918,7 +953,10 @@ ensemble_hte <- function(formula = NULL, treatment = NULL, data = NULL,
918953
ite_rep <- rowMeans(predictions_m[, ite_cols, drop = FALSE], na.rm = TRUE)
919954
} else if (ensemble_strategy == "cv") {
920955
# Cross-validated BLP ensemble
921-
ens_splits <- create_folds(n, M = 1, K = ensemble_folds, stratify_var = splits[[m]])[[1]]
956+
# Cluster by individual_id_vec (when panel) so all observations for the same
957+
# unit stay in the same ensemble fold, mirroring the outer cross-fitting split.
958+
ens_splits <- create_folds(n, M = 1, K = ensemble_folds, stratify_var = splits[[m]],
959+
cluster_id = split_id_vec)[[1]]
922960
dt_ens <- data.table(
923961
Y = Y,
924962
D = D,
@@ -1085,7 +1123,10 @@ ensemble_hte <- function(formula = NULL, treatment = NULL, data = NULL,
10851123
prop_score = prop_score,
10861124
weights = W,
10871125
train_idx = train_idx,
1088-
individual_id = individual_id_vec,
1126+
individual_id = split_id_vec,
1127+
individual_id_name = split_id_name,
1128+
se_cluster_id = cluster_id_vec,
1129+
se_cluster_id_name = cluster_id_name,
10891130
splits = splits,
10901131
n = n,
10911132
n_train = n_train,
@@ -1207,8 +1248,9 @@ print.ensemble_hte_fit <- function(x, ...) {
12071248
cat(" Treatment: ", x$treatment, "\n", sep = "")
12081249
cat(" Covariates: ", ncol(x$X), "\n", sep = "")
12091250
if (!is.null(x$individual_id)) {
1210-
n_individuals <- length(unique(x$individual_id))
1211-
cat(" Panel data: ", n_individuals, " individuals\n", sep = "")
1251+
split_lbl <- if (!is.null(x$individual_id_name)) x$individual_id_name else "individual"
1252+
cat(" Split by: ", split_lbl, " (",
1253+
length(unique(x$individual_id)), " groups)\n", sep = "")
12121254
}
12131255
cat("\n")
12141256
cat("Model specification:\n")
@@ -1229,8 +1271,13 @@ print.ensemble_hte_fit <- function(x, ...) {
12291271
}
12301272
cat(" Covariate scaling: ", scale_status, "\n", sep = "")
12311273
cat(" Hyperparameter tuning: ", tune_status, "\n", sep = "")
1232-
if (!is.null(x$individual_id)) {
1233-
cat(" Standard errors: cluster-robust (at individual level)\n", sep = "")
1274+
cluster_vec <- if (!is.null(x$se_cluster_id)) x$se_cluster_id else x$individual_id
1275+
if (!is.null(cluster_vec)) {
1276+
cl_lbl <- if (!is.null(x$se_cluster_id_name)) x$se_cluster_id_name
1277+
else if (!is.null(x$individual_id_name)) x$individual_id_name
1278+
else "individual"
1279+
cat(" Std. errors: cluster-robust by ", cl_lbl, " (",
1280+
length(unique(cluster_vec)), " clusters)\n", sep = "")
12341281
}
12351282
if (!is.null(x$store_baseline) && x$store_baseline != "none") {
12361283
baseline_desc <- switch(x$store_baseline,

0 commit comments

Comments
 (0)