Estimates an average treatment effect as a difference in means, choosing
the point estimator, the variance, and the degrees of freedom that match
the randomization rather than requiring you to name them. Unit, cluster,
blocked, block-cluster, matched-pair, and matched-pair cluster designs are
recognised, and the design element of the result reports which case
applied.
Usage
difference_in_means(
formula,
data,
blocks,
clusters,
weights,
subset,
se_type = c("default", "none"),
condition1 = NULL,
condition2 = NULL,
ci = TRUE,
alpha = 0.05
)Arguments
- formula
(required) An object of class formula with one variable on the RHS
- data
(optional) A
data.frame- blocks
(optional) A bare (unquoted) name of the block variable
- clusters
(optional) A bare (unquoted) name of the cluster variable
- weights
(optional) The bare (unquoted) name of the weights variable
- subset
(optional) A bare (unquoted) expression specifying a subset
- se_type
(optional)
"default"or"none"- condition1
(optional) Value in treatment for the control condition
- condition2
(optional) Value in treatment for the treatment condition
- ci
(optional) Logical. Whether to compute p-values and confidence intervals.
- alpha
(optional) The significance level, 0.05 by default.
Value
An object of class "difference_in_means", a list holding
coefficients, std.error, df, statistic, p.value, conf.low,
conf.high, term, outcome, condition1, condition2, vcov,
nobs, alpha, and design, a string naming the case that applied:
"Standard", "Blocked", "Small blocks", "Hybrid blocked",
"Matched-pair", "Clustered", "Block-clustered", or
"Matched-pair clustered".
Details
Selects the appropriate point estimate, standard errors, and degrees of freedom for unit randomized, cluster randomized, block randomized, block-cluster randomized, matched-pairs, and matched-pair cluster randomized designs.
Blocks of different sizes. For unit randomized blocks, blocks are
classified by how many units each arm holds rather than by how large the
block is. A block with at least two treated and two control units has its
own Neyman variance. A block with a single treated or single control
unit has no estimable within-block variance, so the variation across such
blocks stands in for it. A design containing both kinds combines the two
parts by squared share of the sample, following Pashley and Miratrix
(2021). The design element of the returned object reports which case
applied: "Blocked", "Matched-pair", "Small blocks", or
"Hybrid blocked".
Two designs are refused, because the variance genuinely cannot be
estimated: exactly one block with a singleton arm, and a set of
different-sized such blocks in which one holds half or more of their
units. Both messages suggest merging blocks or using lm_robust() with
block fixed effects.
If weights are specified, estimation is handed to lm_robust() with HC2
standard errors.
Blocks of clusters. Pashley and Miratrix treat treatment assigned
within blocks, not blocks of clusters, so blocked designs with clusters
use the earlier estimators. Every block must have at least two treated and
two control clusters, unless the design is matched-pair clustered, where
the variance is estimated across blocks. A block with a single treated or
control cluster is refused: its within-block variance is not estimable,
and estimating it anyway understates the standard error by roughly the
block's cluster count.
References
Gerber, Alan S. and Donald P. Green. 2012. Field Experiments: Design, Analysis, and Interpretation. New York: W.W. Norton.
Imai, Kosuke, Gary King, and Clayton Nall. 2009. "The Essential Role of Pair Matching in Cluster-Randomized Experiments." Statistical Science 24(1): 29-53. doi:10.1214/08-STS274 .
Pashley, Nicole E. and Luke W. Miratrix. 2021. "Insights on Variance Estimation for Blocked and Matched Pairs Designs." Journal of Educational and Behavioral Statistics 46(3): 271-296. doi:10.3102/1076998620946272 .
Examples
set.seed(30)
dat <- data.frame(y = rnorm(100), z = rep(0:1, 50))
# Unblocked, unclustered: the Welch-corrected two-sample difference
fit <- difference_in_means(y ~ z, data = dat)
fit
#> Design: Standard
#> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
#> z -0.1864076 0.2119164 -0.8796279 0.3812186 -0.6069631 0.2341479 97.73532
fit$design
#> [1] "Standard"
# Blocked designs use the Neyman variance within each block
dat_bl <- data.frame(
bl = rep(1:10, each = 10),
z = rep(rep(0:1, each = 5), times = 10)
)
dat_bl$y <- rnorm(100) + 0.3 * dat_bl$z
difference_in_means(y ~ z, data = dat_bl, blocks = bl)
#> Design: Blocked
#> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
#> z 0.03347738 0.2101101 0.1593326 0.8738084 -0.3846551 0.4516098 80
# Matched pairs are recognised as such
dat_pr <- data.frame(pr = rep(1:50, each = 2), z = rep(c(0, 1), 50))
dat_pr$y <- rnorm(100) + 0.3 * dat_pr$z
difference_in_means(y ~ z, data = dat_pr, blocks = pr)$design
#> [1] "Matched-pair"
# Blocks of unequal shape, which earlier versions refused, use the
# Pashley and Miratrix (2021) estimators. `design` reports which case
# applied rather than leaving it to be inferred from the block sizes.
dat_hy <- rbind(dat_bl[c("bl", "z", "y")],
transform(dat_pr[c("pr", "z", "y")], bl = pr + 100)[c("bl", "z", "y")])
difference_in_means(y ~ z, data = dat_hy, blocks = bl)$design
#> [1] "Hybrid blocked"
# Clustered assignment
dat_cl <- data.frame(cl = rep(1:20, each = 5))
dat_cl$z <- rep(rep(0:1, each = 5), times = 10)
dat_cl$y <- rnorm(100) + 0.3 * dat_cl$z
difference_in_means(y ~ z, data = dat_cl, clusters = cl)
#> Design: Clustered
#> Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
#> z 0.6286243 0.2326784 2.701688 0.01459924 0.1397852 1.117463 18