Skip to contents

Extends bd-HSIC and DR-DATE/DR-DETT to hierarchical (nested/clustered) data by decomposing the test statistic into within-cluster and between-cluster components.

Usage

hierarchical_test(
  y,
  treatment,
  covariates,
  cluster_id,
  method = c("dr-date", "dr-dett", "bd-hsic"),
  kernel_y = kernel_spec(),
  n_permutations = 500L,
  weight_method = c("equal", "icc", "within_only"),
  seed = NULL,
  verbose = FALSE,
  ...
)

Arguments

y

Numeric vector or matrix. Outcome.

treatment

Treatment variable (binary for DR tests, any for bd-HSIC).

covariates

Numeric matrix of confounders.

cluster_id

Factor or integer vector identifying clusters.

method

Character. "dr-date" (default), "dr-dett", or "bd-hsic".

kernel_y

Kernel specification for outcomes.

n_permutations

Integer. Number of permutations. Default is 500.

weight_method

Character. How to weight within/between components: "equal" (default), "icc" (variance decomposition), or "within_only".

seed

Integer or NULL.

verbose

Logical.

...

Additional arguments passed to the underlying test. Do not pass cross_fit or n_folds: within-cluster sub-tests are always fit in-sample, as the top-level within-cluster permutation supplies the calibration.

Value

An object of class "kernel_test_result" with additional hierarchical component containing within/between statistics.

Details

For clustered data (e.g., patients within hospitals, plots within farms), standard kernel tests may have inflated type I error because observations within the same cluster are not independent.

This function decomposes the test into:

  • Within-cluster: Average of within-cluster test statistics (tests for treatment effects within each cluster).

  • Between-cluster: Test on cluster-level mean embeddings (tests for treatment effects across clusters).

The combined statistic is a weighted sum, with weights determined by weight_method. Permutation is performed within clusters to preserve the hierarchical structure.

See also

Examples

# \donttest{
set.seed(42)
n_clusters <- 20
n_per <- 30
n <- n_clusters * n_per
cluster_id <- rep(1:n_clusters, each = n_per)

# Cluster-level random effects
cluster_effect <- rnorm(n_clusters, sd = 1)[cluster_id]
x <- matrix(rnorm(n * 2), n, 2)
t <- rbinom(n, 1, plogis(0.3 * x[, 1]))
y <- 0.5 * t + cluster_effect + x[, 1] + rnorm(n)

result <- hierarchical_test(y, t, x, cluster_id,
  method = "dr-date",
  n_permutations = 100,
  seed = 1
)
print(result)
#> 
#>    Hierarchical-DRDATE Test
#> 
#> Statistic: 0.0309451 
#> P-value:   0.0099 
#> N:         600 
#> Perms:     100 
#> Kernel Y:  rbf (bw = 1.822)
#> 
# }