A Practical Guide to sreg

0.1 Overview

The sreg package estimates average treatment effects (ATEs) in stratified randomized experiments under covariate-adaptive randomization. It supports:

• individual- and cluster-level treatment assignment;
• one or multiple active treatment arms relative to control;
• experiments without strata and experiments with large strata;
• uniform small-strata designs, including matched pairs and general k-tuples;
• mixed designs containing both small and large strata;
• estimation with or without linear covariate adjustment;
• heteroskedasticity-robust standard errors with an optional HC1 correction;
• simulation through sreg.rgen(); and
• formatted summaries and confidence-interval plots.

The package’s main functions are sreg(), which performs estimation and inference, and sreg.rgen(), which generates simulated experiments. The print() and plot() methods summarize fitted sreg objects.

library(sreg)
#> Package 'sreg' version 2.1.0
#> Type 'citation("sreg")' for citing this R package in publications.

0.2 The basic workflow

A typical analysis has four steps:

1. Identify the randomization unit and the stratification design.
2. Prepare the outcome, treatment, stratum, cluster, and covariate variables.
3. Call sreg() with arguments matching the experimental design.
4. Inspect, print, or plot the returned object.

The core call is:

sreg(
  Y, S = NULL, D,
  G.id = NULL, Ng = NULL,
  X = NULL, HC1 = TRUE,
  small.strata = FALSE, k = NULL
)

0.3 Large-strata designs

When small.strata = FALSE, sreg() applies the large-strata estimator. This is the default. The following example generates an individual-level experiment with four strata, a control arm, and two active treatments.

set.seed(101)

large_data <- sreg.rgen(
  n = 600,
  tau.vec = c(0.3, 0.7),
  n.strata = 4,
  cluster = FALSE
)

head(large_data)
#>          Y S D      x_1      x_2
#> 1 2.154491 1 0 3.527956 2.336471
#> 2 4.238873 3 2 7.534395 2.274160
#> 3 2.552000 2 2 5.528556 1.166914
#> 4 3.726746 2 2 4.390154 2.802584
#> 5 3.882640 2 2 5.534880 2.172677
#> 6 3.392253 3 0 2.128014 1.625707
table(large_data$S, large_data$D)
#>    
#>      0  1  2
#>   1 30 28 28
#>   2 67 67 67
#>   3 69 69 69
#>   4 36 35 35

fit_large_unadjusted <- sreg(
  Y = large_data$Y,
  S = large_data$S,
  D = large_data$D
)

fit_large_unadjusted$tau.hat
#> [1] 0.4844155 0.7490444
fit_large_unadjusted$se.rob
#> [1] 0.1441396 0.1480421

X_large <- large_data[c("x_1", "x_2")]

fit_large <- sreg(
  Y = large_data$Y,
  S = large_data$S,
  D = large_data$D,
  X = X_large
)

print(fit_large)
#> Saturated Model Estimation Results under CAR with linear adjustments
#> Observations: 600 
#> Number of treatments: 2 
#> Number of strata: 4 
#> Setup: large strata 
#> Standard errors: adjusted (HC1) 
#> Treatment assignment: individual level 
#> Covariates used in linear adjustments: x_1, x_2
#> ---
#> Coefficients:
#>       Tau   As.se  T-stat P-value CI.left(95%) CI.right(95%) Significance
#> 1 0.47972 0.10108 4.74581       0      0.28160       0.67783          ***
#> 2 0.84412 0.10318 8.18100       0      0.64189       1.04635          ***
#> ---
#> Signif. codes:  0 `***` 0.001 `**` 0.01 `*` 0.05 `.` 0.1 ` ` 1

0.4 Plotting treatment effects

Calling plot() on an sreg object displays the estimated ATEs and their confidence intervals. Treatment labels, colors, axes, grid lines, and other graphical elements can be customized.

plot(
  fit_large,
  treatment_labels = c("Program A", "Program B"),
  title = "Estimated treatment effects",
  x_axis_title = "ATE relative to control",
  bar_fill = c("#3B82F6", "#14B8A6"),
  point_fill = "white",
  zero_line = TRUE
)

The method invisibly returns a ggplot object, so it can be stored or further modified:

p <- plot(fit_large, treatment_labels = c("Program A", "Program B"))
p + ggplot2::theme_minimal()

See help("plot.sreg") for all graphical arguments.

0.5 Uniform small-strata designs

Set small.strata = TRUE for designs such as matched pairs and k-tuples. In a uniform small-strata design, all strata have the same size. The package observes that common size directly, so the k argument is optional.

set.seed(102)

small_data <- sreg.rgen(
  n = 300,
  tau.vec = c(0.3, 0.7),
  cluster = FALSE,
  small.strata = TRUE,
  k = 3,
  treat.sizes = c(1, 1, 1)
)

table(table(small_data$S))
#> 
#>   3 
#> 100

Here, every stratum is a triplet containing one observation in each of the three treatment arms. Estimation does not require k = 3 because the common size is unambiguous:

fit_small <- sreg(
  Y = small_data$Y,
  S = small_data$S,
  D = small_data$D,
  X = small_data[c("x_1", "x_2")],
  small.strata = TRUE
)

print(fit_small)
#> Saturated Model Estimation Results under CAR with linear adjustments
#> Observations: 300 
#> Number of treatments: 2 
#> Number of strata: 100 
#> Setup: small strata 
#> Strata size (k): 3 
#> Standard errors: adjusted (HC1) 
#> Treatment assignment: individual level 
#> Covariates used in linear adjustments: x_1, x_2
#> ---
#> Coefficients:
#>       Tau   As.se  T-stat P-value CI.left(95%) CI.right(95%) Significance
#> 1 0.19561 0.14924 1.31072 0.18995     -0.09689       0.48810             
#> 2 0.85428 0.13635 6.26530 0.00000      0.58704       1.12153          ***
#> ---
#> Signif. codes:  0 `***` 0.001 `**` 0.01 `*` 0.05 `.` 0.1 ` ` 1

Supplying k = 3 gives the same design choice and adds a validation check. If the supplied value does not match the observed common size, sreg() returns an error.

sreg(
  Y = small_data$Y,
  S = small_data$S,
  D = small_data$D,
  small.strata = TRUE,
  k = 3
)

The small-strata estimators are intended for experiments with a fixed treatment allocation within each stratum. The researcher remains responsible for verifying that the supplied assignments follow the randomization protocol; sreg() does not independently verify every feature of that protocol.

0.6 Mixed small- and large-strata designs

A mixed design contains a recurring small-stratum component and a large-strata component. sreg.rgen() can generate both components in one call. For individual-level designs, n.small is the number of observations assigned to small strata. For cluster-randomized designs, it is the number of clusters assigned to small strata.

The next example uses 80 observations in 20 four-tuples and distributes the remaining 40 observations across four large strata.

set.seed(103)

mixed_data <- sreg.rgen(
  n = 120,
  tau.vec = 0.5,
  cluster = FALSE,
  mixed.strata = TRUE,
  n.small = 80,
  k = 4,
  treat.sizes = c(2, 2),
  n.strata = 4
)

sort(table(mixed_data$S))
#> 
#>  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 24 22 23 
#>  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  6  6 10 18

Because this mixed design uses 4-tuples, pass the known value k = 4 to sreg():

fit_mixed <- suppressWarnings(sreg(
  Y = mixed_data$Y,
  S = mixed_data$S,
  D = mixed_data$D,
  small.strata = TRUE,
  k = 4
))

print(fit_mixed)
#> Saturated Model Estimation Results under CAR
#> Observations: 120 
#> Number of treatments: 1 
#> Number of strata: 24 
#> Setup: mixed design (includes both small and large strata) 
#> Strata size (k, small strata only): 4 
#> Standard errors: adjusted (HC1) 
#> Treatment assignment: individual level 
#> Covariates used in linear adjustments: 
#> ---
#> Coefficients:
#>       Tau   As.se  T-stat P-value CI.left(95%) CI.right(95%) Significance
#> 1 0.43248 0.27423 1.57707 0.11478       -0.105       0.96996             
#> ---
#> Signif. codes:  0 `***` 0.001 `**` 0.01 `*` 0.05 `.` 0.1 ` ` 1

Normally, this call emits an informative warning stating that a mixed design was detected, identifying k = 4 as the small-stratum size, and noting that weighted estimators are used. The warning is suppressed above only to keep the vignette output compact.

0.7 Cluster-randomized experiments

Supply G.id when treatment is assigned at the cluster level. Each cluster must belong to one stratum and receive one treatment. Ng contains cluster sizes. When Ng = NULL, the package uses the number of observed records in each cluster and reports a warning.

set.seed(104)

cluster_data <- sreg.rgen(
  n = 60,
  tau.vec = 0.5,
  n.strata = 4,
  cluster = TRUE
)

c(
  observations = nrow(cluster_data),
  clusters = length(unique(cluster_data$G.id)),
  strata = length(unique(cluster_data$S))
)
#> observations     clusters       strata 
#>         1780           60            4

fit_cluster <- sreg(
  Y = cluster_data$Y,
  S = cluster_data$S,
  D = cluster_data$D,
  G.id = cluster_data$G.id,
  Ng = cluster_data$Ng,
  X = cluster_data[c("x_1", "x_2")]
)

print(fit_cluster)
#> Saturated Model Estimation Results under CAR with linear adjustments
#> Observations: 1780 
#> Clusters: 60 
#> Number of treatments: 1 
#> Number of strata: 4 
#> Setup: large strata 
#> Standard errors: adjusted (HC1) 
#> Treatment assignment: cluster level 
#> Covariates used in linear adjustments: x_1, x_2
#> ---
#> Coefficients:
#>      Tau   As.se  T-stat P-value CI.left(95%) CI.right(95%) Significance
#> 1 0.2872 0.18441 1.55741 0.11937     -0.07423       0.64863             
#> ---
#> Signif. codes:  0 `***` 0.001 `**` 0.01 `*` 0.05 `.` 0.1 ` ` 1
#> [1] TRUE

For cluster-randomized small or mixed designs, stratum size k means the number of distinct clusters per small stratum—not the number of individual observations. Individual-level covariates cannot be used directly for cluster-level adjustment; if supplied covariates vary within cluster, sreg() aggregates them to cluster averages and reports this behavior.

0.8 Experiments without stratification

Set S = NULL when treatment was randomized without strata. This option uses the large-strata pathway with the entire experiment treated as one population:

fit_unstratified <- sreg(
  Y = Y,
  S = NULL,
  D = D,
  X = X,
  small.strata = FALSE
)

Do not combine S = NULL with small.strata = TRUE, because small-strata inference requires stratum identifiers.

0.9 The sreg result object

In addition to its printed summary, an sreg object exposes components useful for reporting and further analysis:

Mixed-design results also contain res.small, res.big, and mixed.design, which provide the component fits and identify the combined estimator.

names(fit_mixed)
#>  [1] "tau.hat"      "se.rob"       "t.stat"       "p.value"      "as.CI"       
#>  [6] "beta.hat"     "CI.left"      "CI.right"     "ols.iter"     "lin.adj"     
#> [11] "data"         "res.small"    "res.big"      "mixed.design" "small.strata"
#> [16] "HC1"
fit_mixed$tau.hat
#> [1] 0.4324797
fit_mixed$se.rob
#> [1] 0.2742298
fit_mixed$mixed.design
#> [1] TRUE

0.10 HC1 correction

HC1 = TRUE, the default, applies the package’s small-sample correction to the variance estimator. Set HC1 = FALSE to obtain the corresponding uncorrected variance estimate:

fit_hc1 <- sreg(Y = Y, S = S, D = D, HC1 = TRUE)
fit_no_hc1 <- sreg(Y = Y, S = S, D = D, HC1 = FALSE)

The choice should follow the inferential specification for the analysis rather than whichever option yields a preferred standard error.

0.11 Working with the built-in empirical data

The package includes AEJapp, replication data associated with Chong et al. (2016). The following prepares the outcome, treatment, and stratum variables used in the package example. The original treatment value 3 is recoded as control (0) to match the required treatment indexing.

data("AEJapp", package = "sreg")

Y_emp <- as.numeric(unclass(AEJapp$gradesq34))
D_raw <- as.numeric(unclass(AEJapp$treatment))
D_emp <- ifelse(D_raw == 3, 0, D_raw)
S_emp <- as.numeric(unclass(AEJapp$class_level))

fit_empirical <- sreg(
  Y = Y_emp,
  S = S_emp,
  D = D_emp
)

fit_empirical$tau.hat
#> [1] -0.05112971  0.40903373
fit_empirical$se.rob
#> [1] 0.2064541 0.2065146

See help("AEJapp") for the data description and source.

0.12 Simulation with sreg.rgen()

sreg.rgen() is intended for examples, Monte Carlo exercises, and learning how different designs are represented. Important arguments include:

For small strata, length(treat.sizes) must equal length(tau.vec) + 1, and its entries must sum to k. In a mixed design, n.small must be divisible by k. Use set.seed() when reproducibility matters.

0.13 Common warnings and errors

• Mixed design detected: the function has split the data into small and large components and is using the combined estimator.
• Too few strata qualify as small: fewer than 25% of strata have the selected small size, or k was omitted for a mixed design with k-tuples larger than three.
• small.strata = FALSE warning: the observed design appears to contain a substantial small-strata component, but the large-strata estimator was requested.
• Missing S: small-strata inference cannot be performed without stratum indicators.
• Missing Ng: observed cluster counts are substituted for cluster sizes.
• Covariates do not vary sufficiently: adjusted estimation may be infeasible, in which case the package may proceed without adjustment or ask the user to reduce the covariate set.
• Too many covariates: some stratum-treatment cells may not support the requested regressions; reduce X or fit the unadjusted estimator.

Warnings about the design should not be suppressed in substantive analyses until the researcher has confirmed that the arguments match the actual randomization procedure.

0.14 Recommended analysis checklist

Before reporting results:

1. Confirm that D uses 0 for control and consecutive positive integers for active treatments.
2. Confirm that S uses consecutive positive integers beginning at 1.
3. Verify whether assignment occurred at the individual or cluster level.
4. For cluster assignment, verify G.id, Ng, and that treatment is constant within cluster.
5. Decide whether the design has large, uniform small, or mixed strata.
6. For mixed general k-tuples, supply the design value through k.
7. Verify the fixed treatment allocation used in small strata.
8. Check whether covariates vary sufficiently for adjustment and whether they are measured at the appropriate level.
9. Review all warnings before interpreting standard errors.
10. Record the package version and use a reproducible seed for simulations.

0.15 References and further help

Use the following commands for function-level documentation:

help("sreg")
help("sreg.rgen")
help("plot.sreg")
help("AEJapp")
citation("sreg")

The package help files and supplementary materials contain the formal references and estimator expressions underlying the supported designs.