# Model-robust and efficient covariate adjustment for cluster-randomized
experiments

Bingkai Wang1 Chan Park1 Dylan S. Small1 and Fan Li2 1The Statistics and Data
Science Department of the Wharton School  University of Pennsylvania
Philadelphia  PA  USA 2Department of Biostatistics and Center for Methods in
Implementation and Prevention Science  Yale University School of Public Health
New Haven  CT  USA

###### Abstract

Cluster-randomized experiments are increasingly used to evaluate interventions
in routine practice conditions, and researchers often adopt model-based
methods with covariate adjustment in the statistical analyses. However, the
validity of model-based covariate adjustment is unclear when the working
models are misspecified, leading to ambiguity of estimands and risk of bias.
In this article, we first adapt two conventional model-based methods,
generalized estimating equations and linear mixed models, with weighted
g-computation to achieve robust inference for cluster-average and individual-
average treatment effects. Furthermore, we propose an efficient estimator for
each estimand that allows for flexible covariate adjustment and additionally
addresses cluster size variation dependent on treatment assignment and other
cluster characteristics. Such cluster size variations often occur post-
randomization and, if ignored, can lead to bias of model-based estimators. For
our proposed estimator, we prove that when the nuisance functions are
consistently estimated by machine learning algorithms, the estimator is
consistent, asymptotically normal, and efficient. When the nuisance functions
are estimated via parametric working models, the estimator is triply-robust.
Simulation studies and analyses of three real-world cluster-randomized
experiments demonstrate that the proposed methods are superior to existing
alternatives.

Keywords: Estimand; Efficiency; Group-randomized trial; Covariate adjustment;
Unequal cluster size.

##  1 Introduction

Cluster-randomized experiments refer to study designs that randomize treatment
at the cluster level; clusters can be villages, hospitals, or worksites
(Murray et al., , 1998; Donner and Klar, , 2000). Cluster randomization is
often used to study group-level interventions or to prevent treatment
contamination, and is increasingly adopted in pragmatic clinical trials
evaluating interventions in routine practice conditions. For example, among
$22$ demonstration projects across different disease areas supported by the
United States National Institute of Health Pragmatic Clinical Trials
Collaboratory since 2012, $16$ studies adopted cluster randomization.

In the analyses of cluster-randomized experiments, covariate adjustment is
essential to address baseline chance imbalances and improve the precision of
treatment effect estimators. However, challenges in covariate adjustment
persist due to the complexity associated with the multilevel data structure
under cluster randomization. First, although model-based methods, including
generalized estimating equations (Liang and Zeger, , 1986) and generalized
linear mixed models (Breslow and Clayton, , 1993), are commonly used to
perform covariate adjustment in cluster-randomized experiments, it remains
unclear whether the treatment effect coefficient corresponds to a clearly-
defined estimand of interest, especially when the working model is
misspecified. Even in the absence of covariates, Wang et al.,  (2022) has
demonstrated that the treatment coefficient from generalized estimating
equations with an exchangeable working correlation structure corresponds to an
ambiguous estimand when the cluster size is informative for causal effects.
Second, participants are frequently sampled after cluster randomization such
that the observed cluster size may depend on treatment assignment and other
cluster attributes. Failure to address such cluster-dependent sampling schemes
can lead to bias of treatment effect estimators from standard techniques. For
example, when the observed cluster size depends on cluster-level covariates,
Bugni et al.,  (2022) demonstrated that the standard difference-in-means
estimator is biased. To date, robust causal inference methods that can
simultaneously maximize the precision gain from covariate adjustment and
address cluster-dependent sampling are unavailable for cluster-randomized
experiments.

In this article, we address the aforementioned challenges in estimating two
classes of causal estimands: the cluster-average treatment effect and
individual-average treatment effect. The former estimand addresses the
question of “what is the expected change in outcome associated with treatment
for a typical cluster with its natural source population?” and gives equal
weight to each cluster, whereas the latter estimand addresses the question of
“what is the expected change in outcome associated with treatment for a
typical individual?” and gives equal weight to each individual. The two
estimands represent treatment effects at different levels and differ when
there is treatment effect heterogeneity by cluster size (Kahan et al., ,
2022). To estimate these estimands, we first adapt generalized estimating
equations and linear mixed models through weighted g-computation, and provide
a set of sufficient conditions to achieve model-robust inference, i.e.,
consistency and asymptotic normality under arbitrary working model
misspecification. These new insights help clarify when conventional multilevel
regression models, which are routinely used as standard practice (Turner et
al., , 2017), provide valid average causal effect estimators even if the
regression model formulation differs from the unknown and potentially complex
data generating process.

The weighted g-computation estimators, however, can be inefficient for
covariate adjustment and become biased when the sufficient conditions for
model robustness do not hold, e.g., when the observed cluster size depends on
treatment and covariates. To address such limitations, we characterize the
efficient influence function for both estimands and propose efficient
estimators that allow for flexible covariate adjustment and address cluster-
dependent sampling. When the nuisance functions are estimated by parametric
models, such as generalized estimating equations or generalized linear mixed
models, the proposed estimators are triply-robust, that is, they are
consistent to the specified causal estimand if two out of the three nuisance
functions are consistently estimated. When the nuisance working models are all
consistently estimated, for example, by machine learning algorithms with
cross-fitting, the proposed estimators achieve the efficiency lower bound
given our causal models.

Our results build on but differ from the existing literature on causal
inference for cluster-randomized experiments. Imai et al.,  (2009) and
Middleton and Aronow,  (2015) proposed cluster-level methods for estimating
the average treatment effects but did not consider covariate adjustment to
improve efficiency. Schochet et al.,  (2021) established the finite-population
central limit theorem for linearly-adjusted estimators, and Su and Ding,
(2021) extended their results by providing a unified theory for a class of
weighted average treatment effect estimands. These prior works considered a
finite-population framework with linear working models and an independence
working correlation structure. In contrast, we study causal effect estimators
under a super-population framework and address robust and efficient estimation
with a wider class of working models. Balzer et al.,  (2021) and Benitez et
al.,  (2021) applied targeted maximum likelihood estimation under hierarchical
structural models for cluster-randomized experiments. However, they did not
consider cluster size variation arising from cluster-dependent sampling that
occurs post-randomization. Bugni et al.,  (2022) developed moment-based
estimators addressing cluster-dependent sampling but did not consider
covariate adjustment. We therefore improve these existing contributions by
proposing efficient estimators to achieve flexible covariate adjustment and
accommodate post-randomization cluster-dependent sampling.

The remainder of the article is organized as follows. In Section 2, we
describe our super-population framework, present the causal estimands of
interest, and structural assumptions for identification. Section 3 adapts
generalized estimating equations and linear mixed models to target our
estimands. In Section 4, we develop the efficient influence functions and
propose our efficient estimators. In Sections 5 and 6, we demonstrate our
theoretical results via simulation studies and re-analyses of three cluster-
randomized trials. Section 7 concludes with a discussion.

##  2 Notation, estimands, and assumptions

We consider a cluster-randomized experiment with $m$ clusters. For each
cluster $i$, we let $N_{i}$ denote the total number of individuals in the
underlying source population, $M_{i}$ be the observed number of individuals
sampled into the study, $A_{i}\in\\{0,1\\}$ be the cluster-level treatment
indicator, and $C_{i}$ be a $q$-dimensional vector of cluster-level
covariates. For each individual $j$ in cluster $i$, we define $Y_{ij}$ as the
outcome, $X_{ij}$ as a $p$-dimensional vector of individual-level covariates,
and $S_{ij}\in\\{0,1\\}$ as the sampling indicator, i.e., recruited into the
experiment. By definition, the observed cluster size
$M_{i}=\sum_{j=1}^{N_{i}}S_{ij}\leq N_{i}$.

We proceed under the potential outcome framework (Neyman et al., , 1990) and
define $Y_{ij}(a)$ as the potential outcome and $S_{ij}(a)$ as the potential
sampling indicator if cluster $i$ were assigned to treatment group
$a\in\\{0,1\\}$. We assume consistency such that
$Y_{ij}=A_{i}Y_{ij}(1)+(1-A_{i})Y_{ij}(0)$ and
$S_{ij}=A_{i}S_{ij}(1)+(1-A_{i})S_{ij}(0)$. As a result, the potential
observed cluster size in a treated cluster, denoted as
$M_{i}(1)=\sum_{j=1}^{N_{i}}S_{ij}(1)$, can be different from its counterpart
in a control cluster, denoted as $M_{i}(0)=\sum_{j=1}^{N_{i}}S_{ij}(0)$.
Defining $Y_{i}(a)=\\{Y_{i1}(a),\ldots,Y_{iN_{i}}(a)\\}\in\mathbb{R}^{N_{i}}$
as the collection of potential outcomes,
$S_{i}(a)=\\{S_{i1}(a),\ldots,S_{iN_{i}}(a)\\}\in\mathbb{R}^{N_{i}}$ as the
collection of potential sampling indicators, and
$\mathrm{X}_{i}=(X_{i1},\ldots,X_{iN_{i}})^{\top}\in\mathbb{R}^{N_{i}\times
p}$, we write the collection of random variables in a cluster as
$W_{i}=\\{Y_{i}(1),Y_{i}(0),S_{i}(1),S_{i}(0),N_{i},C_{i},\mathrm{X}_{i}\\}$.
We introduce the following assumptions on the complete, while not fully
observed, data $\\{(W_{1},A_{1}),\dots,(W_{m},A_{m})\\}$.

######  Assumption 1 (Super-population).

(a) Random variables $W_{1},\dots,W_{m}$ are mutually independent. (b) The
source population size $N_{i}$ follows an unknown distribution
$\mathcal{P}^{N}$ over a finite support on $\mathbb{N}^{+}$. (c) Given
$N_{i}$, $W_{i}$ follows an unknown distribution $\mathcal{P}^{W\mid N}$ with
finite second moments.

######  Assumption 2 (Cluster randomization).

The treatment indicator $A_{i}$ for each cluster is an independent realization
from a Bernoulli distribution $\mathcal{P}^{A}$ with $pr(A=1)=\pi\in(0,1)$.
Furthermore, $(A_{1},\dots,A_{m})$ is independent of $(W_{1},\dots,W_{m})$.

######  Assumption 3 (Cluster-dependent sampling).

For $a\in\\{0,1\\}$, the potential observed cluster size
$M(a)=h_{a}(N,C,\epsilon_{a})$ for some unknown function $h_{a}$ and exogenous
random noise $\epsilon_{a}$ that is independent of
$\\{N,C,\mathrm{X},Y(a)\\}$. Further, for each possible $N$-dimensional binary
vector ${s}$ with $M(a)$ ones,

$pr\\{S(a)={s}\mid
Y(a),M(a),\mathrm{X},N,C\\}=\frac{\\{M(a)\\}!\\{N-M(a)\\}!}{N!}.$

Assumption 1(a) implies that the data vectors from different clusters are
independent, while the outcomes and covariates in the same cluster can be
arbitrarily correlated. Assumption 1(b)-(c) formalize the condition that
$W_{1},\ldots,W_{m}$ are marginally identically distributed according to a
mixture distribution, $\mathcal{P}^{W\mid N}\times\mathcal{P}^{N}$. This
technical condition is useful for handling the varying dimension of $W_{i}$
across clusters. Assumption 2 holds under cluster randomization. Given
Assumptions 1-2, the expectation on $(W_{i},A_{i})$ is taken with respect to
$\mathcal{P}^{W\mid N}\times\mathcal{P}^{N}\times\mathcal{P}^{A}$. Assumption
3 implies that the number of sampled individuals can depend on the assignment
$A$ and cluster characteristics (source population size $N$ and cluster
covariates $C$), but the sampling process is completely random given the
number of sampled individuals and the source population size. This assumption
relaxes the setting of Bugni et al.,  (2022) to allow for arm-specific
sampling. Important special cases of Assumptions 3 include full enrollment (or
focusing on the population of sampled individuals) such that $M(a)=N$; random
sampling with an arm-specific fixed size such that $M(a)=m_{a}$ for some
integer $m_{a}\leq N$; and independent cluster-specific sampling such that
each $S_{ij}(a)$ is independently determined by flipping a cluster-specific
coin. Assumption 3 can be violated when sampling additionally depends on
individual-level covariates, which are generally unobserved for
nonparticipants. This sampling mechanism leads to post-randomization selection
bias (Li et al., , 2022), and we leave this form of selection bias for
separate work.

We define the class of cluster-average treatment effect ($\Delta_{C}$) and
individual-average treatment effect ($\Delta_{I}$) estimands as

$\displaystyle\Delta_{C}=f\left\\{\mu_{C}(1),\mu_{C}(0)\right\\},\quad\Delta_{I}=f\left\\{\mu_{I}(1),\mu_{I}(0)\right\\},$

where $f$ is a pre-specified smooth function determining the scale of effect
measure and, for $a=0,1$,

$\displaystyle\mu_{C}(a)=E\left\\{\frac{\sum_{j=1}^{N_{i}}Y_{ij}(a)}{N_{i}}\right\\},\quad\mu_{I}(a)=\frac{E\left\\{\sum_{j=1}^{N_{i}}Y_{ij}(a)\right\\}}{E(N_{i})}.$

For example, $f(x,y)=x-y$ leads to the difference estimands, $f(x,y)=x/y$
leads to the relative risk estimands, and $f(x,y)=x(1-y)/\\{y(1-x)\\}$ leads
to the odds ratio estimands. The two classes of estimands $\Delta_{C}$ and
$\Delta_{I}$ differ based on the corresponding treatment-specific mean
potential outcomes, $\mu_{C}(a)$ versus $\mu_{I}(a)$. The former represents
the average potential outcome associated with treatment $a$ for a typical
cluster along with its natural source population, while the latter represents
the average potential outcome associated with treatment $a$ for a typical
individual. Depending on the nature of the intervention and the endpoint,
either or both estimands may be relevant in a given cluster-randomized
experiment. These two estimands can be considered as the super-population
analogs to those considered in Su and Ding,  (2021) and Kahan et al.,  (2022)
under the finite-population framework and a generalization of the difference
estimands in Bugni et al.,  (2022). We refer to Kahan et al.,  (2022) for an
in-depth discussion on differentiating these two estimands and their example
applications to cluster-randomized experiments.

Finally, we write the observed data for cluster $i$ as
$\mathcal{O}_{i}=(Y_{i}^{o},M_{i},A_{i},N_{i},C_{i},\mathrm{X}_{i}^{o})$,
where $Y_{i}^{o}=\\{Y_{ij}:S_{ij}=1,j=1,\dots,N_{i}\\}\in\mathbb{R}^{M_{i}}$
and
$\mathrm{X}_{i}^{o}=\\{X_{ij}:S_{ij}=1,j=1,\dots,N_{i}\\}\in\mathbb{R}^{M_{i}\times
p}$. The relationships among these random variables are summarized in Figure
1. The central task is to estimate $\Delta_{C}$ and $\Delta_{I}$ with
$\mathcal{O}_{1},\ldots,\mathcal{O}_{m}$. Of note, we assume that the source
population size $N_{i}$ is available or can be elicited from stakeholders;
this is practically feasible for cluster-randomized experiments conducted
within healthcare delivery systems, schools, or worksites. When the source
population size in each cluster is unknown, we will discuss how to estimate
$\Delta_{C}$ in Section 4, but $\Delta_{I}$ is generally not identifiable
unless we impose $N=M$.

$\mathbf{X}$$\boldsymbol{C}$$N$$A$$\boldsymbol{Y}$$\boldsymbol{Y}^{o}$$M$$\boldsymbol{S}$$\mathbf{X}^{o}$

Figure 1: A graph representation of the causal relationship among random
variables in cluster-randomized experiments. Observed and unobserved nodes
have solid and dashed boundaries, respectively.

##  3 Sufficient conditions for robust inference via model-based methods

###  3.1 Generalized estimating equations

One popular approach to analyze cluster-randomized experiments is through
generalized estimating equations, often specified through the marginal mean
model, $g\\{E(Y_{ij}\mid U_{ij})\\}=U_{ij}^{\top}\beta$, where $g$ is the link
function, $U_{ij}=(1,A_{i},L_{ij}^{\top})^{\top}$ for user-specified
covariates $L_{ij}$ as an arbitrary function of $(N_{i},C_{i},X_{ij})$, and
$\beta=(\beta_{0},\beta_{A},\beta_{L}^{\top})^{\top}$. We assume that $g$ is
the canonical link, e.g., $g(x)=x$ for continuous outcomes and
$g(x)=\log\\{x/(1-x)\\}$ for binary outcomes. The parameters $\beta$ are
estimated by solving the following estimating equations:

$\sum_{i=1}^{m}\mathrm{D}_{i}^{\top}\mathrm{V}_{i}^{-1}(Y_{i}^{o}-\mu_{i}^{o})=0,$
(1)

where $\mu_{i}^{o}=\\{E(Y_{ij}\mid U_{ij}):S_{ij}=1\\}\in\mathbb{R}^{M_{i}}$,
$\mathrm{D}_{i}={\partial\mu_{i}^{o}}/{\partial\beta^{\top}}$,
$\mathrm{V}_{i}=\mathrm{Z}_{i}^{1/2}\mathrm{R}_{i}(\rho)\mathrm{Z}_{i}^{1/2}$
with $\mathrm{R}_{i}(\rho)\in\mathbb{R}^{M_{i}\times M_{i}}$ being a working
correlation matrix and
$\mathrm{Z}_{i}=\mathrm{diag}\\{v(Y_{ij}):S_{ij}=1\\}\in\mathbb{R}^{M_{i}\times
M_{i}}$ for some known variance function $v$. We consider two working
correlation structures that are commonly used for analyzing cluster-randomized
experiments: the independence correlation,
$\mathrm{R}_{i}(\rho)=\mathrm{I}_{M_{i}}$, and the exchangeable correlation,
$\mathrm{R}_{i}(\rho)=(1-\rho)\mathrm{I}_{M_{i}}+\rho
1_{M_{i}}1_{M_{i}}^{\top}$, where
$\mathrm{I}_{M_{i}}\in\mathbb{R}^{M_{i}\times M_{i}}$ is the identity matrix
and $1_{M_{i}}\in\mathbb{R}^{M_{i}}$ is a vector of ones. The estimator for
$\beta$ is denoted by $\widehat{\beta}$, and the correlation parameter $\rho$
is estimated by a moment estimator $\widehat{\rho}$ described in Example 4 of
Liang and Zeger,  (1986). For analyzing cluster-randomized experiments, a
conventional practice is to directly use the coefficient $\widehat{\beta}_{A}$
along with robust sandwich variance for inference. However, this practice can
lead to an ambiguous treatment effect estimand even in the absence of
covariate adjustment (Wang et al., , 2022).

To estimate our target estimands $\Delta_{C}$ and $\Delta_{I}$, we propose a
simple weighted g-computation approach, defined as
$\widehat{\Delta}_{C}^{\textrm{GEE-g}}=f\left\\{\widehat{\mu}_{C}^{\textrm{GEE-g}}(1),\widehat{\mu}_{C}^{\textrm{GEE-g}}(0)\right\\}$
and
$\widehat{\Delta}_{I}^{\textrm{GEE-g}}=f\left\\{\widehat{\mu}_{I}^{\textrm{GEE-g}}(1),\widehat{\mu}_{I}^{\textrm{GEE-g}}(0)\right\\}$,
where, for $a=0,1$,

$\displaystyle\widehat{\mu}_{C}^{\textrm{GEE-g}}(a)$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}\frac{1}{M_{i}}\sum_{j:S_{ij}=1}g^{-1}\left(\widehat{\beta}_{0}+\widehat{\beta}_{A}a+\widehat{\beta}_{L}^{\top}L_{ij}\right),$
$\displaystyle\widehat{\mu}_{I}^{\textrm{GEE-g}}(a)$
$\displaystyle=\frac{1}{\sum_{i=1}^{m}N_{i}}\sum_{i=1}^{m}\frac{N_{i}}{M_{i}}\sum_{j:S_{ij}=1}g^{-1}\left(\widehat{\beta}_{0}+\widehat{\beta}_{A}a+\widehat{\beta}_{L}^{\top}L_{ij}\right).$

The variance of $\widehat{\Delta}_{C}^{\textrm{GEE-g}}$ can be estimated by
the robust sandwich variance estimator, defined as

$\widehat{V}_{C}^{\textrm{GEE-g}}=m^{-2}\nabla\widehat{f}_{C}^{\textrm{GEE-g}}{}^{\top}\left(\sum_{i=1}^{m}\widehat{IF}_{C,i}^{\textrm{GEE-g}}\widehat{IF}_{C,i}^{\textrm{GEE-g}}{}^{\top}\right)\nabla\widehat{f}_{C}^{\textrm{GEE-g}},$

where $\nabla\widehat{f}_{C}^{\textrm{GEE-g}}$ is the gradient of $f$
evaluated at
$\left\\{\widehat{\mu}_{C}^{\textrm{GEE-g}}(1),\widehat{\mu}_{C}^{\textrm{GEE-g}}(0)\right\\}$
and

$\widehat{IF}_{C,i}^{\textrm{GEE-g}}=\left[\begin{array}[]{c}\displaystyle\left\\{\frac{1}{m}\sum_{l=1}^{m}\frac{A_{l}M_{l}}{1+(M_{l}-1)\widehat{\rho}}\right\\}^{-1}\frac{A_{i}}{1+(M_{i}-1)\widehat{\rho}}\\{\overline{Y}_{i}^{o}-\overline{\mu}_{i}(1,\widehat{\beta})\\}+\overline{\mu}_{i}(1,\widehat{\beta})-\widehat{\mu}_{C}^{\textrm{GEE-g}}(1)\\\
\displaystyle\left\\{\frac{1}{m}\sum_{l=1}^{m}\frac{(1-A_{l})M_{l}}{1+(M_{l}-1)\widehat{\rho}}\right\\}^{-1}\frac{1-A_{i}}{1+(M_{i}-1)\widehat{\rho}}\\{\overline{Y}_{i}^{o}-\overline{\mu}_{i}(0,\widehat{\beta})\\}+\overline{\mu}_{i}(0,\widehat{\beta})-\widehat{\mu}_{C}^{\textrm{GEE-g}}(0)\end{array}\right],$

is the estimated influence function for cluster $i$,
$\overline{Y}_{i}^{o}=M_{i}^{-1}S_{i}^{\top}Y_{i}$, and
$\overline{\mu}_{i}(a,\widehat{\beta})=M_{i}^{-1}S_{i}^{\top}\mu_{i}(a,\widehat{\beta})$
with $\mu_{i}(a,\widehat{\beta})$ being the estimate of
$\\{g^{-1}(U_{ij}^{\top}\beta):j=1,\dots,N_{i}\\}$ given $A=a$ for
$a\in\\{0,1\\}$. Similarly, the variance of
$\widehat{\Delta}_{I}^{\textrm{GEE-g}}$ is estimated by

$\widehat{V}_{I}^{\textrm{GEE-g}}=(\sum_{i=1}^{m}N_{i})^{-2}\nabla\widehat{f}_{I}^{\textrm{GEE-g}}{}^{\top}\left(\sum_{i=1}^{m}N_{i}^{2}\widehat{IF}_{I,i}^{\textrm{GEE-g}}\widehat{IF}_{I,i}^{\textrm{GEE-g}}{}^{\top}\right)\nabla\widehat{f}_{I}^{\textrm{GEE-g}}$

where $\widehat{f}_{I}^{\textrm{GEE-g}}$ and
$\widehat{IF}_{I,i}^{\textrm{GEE-g}}$ are defined in the same way as
$\widehat{f}_{C}^{\textrm{GEE-g}}$ and $\widehat{IF}_{C,i}^{\textrm{GEE-g}}$
except that the subscript is substituted by $I$ in all estimated parameters.

To proceed, we additionally make the following assumption, which implies that
the observed cluster size within each arm has only exogenous randomness.

######  Assumption 4 (Arm-specific random observed cluster sizes).

For $a\in\\{0,1\\}$, $M(a)$ is independent of $\\{N,C,\mathrm{X},Y(a)\\}$.

Theorem 1 below articulates several model specifications of generalized
estimating equations such that $\widehat{\Delta}_{I}^{\textrm{GEE-g}}$ and
$\widehat{\Delta}_{C}^{\textrm{GEE-g}}$ are consistent and asymptotically
normal, leading to valid statistical inference for well-defined estimands. The
regularity conditions required are moment and continuity assumptions similar
to those invoked in Theorem 5.31 of van der Vaart,  (1998) for M-estimators.
Particularly, model specifications S2-S4 in Theorem 1 make no assumption on
the underlying distribution for the potential outcomes, and hence allow for
model-robust estimation of both estimands.

######  Theorem 1.

(a) Under Assumptions 1-4 and regularity conditions in the Supplementary
Material,
$\left(\widehat{V}_{C}^{\textrm{GEE-g}}\right)^{-1/2}(\widehat{\Delta}_{C}^{\textrm{GEE-g}}-{\Delta}_{C})\xrightarrow{d}\mathcal{N}(0,1)$
if any of the following conditions holds: (S1) the mean model
$g\\{E(Y_{ij}\mid U_{ij})\\}=U_{ij}^{\top}\beta$ is correctly specified; (S2)
an independence working correlation structure is used; (S3) $g$ is the
identity link function and the working variance is constant with $v(Y_{ij}\mid
U_{ij})=\sigma^{2}$; (S4) $L_{ij}$ is a function of cluster-level covariates
$(N_{i},C_{i})$ but involves no individual-level covariate. (b) If the
estimating equations (1) are additionally weighted by the source population
size $N_{i}$ for each cluster $i$, then
$\left(\widehat{V}_{I}^{\textrm{GEE-g}}\right)^{-1/2}(\widehat{\Delta}_{I}^{\textrm{GEE-g}}-{\Delta}_{I})\xrightarrow{d}\mathcal{N}(0,1)$
if any of the above conditions (S1-S4) holds.

######  Remark 1.

Assumption 4 can be relaxed under model specification (S1) or (S2) if each
cluster is further weighted by $1/{M_{i}}$ in Equations (1) before applying
the weighted g-computation formulas. In other cases, Assumption 4 is generally
required for valid inference of the g-computation estimators.

###  3.2 (Generalized) linear mixed models

Generalized linear mixed models are another popular method for analyzing
cluster-randomized experiments. If we write $b_{i}$ as the random effect for
cluster $i$, a typical generalized linear mixed model applied to cluster-
randomized experiment often includes the following assumptions: (i)
$b_{1},\dots,b_{m}$ are independent, identically distributed from
$\mathcal{N}(0,\tau^{2})$, where $\tau^{2}$ is an unknown variance component;
(ii) $Y_{i1},\dots,Y_{iN_{i}}$ are conditionally independent given
$(U_{i1},\dots,U_{iN_{i}},b_{i})$; and (iii) the distribution
$Y_{ij}\mid(U_{i1},\dots,U_{iN_{i}},b_{i})$ is a member of the exponential
family with $E(Y_{ij}\mid U_{i1},\dots,U_{iN_{i}},b_{i})=E(Y_{ij}\mid
U_{ij},b_{i})=g^{-1}(U_{ij}^{\top}\alpha+b_{i})$, where $g$ is the canonical
link and $\alpha=(\alpha_{0},\alpha_{A},\alpha_{L}^{\top})^{\top}$. Given the
above assumptions, the marginal likelihood of observed data can be written as
$\prod_{i=1}^{m}\int\prod_{j:S_{ij}=1}p(Y_{ij}\mid
U_{ij},b_{i};\alpha)\phi(b_{i};\tau^{2})\textrm{d}b_{i}$, where $p(Y_{ij}\mid
U_{ij},b_{i};\alpha)$ is the conditional density of $Y_{ij}$ and
$\phi(b_{i};\tau^{2})$ is the normal density of $b_{i}$. In general, the
estimation of model parameters proceeds by maximizing this likelihood
function, and a common practice is to consider $\alpha_{A}$ as the treatment
effect parameter. However, the interpretation of $\alpha$ is conditional on
random effects, and if the mixed model is misspecified, $\alpha$ lacks direct
connection to our marginal estimands, with an important exception that we
detail below.

An interesting special case where the mixed model provides model-robust
inference in cluster-randomized experiments is when a linear mixed model is
considered as the working model. In this case, we define the weighted
g-computation estimators as
$\widehat{\Delta}_{C}^{\textrm{LMM-g}}=f\left\\{\widehat{\mu}_{C}^{\textrm{LMM-g}}(1),\widehat{\mu}_{C}^{\textrm{LMM-g}}(0)\right\\}$
and
$\widehat{\Delta}_{I}^{\textrm{LMM-g}}=f\left\\{\widehat{\mu}_{I}^{\textrm{LMM-g}}(1),\widehat{\mu}_{I}^{\textrm{LMM-g}}(0)\right\\}$,
where

$\displaystyle\widehat{\mu}_{C}^{\textrm{LMM-g}}(a)=$
$\displaystyle\widehat{\alpha}_{0}+\widehat{\alpha}_{A}a+\widehat{\alpha}_{L}^{\top}\frac{\sum_{i=1}^{m}\overline{L}^{o}_{i}}{m},$
$\displaystyle\widehat{\mu}_{I}^{\textrm{LMM-g}}(a)=$
$\displaystyle\widehat{\alpha}_{0}+\widehat{\alpha}_{A}a+\widehat{\alpha}_{L}^{\top}\frac{\sum_{i=1}^{m}N_{i}\overline{L}^{o}_{i}}{\sum_{i=1}^{m}N_{i}}$

with $\overline{L}^{o}_{i}={M_{i}}^{-1}\sum_{j=1}^{N_{i}}S_{ij}L_{ij}$ being
the average covariate value for cluster $i$ among the observed individuals.
Furthermore, when the interest lies in the difference estimands with
$f(x,y)=x-y$, we simply have
$\widehat{\Delta}_{C}^{\textrm{LMM-g}}=\widehat{\alpha}_{A}$ and
$\widehat{\Delta}_{I}^{\textrm{LMM-g}}=\widehat{\alpha}_{A}$ such that the
treatment effect estimator from the linear mixed model can be taken as the
average causal effect estimator, but these two coefficients are estimated via
different weights applied to the log-likelihood as we explain in Theorem 2.
The variance estimate of $\widehat{\Delta}_{C}^{\textrm{LMM-g}}$, denoted as
$\widehat{V}_{C}^{\textrm{LMM-g}}$, is constructed in the same way as
$\widehat{V}_{C}^{\textrm{GEE-g}}$ with $g$ being the identity link and
$\\{\widehat{\mu}_{C}^{\textrm{GEE-g}}(1),\widehat{\mu}_{C}^{\textrm{GEE-g}}(0),\widehat{\beta},\widehat{\rho}\\}$
substituted by
$\\{\widehat{\mu}_{C}^{\textrm{LMM-g}}(1),\widehat{\mu}_{C}^{\textrm{LMM-g}}(0),\widehat{\alpha},\widehat{\tau}^{2}(\widehat{\tau}^{2}+\widehat{\sigma}^{2})^{-1}\\}$;
for $\widehat{\Delta}_{I}^{\textrm{LMM-g}}$, the variance estimate
$\widehat{V}_{I}^{\textrm{LMM-g}}$ can be constructed in a similar way.

Theorem 2 shows that $\widehat{\Delta}_{C}^{\textrm{LMM-g}}$ and
$\widehat{\Delta}_{I}^{\textrm{LMM-g}}$ are asymptotically valid if the
observed cluster size is only subject to exogenous variation within each arm,
even when the linear mixed model is arbitrarily misspecified.

######  Theorem 2.

Under Assumptions 1-4 and regularity conditions provided in the Supplementary
Material,
$\left(\widehat{V}_{C}^{\textrm{LMM-g}}\right)^{-1/2}(\widehat{\Delta}_{C}^{\textrm{LMM-g}}-{\Delta}_{C})\xrightarrow{d}\mathcal{N}(0,1)$.
If each cluster is weighted by $N_{i}$ in the log-likelihood function of the
working linear mixed model, then
$\left(\widehat{V}_{I}^{\textrm{GEE-g}}\right)^{-1/2}(\widehat{\Delta}_{I}^{\textrm{LMM-g}}-{\Delta}_{I})\xrightarrow{d}\mathcal{N}(0,1)$.

Theorem 2 can be considered as the counterpart of Theorem 1 for linear mixed
models in cluster-randomized experiments. Compared to generalized estimating
equations with an identity link function and working exchangeable covariance
structure, linear mixed models yield similar estimating equations, but the
estimators for nuisance parameters are different, i.e., moment estimation for
$\rho$ in generalized estimating equations versus maximum likelihood
estimation for $\sigma^{2}$, $\tau^{2}$ in linear mixed models, leading to
slightly different asymptotic variances. In the special case that
$M_{i}(1)=M_{i}(0)$, $N_{i}$ is constant, $f(x,y)=x-y$, and $Y_{ij}$ are
marginally identically distributed, Theorem 2 reduces to Theorem 1 of Wang et
al.,  (2021).

Finally, while generalized estimating equations and generalized linear mixed
models are the most widely used approach for analyzing cluster-randomized
experiments, there are two alternative estimators, the augmented generalized
estimating equations (Stephens et al., , 2012) and targeted maximum likelihood
estimation (Balzer et al., , 2021; Benitez et al., , 2021) that can also
provide robust estimation under Assumptions 1-4 when certain aspects of
working models are misspecified. We provide discussions of those approaches in
the Supplementary Material.

##  4 Efficient estimation

###  4.1 Efficient influence functions

For model-based methods to provide model-robust inference, Theorems 1 and 2
require the observed cluster size to be independent of cluster characteristics
$(N,C)$. To address this limitation, we develop more principled estimators to
simultaneously optimize covariate adjustment and address variable observed
cluster sizes arising from post-randomization cluster-dependent sampling
schemes. The proposed estimators are motivated by the efficient influence
function, which is a non-parametric functional of observed data that
characterizes the target estimand (Hines et al., , 2022). With the efficient
influence function, one can derive the efficiency bound for an estimand, i.e.,
the asymptotic variance of all regular and asymptotically linear estimators
over the underlying causal model is lower bounded by the variance of the
efficient influence function. More importantly, recent advances in causal
inference (e.g. van der Laan et al., , 2011; Chernozhukov et al., , 2018)
showed how to use the efficient influence functions to construct an efficient
estimator, i.e., achieving the efficiency bound, by incorporating machine
learning algorithms. Theorem 3 provides the efficient influence functions for
$\mu_{C}(a)$ and $\mu_{I}(a)$, from which the efficient influence functions
for $\Delta_{C}$ and $\Delta_{I}$ can be obtained by the chain rule.

######  Theorem 3.

(a) Given Assumptions 1-3, the efficient influence function for $\mu_{C}(a)$
is

$\displaystyle\text{EIF}_{C}(a)$
$\displaystyle=\frac{I\\{A=a\\}}{\pi^{a}(1-\pi)^{1-a}}\left\\{\overline{Y}^{o}-E\left(\overline{Y}^{o}\mid
A=a,\mathrm{X}^{o},M,N,C\right)\right\\}$ $\displaystyle\quad+\frac{pr(A=a\mid
M,N,C)}{\pi^{a}(1-\pi)^{1-a}}\left\\{E\left(\overline{Y}^{o}\mid
A=a,\mathrm{X}^{o},M,N,C\right)-E\left(\overline{Y}^{o}\mid
A=a,N,C\right)\right\\}$ $\displaystyle\quad+E\left(\overline{Y}^{o}\mid
A=a,N,C\right)-\mu_{C}(a),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ a\in\\{0,1\\},$

where $I\\{A=a\\}$ is an indicator function of $A=a$ and
$\overline{Y}^{o}=M^{-1}\sum_{j=1}^{N}S_{\cdot j}Y_{\cdot j}$ refers to the
average outcome among the observed participants in each cluster. (b) Given
Assumptions 1-3, the efficient influence function for $\mu_{I}(a)$ is

$\displaystyle\text{EIF}_{I}(a)=\frac{N}{E(N)}\left\\{\text{EIF}_{C}(a)+\mu_{C}(a)-\mu_{I}(a)\right\\},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ a\in\\{0,1\\}.$

The efficient influence functions in Theorem 3 involve three nuisance
functions, which we denote by $\eta_{a}^{*}=E\left(\overline{Y}^{o}\mid
A=a,\mathrm{X}^{o},M,N,C\right)$, $\zeta_{a}^{*}=E\left(\overline{Y}^{o}\mid
A=a,N,C\right)$, and $\kappa_{a}^{*}=pr(A=a\mid M,N,C)$. Compared to
$\text{EIF}_{C}(a)$, $\text{EIF}_{I}(a)$ additionally includes a weight,
$N/E(N)$, to target the individual-average treatment effect. When $M$ is
identical to the source population size $N$, the number of sampled individuals
can be treated as a pre-randomization variable, and the efficient influence
function for $\mu_{C}(a)$ reduces to that in Balzer et al.,  (2019), where the
only nuisance function is $\eta_{a}^{*}$. In more general settings as we
consider here, two additional nuisance functions, $\zeta_{a}^{*}$ and
$\kappa_{a}^{*}$, are needed to leverage the post-randomization information
$\mathrm{X}^{o}$ for improving efficiency. Next, we construct new estimators
for $\Delta_{C}$ and $\Delta_{I}$ based on the efficient influence functions.

###  4.2 Efficient estimators based on the efficient influence functions

Based on Theorem 3, we propose the following estimator for $\mu_{C}(a)$ and
$\mu_{I}(a)$:

$\displaystyle\widehat{\mu}_{C}^{\textrm{Eff}}(a)$
$\displaystyle=\frac{1}{m}\sum_{i=1}^{m}D_{i}(a,\widehat{h}_{a}),$
$\displaystyle\widehat{\mu}_{I}^{\textrm{Eff}}(a)$
$\displaystyle=\frac{1}{\sum_{i=1}^{m}N_{i}}\sum_{i=1}^{m}N_{i}D_{i}(a,\widehat{h}_{a}),$

where
$\widehat{h}_{a}=(\widehat{\eta}_{a},\widehat{\zeta}_{a},\widehat{\kappa}_{a})$
are user-specified estimators for nuisance functions
$h^{*}_{a}=(\eta_{a}^{*},\zeta_{a}^{*},\kappa_{a}^{*})$ and

$\displaystyle D_{i}(a,\widehat{h}_{a})$
$\displaystyle=\frac{I\\{A_{i}=a\\}}{\pi^{a}(1-\pi)^{1-a}}\left\\{\overline{Y}_{i}^{o}-\widehat{\eta}_{a}(\mathrm{X}_{i}^{o},M_{i},N_{i},C_{i})\right\\}$
$\displaystyle\quad+\frac{\widehat{\kappa}_{a}(M_{i},N_{i},C_{i})}{\pi^{a}(1-\pi)^{1-a}}\left\\{\widehat{\eta}_{a}(\mathrm{X}_{i}^{o},M_{i},N_{i},C_{i})-\widehat{\zeta}_{a}(N_{i},C_{i})\right\\}+\widehat{\zeta}_{a}(N_{i},C_{i}).$

Then the target estimands defined in Section 2 are estimated by
$\widehat{\Delta}_{C}^{\textrm{Eff}}=f\left\\{\widehat{\mu}_{C}^{\textrm{Eff}}(1),\widehat{\mu}_{C}^{\textrm{Eff}}(0)\right\\}$
and
$\widehat{\Delta}_{I}^{\textrm{Eff}}=f\left\\{\widehat{\mu}_{I}^{\textrm{Eff}}(1),\widehat{\mu}_{I}^{\textrm{Eff}}(0)\right\\}$.
Among the many possibilities for estimating these nuisance functions, we
primarily consider the following two approaches.

The first approach considers parametric working models, that are,
$\widehat{\eta}_{a}=\eta_{a}(\widehat{\theta}_{\eta,a}),\widehat{\zeta}_{a}=\zeta_{a}(\widehat{\theta}_{\zeta,a}),\widehat{\kappa}_{a}=\kappa_{a}(\widehat{\theta}_{\kappa,a})$
for pre-specified functions $\eta_{a},\zeta_{a},\kappa_{a}$ with finite-
dimensional parameters $\theta_{\eta,a},\theta_{\zeta,a},\theta_{\kappa,a}$.
In this case, we use superscript “Eff-PM”, e.g.,
$\widehat{\Delta}_{C}^{\textrm{Eff-PM}}$, to highlight the role of parametric
models in estimation. Example working models include generalized estimating
equations, generalized linear mixed models, penalized regression for variable
selection, among others; in all cases, the associated parameters are estimated
by solving estimating equations, and we assume that the estimating equations
satisfy regularity conditions provided in the Supplementary Material such that
the nuisance parameter estimators are asymptotically linear. With this method,
the parameters
$\theta=\\{\mu_{C}(a),\theta_{\eta,a},\theta_{\zeta,a},\theta_{\kappa,a}:a=0,1\\}$,
can be viewed as jointly estimated by solving estimating equations
$\sum_{i=1}^{m}\psi_{C}(\mathcal{O}_{i};\theta)=0$ for a known function
$\psi_{C}$, and the sandwich variance estimator $\widehat{V}_{C}^{\textrm{Eff-
PM}}$ for $\widehat{\Delta}_{C}^{\textrm{Eff-PM}}$ is computed as

$\nabla\widehat{f}_{C}^{\top}\left\\{\sum_{i=1}^{m}\frac{\partial}{\partial\theta}\psi_{C}(\mathcal{O}_{i};\theta)|_{\theta=\widehat{\theta}}\right\\}^{-1}\left(\sum_{i=1}^{m}\psi_{C}(\mathcal{O}_{i};\widehat{\theta})\psi_{C}(\mathcal{O}_{i};\widehat{\theta}){}^{\top}\right)\left\\{\sum_{i=1}^{m}\frac{\partial}{\partial\theta}\psi_{C}(\mathcal{O}_{i};\theta)|_{\theta=\widehat{\theta}}\right\\}^{-1}\nabla\widehat{f}_{C},$

where $\nabla\widehat{f}_{C}$ is the gradient of $f$ over $\theta$ evaluated
at $\widehat{\theta}$. Following a similar procedure we can compute the
sandwich variance estimator $\widehat{V}_{I}^{\textrm{Eff-PM}}$ for
$\widehat{\Delta}_{I}^{\textrm{Eff-PM}}$. For the subsequent technical
discussions, we denote the probability limit of
$(\widehat{\theta}_{\eta,a},\widehat{\theta}_{\zeta,a},\widehat{\theta}_{\kappa,a})$
as
$(\underline{\theta}_{\eta,a},\underline{\theta}_{\zeta,a},\underline{\theta}_{\kappa,a})$.

The second approach uses machine learning algorithms with cross-fitting to
estimate all nuisance functions. For this case, we use superscript “Eff-ML”,
e.g., $\widehat{\Delta}_{C}^{\textrm{Eff-ML}}$, to indicate using machine
learning algorithms. We assume that each nuisance function estimator is
consistent to the truth at an $m^{{1}/{4}}$ rate such that
$\widehat{h}_{a}-h_{a}^{*}=o_{p}(m^{-1/4})$ in $L_{2}(\mathcal{P})$-norm. This
rate can be achieved by many methods such as random forests (Wager and
Walther, , 2015), neural networks (Farrell et al., , 2021), and boosting (Luo
and Spindler, , 2016). In addition, we assume a regularity condition that
$\widehat{\kappa}_{a}$ and $E\\{(\eta_{a}^{*})^{2}\mid M,N,C\\}$ are uniformly
bounded, a similar condition invoked in Chernozhukov et al.,  (2018) for
controlling the remainder term and consistently estimating the variance. For
cross-fitting, we randomly partition $m$ clusters into $K$ parts with roughly
equal sizes (the size difference is at most 1), and denote
$\mathcal{O}^{*}_{k}$ as the $k$-th part and
$\mathcal{O}^{*}_{-k}=\bigcup_{k^{\prime}\in\\{1,\ldots,K\\}\setminus\\{k\\}}\mathcal{O}^{*}_{k^{\prime}}$.
For each $k$, we compute $\widehat{h}_{a,k}$, which is the nuisance function
trained on $\mathcal{O}^{*}_{-k}$ and evaluated at $\mathcal{O}^{*}_{k}$, and
the estimated nuisance function $\widehat{h}_{a}$ evaluated on all clusters is
then the combination of $\widehat{h}_{a,k}$ for all $k=1,\dots,K$. In
practice, we recommend choosing $K$ such that $m/K\geq 10$. For variance
estimation, we propose the following consistent estimators

$\displaystyle\widehat{V}_{C}^{\textrm{Eff-ML}}$
$\displaystyle=\frac{1}{m^{2}}\sum_{k=1}^{K}\sum_{i\in\mathcal{O}_{k}}\left[\sum_{a=0}^{1}\dot{f}_{a,C}^{\textrm{Eff-
ML}}\left\\{D_{i}(\widehat{h}_{a,k})-\frac{1}{|\mathcal{O}_{k}|}\sum_{l\in\mathcal{O}_{k}}D_{l}(\widehat{h}_{a,k})\right\\}\right]^{2}$
$\displaystyle\widehat{V}_{I}^{\textrm{Eff-ML}}$
$\displaystyle=\frac{1}{(\sum_{i=1}^{m}N_{i})^{2}}\sum_{k=1}^{K}\sum_{i\in\mathcal{O}_{k}}\left[\sum_{a=0}^{1}\dot{f}_{a,I}^{\textrm{Eff-
ML}}\left\\{N_{i}D_{i}(\widehat{h}_{a,k})-\frac{1}{|\mathcal{O}_{k}|}\sum_{l\in\mathcal{O}_{k}}N_{l}D_{l}(\widehat{h}_{a,k})\right\\}\right]^{2}$

where $\dot{f}_{a,C}^{\textrm{Eff-ML}}$ is the partial derivative of
$\mu_{C}(a)$ evaluated at {$\widehat{\mu}_{C}^{\textrm{Eff-
ML}}(1),\widehat{\mu}_{C}^{\textrm{Eff-ML}}(0)\\}$,
$\dot{f}_{a,I}^{\textrm{Eff-ML}}$ is the partial derivative of $\mu_{I}(a)$
evaluated at {$\widehat{\mu}_{I}^{\textrm{Eff-
ML}}(1),\widehat{\mu}_{I}^{\textrm{Eff-ML}}(0)\\}$, and $|\mathcal{O}_{k}|$ is
the size of part $k$ in crossing fitting.

Theorem 4 summarizes the asymptotic behaviors of our proposed estimators under
both strategies for estimating the nuisance functions.

######  Theorem 4.

Given Assumptions 1-3 and above regularity conditions,

  1. (a)

if $\kappa_{a}({\underline{\theta}}_{\kappa,a})=\kappa_{a}^{*}$ or
$E\\{\eta_{a}({\underline{\theta}}_{\eta,a})\mid
M,N,C\\}=\zeta_{a}({\underline{\theta}}_{\zeta,a})$, then
$\left(\widehat{V}_{C}^{\textrm{Eff-
PM}}\right)^{-1/2}(\widehat{\Delta}_{C}^{\textrm{Eff-
PM}}-\Delta_{C})\xrightarrow{d}\mathcal{N}(0,1)$ and
$\left(\widehat{V}_{I}^{\textrm{Eff-
PM}}\right)^{-1/2}(\widehat{\Delta}_{I}^{\textrm{Eff-
PM}}-\Delta_{I})\xrightarrow{d}\mathcal{N}(0,1)$;

  2. (b)

if $\widehat{h}_{a}-h_{a}^{*}=o_{p}(m^{-1/4})$ in $L_{2}(\mathcal{P})$-norm,
then $\left(\widehat{V}_{C}^{\textrm{Eff-
ML}}\right)^{-1/2}(\widehat{\Delta}_{C}^{\textrm{Eff-
ML}}-\Delta_{C})\xrightarrow{d}\mathcal{N}(0,1)$ and
$\left(\widehat{V}_{I}^{\textrm{Eff-
ML}}\right)^{-1/2}(\widehat{\Delta}_{I}^{\textrm{Eff-
ML}}-\Delta_{I})\xrightarrow{d}\mathcal{N}(0,1)$. Furthermore,
$\widehat{V}_{C}^{\textrm{Eff-ML}}$ and $\widehat{V}_{I}^{\textrm{Eff-ML}}$
achieve the efficiency lower bounds asymptotically.

For the proposed estimators with parametric working models, Theorem 4 implies
that the estimator is consistent if $\kappa_{a}$ is correctly specified, or
the working models $\eta_{a}$ and $\zeta_{a}$ are equivalent conditioning on
$M,N,C$. In the special case that all cluster-level covariates are discrete,
the latter condition can be satisfied by setting

$\eta_{a}(\underline{\theta}_{\eta,a})=\zeta_{a}+\underline{\beta}^{\top}\left(\overline{X}^{o}-\sum_{n}\sum_{c}I\\{N=n,C=c\\}\underline{\theta}_{n,c}\right),$

where $\underline{\beta}$ is an arbitrary $p$-dimensional vector and
$\underline{\theta}_{n,c}=E\\{\overline{X}^{o}\mid N=n,C=c\\}$; in this
special case, the resulting estimators are in fact robust to arbitrary working
model misspecification. In more general cases with non-discrete cluster-level
covariates, the proposed estimators are at least triply-robust. That is, they
are consistent to their respective target estimands as long as two out of the
three nuisance functions are correctly modeled, since $E(\eta^{*}_{a}\mid
M,N,C)=\zeta^{*}_{a}$ as proved in the Supplementary Material.

For the proposed estimators using machine learning algorithms, efficiency can
be achieved when all nuisance parameters are consistently estimated, leading
to higher asymptotic precision than potentially misspecified parametric
working models. For modeling $\eta_{a}^{*}$, since $\mathrm{X}^{o}$ is a
matrix and its dimension may change across clusters, a feasible practical
strategy is to fit $Y_{ij}$ on $(X_{ij},M_{i},N_{i},C_{i})$ and pre-specified
summary statistics of $\mathrm{X}_{i}^{o}$ with fixed dimensions, e.g.,
$\overline{X}_{i}^{o}$, and then compute the cluster-average of predictions as
the model fit. Alternatively, one can directly model $\overline{Y}_{i}^{o}$ on
$(M_{i},N_{i},C_{i})$ and functions of $\mathrm{X}_{i}^{o}$, which can be
potentially high-dimensional in order to capture its higher-order association
with $\overline{Y}_{i}^{o}$.

######  Remark 2.

In practice, each of the two approaches for nuisance function estimation has
its pros and cons. The machine learning methods are asymptotically more
precise, but the efficiency gain typically requires at least a moderate number
of clusters. In addition, although the cross-fitting procedure yields the
desired convergence property, it may be prone to finite-sample bias especially
when the sample size is limited (Hines et al., , 2022). This finite-sample
bias can be potentially alleviated by using parsimonious parametric working
models, thereby trading precision for better finite-sample performance
characteristics.

######  Remark 3.

If the observed cluster size $M<N$, our proposed estimators require an
accurate estimate of the source population size for efficient estimation. If
$N$ is not fully available for all clusters, one can still use the observed
data $(\overline{Y}_{i}^{o},A_{i},C_{i})$ to infer $\Delta_{C}$, and the
efficient influence function for $\mu_{C}(a)$ becomes

$\frac{I\\{A=a\\}}{\pi^{a}(1-\pi)^{1-a}}(\overline{Y}^{o}-\zeta_{a}^{*})+\zeta_{a}^{*}-\mu_{C}(a),$

based on which an estimator for $\Delta_{C}$ can be constructed. For example,
we can apply our proposed estimator with $\widehat{\eta}_{a}$ set to be equal
to $\widehat{\zeta}_{a}$ and get a consistent estimator for $\Delta_{C}$ even
when $\widehat{\zeta}_{a}$ is incorrectly specified. The individual-average
treatment effect $\Delta_{I}$, however, is generally not identifiable without
observing $N$.

##  5 Simulation experiments

###  5.1 Simulation design

We conducted two simulation experiments to compare different methods for
analyzing cluster-randomized experiments. The first simulation study focused
on estimating the difference estimands of the cluster-average and individual-
average treatment effect for continuous outcomes, while the second study
focused on the relative risk estimands for binary outcomes. In each
experiment, we considered a small, $m=30$, or large, $m=100$, number of
clusters, and random observed cluster sizes, i.e., $M_{i}$ is independent of
other variables as a special case of Assumption 4, or cluster-dependent
observed cluster sizes, i.e., $M_{i}$ depends on treatment and cluster-level
covariates. Combinations of these specifications are labeled as scenarios 1-4
in Table 1 and Table 2.

In the first simulation experiment, we let $(N_{i},C_{i1},C_{i2})$,
$i=1,\dots,m$ be independent draws from distribution
$\mathcal{P}^{N}\times\mathcal{P}^{C_{1}\mid N}\times\mathcal{P}^{C_{2}\mid
N,C_{1}}$, where $\mathcal{P}^{N}$ is uniform over support $\\{10,50\\}$,
$\mathcal{P}^{C_{1}\mid N}=\mathcal{N}(N/10,4)$, and $\mathcal{P}^{C_{2}\mid
N,C_{1}}=\mathcal{B}[\textrm{expit}\\{\log(N/10)C_{1}\\}]$ is a Bernoulli
distribution with $\textrm{expit}(x)=(1+e^{-x})^{-1}$. Next, for each
individual in the source population, we generated the individual-level
covariates from $X_{ij1}\sim\mathcal{B}\left({N_{i}}/{50}\right)$,
$X_{ij2}\sim\mathcal{N}\left\\{{\sum_{j=1}^{N_{i}}X_{ij1}}(2C_{i2}-1)/{N_{i}},9\right\\}$,
and potential outcomes from
$Y_{ij}(1)\sim\mathcal{N}\left\\{N_{i}/5+N_{i}\sin(C_{i1})(2C_{i2}-1)/30+5e^{X_{ij1}}\mid
X_{ij2}\mid,1\right\\}$ and
$Y_{ij}(0)\sim\mathcal{N}\left\\{\gamma_{i}+N_{i}\sin(C_{i1})(2C_{i2}-1)/30+5e^{X_{ij1}}\mid
X_{ij2}\mid,1\right\\}$ where $\gamma_{i}\sim\mathcal{N}(0,1)$ is a cluster-
level random intercept to allow for residual intracluster correlations.
Furthermore, we set $M_{i}(1)=M_{i}(0)=9+\mathcal{B}(0.5)$ for the random
observed cluster size scenario, and $M_{i}(1)=N_{i}/5+5C_{i2}$,
$M_{i}(0)=3I\\{N_{i}=50\\}+3$ for the cluster-dependent observed cluster size
scenario. Then, we independently sampled $A_{i}\sim\mathcal{B}(0.5)$ and
defined $Y_{ij}=A_{i}Y_{ij}(1)+(1-A_{i})Y_{ij}(0)$ and
$M_{i}=A_{i}M_{i}(1)+(1-A_{i})M_{i}(0)$. Finally, for each cluster, we
uniformly sampled without replacement $M_{i}$ individuals, for whom
$S_{ij}=1$, and the observed data in each simulation replicate are
$\\{Y_{ij},A_{i},M_{i},C_{i1},C_{i2},X_{ij1},X_{ij2}:S_{ij}=1,i=1,\dots,m,j=1,\dots,N_{i}\\}$.
For the second simulation study, the data were generated following the first
simulation study, except that the potential outcomes were drawn from the
following Bernoulli distributions:
$Y_{ij}(1)\sim\mathcal{B}\left[\textrm{expit}\left\\{-N_{i}/20+N_{i}\sin(C_{i1})(2C_{i2}-1)/30+1.5e^{X_{ij1}}\mid
X_{ij2}\mid^{1/2}\right\\}\right]$ and
$Y_{ij}(0)\sim\mathcal{B}\left[\textrm{expit}\left\\{\gamma_{i}+N_{i}\sin(C_{i1})(2C_{i2}-1)/30+1.5(2X_{ij1}-1)\mid
X_{ij2}\mid^{1/2}\right\\}\right]$.

We compared the following methods. The unadjusted method (Bugni et al., ,
2022) is equivalent to our proposed method setting
$\widehat{\eta}_{a}=\widehat{\zeta}_{a}$ to be a constant. Generalized
estimating equations with weighted g-computation were implemented as described
in Section 3.1 with an exchangeable working correlation for continuous
outcomes and an independence working correlation for binary outcomes. Linear
mixed models with weighted g-computation were implied as described in Section
3.2 for both continuous and binary outcomes. For our proposed methods, we
considered both parametric working models and machine learning algorithms,
where the former used generalized linear models for all nuisance functions,
and the latter exploited `SuperLearner` (van der Laan et al., , 2007) for
model-fitting with generalized linear models, regression trees, and neural
networks. Across the comparators, individual-level regression working models
adjusted for covariates $(N_{i},C_{i1},C_{i2},X_{ij1},X_{ij2})$, and cluster-
level working models adjusted for covariates $(N_{i},C_{i1},C_{i2})$. We used
the proposed variance estimator with the degree-of-freedom adjustment to
improve the small-sample performance. For each scenario of both simulation
experiments, we randomly generated $10,000$ data sets and tested the above
methods on each data set.

###  5.2 Simulation Results

Table 1 summarizes the simulation results for continuous outcomes. Since the
outcome distribution varies by the source population size, the cluster-average
treatment effect $\Delta_{C}=6$, which differs from the individual-average
treatment effect $\Delta_{I}=8.67$. Across all scenarios, the proposed methods
with no covariate adjustment, parametric working models, or machine learning
algorithms have negligible bias and nominal coverage, while the model-based
methods, i.e., generalized estimating equations and linear mixed models, show
bias and under-coverage if the observed cluster sizes are cluster-dependent.

Table 1: Results in the first simulation experiment with continuous outcomes.

|  |  Cluster-average treatment effect $\Delta_{C}=6$ |  |  Individual-average treatment effect $\Delta_{I}=8.67$
---|---|---|---|---
Setting | Method | Bias | ESE | ASE | CP |  | Bias | ESE | ASE | CP
Scenario 1: Small $m$ with random observed cluster sizes | Unadjusted | $-0.08$ | $4.85$ | $4.88$ | $0.95$ |  | $-0.15$ | $4.60$ | $4.28$ | $0.93$
GEE-g | $0.13$ | $2.77$ | $2.45$ | $0.92$ |  | $-0.13$ | $3.64$ | $3.10$ | $0.89$
LMM-g | $0.12$ | $2.77$ | $2.78$ | $0.96$ |  | $-0.13$ | $3.63$ | $2.94$ | $0.90$
Eff-PM | $-0.03$ | $3.23$ | $2.66$ | $0.92$ |  | $-0.10$ | $3.89$ | $3.38$ | $0.91$
Eff-ML | $0.03$ | $2.24$ | $2.05$ | $0.95$ |  | $0.03$ | $2.59$ | $2.49$ | $0.95$
Scenario 2: Small $m$ with cluster-dependent observed cluster sizes | Unadjusted | $0.17$ | $5.33$ | $5.31$ | $0.95$ |  | $-0.18$ | $4.82$ | $4.48$ | $0.93$
GEE-g | $1.98$ | $3.75$ | $3.25$ | $0.86$ |  | $0.87$ | $4.53$ | $3.66$ | $0.88$
LMM-g | $1.74$ | $3.63$ | $3.69$ | $0.94$ |  | $0.79$ | $4.26$ | $3.83$ | $0.93$
Eff-PM | $0.08$ | $3.87$ | $3.47$ | $0.93$ |  | $0.09$ | $4.32$ | $3.97$ | $0.93$
Eff-ML | $-0.21$ | $3.42$ | $3.47$ | $0.95$ |  | -$0.06$ | $4.20$ | $3.84$ | $0.93$
Scenario 3: Large $m$ with random observed cluster sizes | Unadjusted | $0.01$ | $2.62$ | $2.65$ | $0.96$ |  | $0.06$ | $2.68$ | $2.63$ | $0.95$
GEE-g | $0.04$ | $1.40$ | $1.38$ | $0.95$ |  | $-0.01$ | $1.89$ | $1.83$ | $0.94$
LMM-g | $0.04$ | $1.40$ | $1.42$ | $0.96$ |  | $-0.01$ | $1.89$ | $1.55$ | $0.90$
Eff-PM | $0.00$ | $1.40$ | $1.38$ | $0.95$ |  | $-0.01$ | $1.93$ | $1.92$ | $0.95$
Eff-ML | $0.04$ | $0.70$ | $0.71$ | $0.95$ |  | $0.00$ | $0.77$ | $0.81$ | $0.97$
Scenario 4: Large $m$ with cluster-dependent observed cluster sizes | Unadjusted | $-0.04$ | $2.94$ | $2.89$ | $0.95$ |  | $-0.06$ | $2.57$ | $2.48$ | $0.94$
GEE-g | $1.80$ | $1.89$ | $1.82$ | $0.82$ |  | $0.74$ | $2.21$ | $2.12$ | $0.91$
LMM-g | $1.72$ | $1.87$ | $1.89$ | $0.85$ |  | $0.72$ | $2.20$ | $2.00$ | $0.91$
Eff-PM | $0.03$ | $1.89$ | $1.83$ | $0.94$ |  | $0.01$ | $2.21$ | $2.16$ | $0.94$
Eff-ML | $0.01$ | $1.91$ | $1.83$ | $0.94$ |  | $0.00$ | $2.23$ | $2.13$ | $0.94$

Unadjusted: the unadjusted estimator. GEE-g: generalized estimating equations
with weighted g-computation. LMM-g: linear mixed models with weighted
g-computation. Eff-PM: our proposed method with parametric working models.
Eff-ML: our proposed method with machine learning algorithms. ESE: empirical
standard error. ASE: average of estimated standard error. CP: coverage
probability based on $t$-distribution.

For scenarios 1 and 3, since the observed cluster size is completely random,
the model-based estimators perform well, which confirms the theoretical
results in Section 3. Among all methods, our method with machine learning
algorithms has the highest precision: its variance is $35\sim 93\%$ and
$49\sim 89\%$ smaller than the other methods for estimating the $\Delta_{C}$
and $\Delta_{I}$, respectively, demonstrating its potential to flexibly
leverage baseline covariates for improving study power. Our estimator that
uses parametric models for covariate adjustment has comparable precision to
model-based methods and is more efficient than the unadjusted estimator.

For scenarios 2 and 4, more individuals are enrolled in treated clusters with
a larger source population, leading to bias for methods that utilized
individual-level data without adjusting for this cluster-dependent sampling
scheme. Specifically, model-based methods have bias ranging from $0.72$ to
$1.98$, and $1\%$ to $13\%$ under coverage. In contrast to the issues caused
by adapting model-based methods, our proposed methods show both validity and
precision. The validity is reflected by the negligible bias and nominal
coverage, and the precision is borne out by their smaller empirical variance
than the unadjusted estimator.

Table 2 summarizes the simulation results for binary outcomes, and the
patterns are generally similar to those for continuous outcomes. In
particular, the proposed methods remain valid across scenarios, but the
advantage of machine-learning algorithms over parametric modeling is less
obvious. Finally, comparing a small versus large number of clusters, all
methods perform less stably with a smaller number of clusters, as expected.
When $m=30$, methods with covariate adjustment tend to underestimate the true
standard error, causing $0\sim 5\%$ under coverage. When $m$ increases to
$100$, the estimated standard errors match the empirical standard error,
thereby implying the validity of our variance estimator.

Table 2: Results in the second simulation experiment with binary outcomes.

|  |  Cluster-average treatment effect $\Delta_{C}=1.54$ |  |  Individual-average treatment effect $\Delta_{I}=1.18$
---|---|---|---|---
Setting | Method | Bias | ESE | ASE | CP |  | Bias | ESE | ASE | CP
Scenario 1: Small $m$ with random observed cluster sizes | Unadjusted | $0.04$ | $0.25$ | $0.24$ | $0.94$ |  | $0.03$ | $0.15$ | $0.13$ | $0.93$
GEE-g | $0.01$ | $0.23$ | $0.19$ | $0.92$ |  | $0.03$ | $0.22$ | $0.13$ | $0.93$
LMM-g | $0.01$ | $0.22$ | $0.19$ | $0.92$ |  | $0.03$ | $0.16$ | $0.13$ | $0.93$
Eff-PM | $0.01$ | $0.21$ | $0.20$ | $0.93$ |  | $0.02$ | $0.16$ | $0.15$ | $0.95$
Eff-ML | $0.03$ | $0.22$ | $0.20$ | $0.94$ |  | $0.03$ | $0.15$ | $0.14$ | $0.95$
Scenario 2: Small $m$ with cluster-dependent observed cluster sizes | Unadjusted | $0.05$ | $0.30$ | $0.29$ | $0.94$ |  | $0.03$ | $0.17$ | $0.14$ | $0.93$
GEE-g | $-0.18$ | $0.30$ | $0.20$ | $0.69$ |  | $-0.05$ | $0.16$ | $0.17$ | $0.89$
LMM-g | $-0.12$ | $0.24$ | $0.20$ | $0.78$ |  | $-0.04$ | $0.15$ | $0.16$ | $0.90$
Eff-PM | $0.04$ | $0.28$ | $0.26$ | $0.92$ |  | $0.03$ | $0.17$ | $0.17$ | $0.95$
Eff-ML | $0.06$ | $0.30$ | $0.30$ | $0.93$ |  | $0.05$ | $0.19$ | $0.17$ | $0.95$
Scenario 3: Large $m$ with random observed cluster sizes | Unadjusted | $0.01$ | $0.13$ | $0.13$ | $0.95$ |  | $0.01$ | $0.07$ | $0.07$ | $0.95$
GEE-g | $0.00$ | $0.11$ | $0.10$ | $0.95$ |  | $0.00$ | $0.07$ | $0.07$ | $0.95$
LMM-g | $0.00$ | $0.11$ | $0.10$ | $0.94$ |  | $0.01$ | $0.07$ | $0.07$ | $0.95$
Eff-PM | $0.01$ | $0.10$ | $0.10$ | $0.95$ |  | $0.00$ | $0.07$ | $0.08$ | $0.97$
Eff-ML | $0.01$ | $0.10$ | $0.10$ | $0.95$ |  | $0.01$ | $0.06$ | $0.07$ | $0.97$
Scenario 4: Large $m$ with cluster-dependent observed cluster sizes | Unadjusted | $0.01$ | $0.15$ | $0.15$ | $0.95$ |  | $0.01$ | $0.08$ | $0.08$ | $0.94$
GEE-g | $-0.20$ | $0.10$ | $0.10$ | $0.44$ |  | $-0.06$ | $0.07$ | $0.09$ | $0.88$
LMM-g | $-0.14$ | $0.11$ | $0.10$ | $0.67$ |  | $-0.05$ | $0.07$ | $0.08$ | $0.89$
Eff-PM | $0.01$ | $0.13$ | $0.13$ | $0.94$ |  | $0.00$ | $0.07$ | $0.08$ | $0.97$
Eff-ML | $0.01$ | $0.14$ | $0.13$ | $0.94$ |  | $0.01$ | $0.08$ | $0.08$ | $0.97$

Unadjusted: the unadjusted estimator. GEE-g: generalized estimating equations
with weighted g-computation. LMM-g: linear mixed models with weighted
g-computation. Eff-PM: our proposed method with parametric working models.
Eff-ML: our proposed method with machine learning algorithms. ESE: empirical
standard error. ASE: average of estimated standard error. CP: coverage
probability based on $t$-distribution.

In the Supplementary Material, we provide additional simulation results for
augmented generalized estimating equations and targeted maximum likelihood
estimation under the same settings. Augmented generalized estimating equations
were adapted with weighted g-computation to target our estimands and showed
similar performance to generalized estimating equations with weighted
g-computation. The target maximum likelihood estimator was unbiased when it
only adjusted for cluster-level covariates, while it had bias for $\Delta_{I}$
if the sampling was cluster-dependent. In most scenarios, both methods were
less precise than our proposed method coupled with machine learning estimators
for the nuisance functions.

##  6 Data applications

###  6.1 Three cluster-randomized experiments

The Work, Family, and Health Study (WFHS) is a cluster-randomized experiment
designed to reduce work-family conflict and improve the health and well-being
of employees (Work, Family, and Health Study , 2018, WFHS). Fifty-six study
groups (clusters) were equally randomized to receive a workplace intervention
(treatment) or not (control) with each cluster including 3–50 employees. The
observed cluster size has mean $11.77$ and standard error $7.47$. We focused
on the control over work hours outcome at the 6-month follow-up, a continuous
outcome measure ranging from 1 to 5. We adjusted for the following covariates:
cluster sizes and group job functions (core or supporting) at the cluster
level, and baseline value of control over work hours and mental health score
at the individual level.

The Pain Program for Active Coping and Training study (PPACT), supported by
National Institute of Health Pragmatic Clinical Trials Collaboratory, is a
pragmatic, cluster-randomized experiment evaluating the effectiveness of a
care–based cognitive behavioral therapy intervention for treating long-term
opioid users with chronic pain (DeBar et al., , 2022). One-hundred-six
clusters of primary care providers (clusters) were equally randomized to
receive the intervention or usual care. Each cluster contained 1-10
participants enrolled via phone screening. The observed cluster size has mean
$2.03$ and standard error $4.9$. The primary outcome was the PEGS (pain
intensity and interference with enjoyment of life, general activity, and
sleep) score at 12 months, a continuous scale assessing pain impact as a
composite of pain intensity and interference. We adjusted for the cluster
sizes and 12 individual-level baseline variables including age, gender,
disability, smoking status, body mass index, alcohol abuse, drug abuse,
comorbidity, depression, number of pain types, average morphine dose, and
heavy opioid usage.

The “Improving rational use of artemisinin combination therapies through
diagnosis-dependent subsidies” cluster-randomized trial, referred to as
“ACTS”, is a cluster-randomized experiment completed in western Kenya
(Prudhomme O’Meara et al., , 2018). Thirty-two community clusters were
randomized in a 1:1 ratio to two arms: malaria rapid diagnostic tests with
vouchers of artemisinin combination therapies provided for positive test
results (treatment) versus standard package (control). The primary outcome was
an indicator of receiving a malaria diagnostic test among fevers in the past
four weeks at 12 months. Using survey sampling, the observed cluster sizes
ranged from 39 to 129, with mean $56.62$ and standard error $16.09$. As the
source population size was unknown, we focused on inference of the observed
population in ACTS. We adjusted for cluster sizes and existence of health
facilities at the cluster level, and gender and wealth index at the individual
level.

The above three cluster-randomized experiments span social science, chronic
pain treatment, and infectious disease control with different cluster sizes,
outcome types, and number of covariates. By re-analyzing these data sets, we
aim to illustrate our methods in multiple real-world settings.

###  6.2 Results of data analysis

For each data set, we estimated the cluster-average and individual-average
treatment effects on the difference scale, using the unadjusted estimator,
generalized estimating equations, linear mixed models, and our proposed method
with parametric working models and machine learning algorithms. We first set
$M=N$ for each data analysis, corresponding to the source population analysis
for WFHS and the enrolled population analysis for PPACT and ACTS. For each
estimator, we reported the point estimates, 95% confidence interval based on
$t$-distribution, and proportion variance reduction compared to the unadjusted
estimator (PVR).

Table 3 summarizes the results of our data analyses. Across all studies and
both estimands, the unadjusted analysis has the widest confidence interval,
and covariate adjustment can bring variance reduction as high as 65%. While
weighted g-computation estimators have higher PVR than our proposed estimators
in WFHS and PPACT for estimating the cluster-average treatment effect, they
may be biased since the cluster size can be informative for the outcome,
thereby violating Assumption 4; in contrast, our proposed methods can remove
such bias and are hence more reliable. In the analysis of these three cluster-
randomized experiments, machine learning algorithms did not show an apparent
advantage over parametric working models.

Table 3: Results of full data analyses.

|  | Cluster-average treatment effect |  | Individual-average treatment effect
---|---|---|---|---
Study | Method | Estimate | 95% C.I. | PVR |  | Estimate | 95% C.I. | PVR
WFHS | Unadjusted | 0.17 | (0.00, 0.34) | - |  | 0.17 | (0.03, 0.31) | -
GEE-g | 0.21 | (0.11, 0.31) | 65% |  | 0.22 | (0.11, 0.32) | 44%
LMM-g | 0.21 | (0.10, 0.31) | 62% |  | 0.22 | (0.11, 0.32) | 44%
Eff-PM | 0.22 | (0.11, 0.33) | 58% |  | 0.21 | (0.12, 0.30) | 58%
Eff-ML | 0.22 | (0.10, 0.33) | 57% |  | 0.20 | (0.11, 0.30) | 55%
PPACT | Unadjusted | -0.83 | (-1.31, -0.35) | - |  | -0.70 | (-1.12, -0.28) | -
GEE-g | -0.53 | (-0.83, -0.22) | 60% |  | -0.53 | (-0.83, -0.22) | 47%
LMM-g | -0.53 | (-0.84, -0.21) | 56% |  | -0.52 | (-0.81, -0.22) | 52%
Eff-PM | -0.59 | (-0.98, -0.21) | 35% |  | -0.51 | (-0.81, -0.21) | 49%
Eff-ML | -0.83 | (-1.21, -0.44) | 35% |  | -0.71 | (-1.01, -0.40) | 48%
ACTS | Unadjusted | 0.08 | (-0.00, 0.16) | - |  | 0.08 | (-0.00, 0.16) | -
GEE-g | 0.07 | (-0.01, 0.15) | 9% |  | 0.07 | (-0.01, 0.15) | 9%
LMM-g | 0.07 | (-0.00, 0.15) | 12% |  | 0.07 | (-0.00, 0.15) | 12%
Eff-PM | 0.08 | (0.00, 0.15) | 20% |  | 0.07 | (-0.00, 0.15) | 18%
Eff-ML | 0.08 | (-0.00, 0.15) | 12% |  | 0.07 | (-0.00, 0.15) | 18%

Unadjusted: the unadjusted estimator. GEE-g: generalized estimating equations
with weighted g-computation. LMM-g: linear mixed models with weighted
g-computation. Eff-PM: our proposed method with parametric working models.
Eff-ML: our proposed method with machine learning algorithms. ESE: empirical
standard error. 95% C.I.: 95% confidence interval based on $t$-distribution.
PVR: proportional variance reduction compared to the unadjusted estimator.

To further illustrate our methods under the setting of cluster-dependent
sampling, we performed a simulation study based on the WFHS data. In each
simulated data,
$\\{(\widetilde{Y}_{i}(0),\widetilde{N}_{i},\widetilde{C}_{i},\widetilde{\mathrm{X}}_{i}),i=1,\dots,56\\}$
were sampled from the empirical distribution of observed data
$\\{(Y_{i},N_{i},C_{i},\widetilde{\mathrm{X}}_{i}):i=1,\dots,56\\}$. Then, we
randomly generated the treatment variable
$\widetilde{A}_{i}\sim\mathcal{B}(0.5)$ and
$\widetilde{Y}_{ij}(1)=\widetilde{Y}_{ij}(0)+N_{i}/3+\mathcal{N}(0,1)$. Next,
we set
$\widetilde{M}_{i}=\widetilde{A}_{i}\widetilde{C}_{i}(\widetilde{N}_{i}-2)+2$,
where $\widetilde{C}_{i}$ is the cluster-level covariate of group job
functions (core for $\widetilde{C}_{i}=1$ and supporting for
$\widetilde{C}_{i}=0$). Finally, for each cluster, we uniformly sampled
without replacement $\widetilde{M}_{i}$ individuals, for whom
$\widetilde{S}_{ij}=1$, and the observed data are
$\\{\widetilde{Y}_{ij},\widetilde{A}_{i},\widetilde{M}_{i},\widetilde{C}_{i1},\widetilde{\mathrm{X}}_{i}:\widetilde{S}_{ij}=1,i=1,\dots,m,j=1,\dots,\widetilde{N}_{i}\\}$
with
$\widetilde{Y}_{ij}=\widetilde{A}\widetilde{Y}_{ij}(1)+(1-\widetilde{A})\widetilde{Y}_{ij}(0)$.
We randomly generated 1,000 data sets and estimated both estimands as in
Section 5.

Table 4 gives the results for the cluster-dependent sampling analysis, which
mimics the real-world setting, i.e., the distribution of outcome and
covariates under control is based on real data and is unspecified. Consistent
with our simulation results in Section 5, the cluster-dependent sampling leads
to bias of model-based methods for both estimands, while our proposed methods
have negligible bias and nominal coverage. Furthermore, our proposed methods
implement flexible covariate adjustment via the working nuisance models and
clearly demonstrate higher precision than the unadjusted estimator.

Table 4: Results from the cluster-dependent sampling analyses based on $1,000$
simulation draws using the WFHS data set.

|  Cluster-average treatment effect ${\Delta}_{C}=3.92$ |  |  Individual-average treatment effect ${\Delta}_{I}=5.37$
---|---|---|---
Method | Bias | ESE | ASE | CP |  | Bias | ESE | ASE | CP
Unadjusted | 0.02 | 0.50 | 0.48 | 0.93 |  | 0.04 | 0.98 | 0.83 | 0.80
GEE-g | -0.71 | 0.40 | 0.40 | 0.56 |  | -1.28 | 0.74 | 0.55 | 0.39
LMM-g | -0.66 | 0.41 | 0.34 | 0.54 |  | -1.15 | 0.83 | 0.41 | 0.33
Eff-PM | 0.00 | 0.34 | 0.35 | 0.95 |  | 0.05 | 0.65 | 0.65 | 0.91
Eff-ML | 0.00 | 0.38 | 0.36 | 0.94 |  | 0.05 | 0.66 | 0.65 | 0.90

Unadjusted: the unadjusted estimator. GEE-g: generalized estimating equations
with weighted g-computation. LMM-g: linear mixed models with weighted
g-computation. Eff-PM: our proposed method with parametric working models.
Eff-ML: our proposed method with machine learning algorithms. ESE: empirical
standard error. ASE: average of estimated standard error. CP: coverage
probability based on $t$-distribution.

##  7 Concluding remarks

Our contributions to the emerging causal inference literature on cluster-
randomized experiments are two-fold. First, we clarified sufficient conditions
under which two mainstream model-based regression estimators are robust for
identifying the cluster-average treatment effect and individual-average
treatment effect, which are two important classes of estimands that can differ
in magnitude due to informative cluster size. Second, we derived the efficient
influence functions and proposed efficient estimators for the two classes of
estimands that allow for principled covariate adjustment and additionally
address the cluster-dependent sampling. The new estimators open the door for
flexibly leveraging parametric working models or machine learning algorithms
to model the potentially complex data generation mechanisms without
compromising the validity of inference in cluster-randomized experiments.

A common objective for covariate adjustment in cluster-randomized experiments
is to address chance imbalance and improve precision (Su and Ding, , 2021).
However, how to best select the optimal set and functional forms of covariates
is an open problem that remains to be addressed in future research. For
cluster-randomized experiments, this problem may be more challenging due to
the unknown intracluster correlations of the outcomes and covariates within
each cluster. In many cases where the investigators can only include a limited
number of clusters, there will be a trade-off between the loss of degrees of
freedom by adjusting for weakly prognostic covariates and the potential
asymptotic power gain by including more baseline variables. While the proposed
estimators may be a useful vehicle to incorporate variable selection
techniques in the working models, we maintain the recommendation of pre-
specifying prognostic covariates for adjustment in the design stage based on
subject-matter knowledge for practical applications.

Throughout the article, we have defined the cluster-average treatment effect
and individual-level treatment effect estimands as a function of the source
population size $N_{i}$, which can differ from the observed cluster size
$M_{i}$. Therefore, accurate identification of these estimands require
knowledge of the source population size. Conceptually, this source population
represents the set of eligible participants in each cluster that may be
recruited had the investigator obtained infinite financial and logistical
resources, and is precisely the set of individuals that the intervention is
designed to target. The availability of $N_{i}$ typically depends on the types
of clusters and the resource of the study. For example, the source population
size can be readily available if schools, worksites and villages are
randomized, whereas the source population size may be estimated from
historical data if clinics or hospitals are randomized. In this latter case,
one may set $N_{i}=M_{i}$ to perform an enrolled population analysis, which
implicitly assumes equivalence between the observed population and the source
population in each cluster. In the case where $N_{i}>M_{i}$ but no information
of $N_{i}$ is available, we stated in Remark 3 that only the cluster-average
treatment effect is identifiable. In future work, it would be interesting
establish alternative conditions to identify the individual-level treatment
effect in the absence of knowledge on source population size.

For developing the asymptotic properties of the model-robust estimators, we
have focused on simple randomization, whereas stratified randomization and
cluster rerandomization are useful alternative strategies to address baseline
imbalance and further improve the study power. For example, Wang et al.,
(2021) proved that stratified randomization leads to more efficient average
treatment effect estimators under a working linear mixed model in cluster-
randomized experiments. Lu et al.,  (2022) recently developed the asymptotic
theory for cluster rerandomization under a finite-population framework. It
would be useful to further extend our results to accommodate stratified
cluster randomization and cluster rerandomization.

## Acknowledgement

Research in this article was partially supported by a Patient-Centered
Outcomes Research Institute Award® (PCORI® Award ME-2020C3-21072). The
statements presented in this article are solely the responsibility of the
authors and do not necessarily represent the official views of PCORI®, or its
Board of Governors or Methodology Committee. The authors thank Drs. Gloria
Coronado and William Vollmer from the Center for Health Research, Kaiser
Permanente, for sharing the de-identified data from the STOP CRC study.

## Supplementary material

Supplementary material includes regularity conditions for Theorems 1-4,
proofs, consistent variance estimators, and review of other methods for
cluster-randomized experiments.

## References

  * Balzer et al., (2021) Balzer, L. B., van der Laan, M., Ayieko, J., Kamya, M., Chamie, G., Schwab, J., Havlir, D. V., and Petersen, M. L. (2021).  Two-Stage TMLE to reduce bias and improve efficiency in cluster randomized trials.  Biostatistics.  kxab043.
  * Balzer et al., (2019) Balzer, L. B., Zheng, W., van der Laan, M. J., and Petersen, M. L. (2019).  A new approach to hierarchical data analysis: targeted maximum likelihood estimation for the causal effect of a cluster-level exposure.  Statistical Methods in Medical Research, 28(6):1761–1780.
  * Benitez et al., (2021) Benitez, A., Petersen, M. L., van der Laan, M. J., Santos, N., Butrick, E., Walker, D., Ghosh, R., Otieno, P., Waiswa, P., and Balzer, L. B. (2021).  Comparative methods for the analysis of cluster randomized trials.  arXiv preprint arXiv:2110.09633.
  * Breslow and Clayton, (1993) Breslow, N. E. and Clayton, D. G. (1993).  Approximate inference in generalized linear mixed models.  Journal of the American statistical Association, 88(421):9–25.
  * Bugni et al., (2022) Bugni, F., Canay, I., Shaikh, A., and Tabord-Meehan, M. (2022).  Inference for cluster randomized experiments with non-ignorable cluster sizes.  arXiv preprint arXiv:2204.08356.
  * Chernozhukov et al., (2018) Chernozhukov, V., Chetverikov, D., Demirer, M., Duflo, E., Hansen, C., Newey, W., and Robins, J. (2018).  Double/debiased machine learning for treatment and structural parameters.  The Econometrics Journal, 21(1):C1–C68.
  * DeBar et al., (2022) DeBar, L., Mayhew, M., Benes, L., Bonifay, A., Deyo, R. A., Elder, C. R., Keefe, F. J., Leo, M. C., McMullen, C., Owen-Smith, A., et al. (2022).  A primary care–based cognitive behavioral therapy intervention for long-term opioid users with chronic pain: a randomized pragmatic trial.  Annals of Internal Medicine, 175(1):46–55.
  * Donner and Klar, (2000) Donner, A. and Klar, N. (2000).  Design and Analysis of Cluster Randomization Trials in Health Research.  London: Arnold.
  * Farrell et al., (2021) Farrell, M. H., Liang, T., and Misra, S. (2021).  Deep neural networks for estimation and inference.  Econometrica, 89(1):181–213.
  * Hines et al., (2022) Hines, O., Dukes, O., Diaz-Ordaz, K., and Vansteelandt, S. (2022).  Demystifying statistical learning based on efficient influence functions.  The American Statistician, pages 1–13.
  * Imai et al., (2009) Imai, K., King, G., and Nall, C. (2009).  The essential role of pair matching in cluster-randomized experiments, with application to the Mexican universal health insurance evaluation.  Statistical Science, 24(1):29–53.
  * Kahan et al., (2022) Kahan, B. C., Li, F., Copas, A. J., and Harhay, M. O. (2022).  Estimands in cluster-randomized trials: choosing analyses that answer the right question.  International Journal of Epidemiology.
  * Li et al., (2022) Li, F., Tian, Z., Bobb, J., Papadogeorgou, G., and Li, F. (2022).  Clarifying selection bias in cluster randomized trials.  Clinical Trials, 19(1):33–41.
  * Liang and Zeger, (1986) Liang, K.-Y. and Zeger, S. L. (1986).  Longitudinal data analysis using generalized linear models.  Biometrika, 73(1):13–22.
  * Lu et al., (2022) Lu, X., Liu, T., Liu, H., and Ding, P. (2022).  Design-based theory for cluster rerandomization.  Biometrika.
  * Luo and Spindler, (2016) Luo, Y. and Spindler, M. (2016).  High-dimensional $l_{2}$boosting: Rate of convergence.  arXiv preprint arXiv:1602.08927.
  * Middleton and Aronow, (2015) Middleton, J. A. and Aronow, P. M. (2015).  Unbiased estimation of the average treatment effect in cluster-randomized experiments.  Statistics, Politics and Policy, 6(1-2):39–75.
  * Murray et al., (1998) Murray, D. M. et al. (1998).  Design and Analysis of Group-Randomized Trials, volume 29.  Oxford University Press, USA.
  * Neyman et al., (1990) Neyman, J. S., Dabrowska, D. M., and Speed, T. (1990).  On the application of probability theory to agricultural experiments. Essay on principles. Section 9.  Statistical Science, pages 465–472.
  * Prudhomme O’Meara et al., (2018) Prudhomme O’Meara, W., Menya, D., Laktabai, J., Platt, A., Saran, I., Maffioli, E., Kipkoech, J., Mohanan, M., and Turner, E. L. (2018).  Improving rational use of acts through diagnosis-dependent subsidies: Evidence from a cluster-randomized controlled trial in western kenya.  PLoS medicine, 15(7):e1002607.
  * Schochet et al., (2021) Schochet, P. Z., Pashley, N. E., Miratrix, L. W., and Kautz, T. (2021).  Design-based ratio estimators and central limit theorems for clustered, blocked RCTs.  Journal of the American Statistical Association, pages 1–22.
  * Stephens et al., (2012) Stephens, A. J., Tchetgen Tchetgen, E. J., and Gruttola, V. D. (2012).  Augmented generalized estimating equations for improving efficiency and validity of estimation in cluster randomized trials by leveraging cluster-level and individual-level covariates.  Statistics in Medicine, 31(10):915–930.
  * Su and Ding, (2021) Su, F. and Ding, P. (2021).  Model-assisted analyses of cluster-randomized experiments.  Journal of the Royal Statistical Society, Series B, 83(5):994–1015.
  * Turner et al., (2017) Turner, E. L., Prague, M., Gallis, J. A., Li, F., and Murray, D. M. (2017).  Review of recent methodological developments in group-randomized trials: Part 2-analysis.  American Journal of Public Health, 107(7):1078–1086.
  * van der Laan et al., (2007) van der Laan, M. J., Polley, E. C., and Hubbard, A. E. (2007).  Super learner.  Statistical applications in genetics and molecular biology, 6(1).
  * van der Laan et al., (2011) van der Laan, M. J., Rose, S., et al. (2011).  Targeted Learning: Causal Inference for Observational and Experimental Data, volume 10.  Springer.
  * van der Vaart, (1998) van der Vaart, A. (1998).  Asymptotic Statistics.  Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.
  * Wager and Walther, (2015) Wager, S. and Walther, G. (2015).  Adaptive concentration of regression trees, with application to random forests.  arXiv preprint arXiv:1503.06388.
  * Wang et al., (2021) Wang, B., Harhay, M. O., Small, D. S., Morris, T. P., and Li, F. (2021).  On the mixed-model analysis of covariance in cluster-randomized trials.  arXiv preprint arXiv:2112.00832.
  * Wang et al., (2022) Wang, X., Turner, E. L., Li, F., Wang, R., Moyer, J., Cook, A. J., Murray, D. M., and Heagerty, P. J. (2022).  Two weights make a wrong: Cluster randomized trials with variable cluster sizes and heterogeneous treatment effects.  Contemporary Clinical Trials, 114:106702.
  * Work, Family, and Health Study , 2018 (WFHS) Work, Family, and Health Study (WFHS) (2018).  Work, family and health network.  Inter-university Consortium for Political and Social Research [distributor].