Improving the Estimation of Site-Specific Effects and their Distribution in Multisite Trials

Suppose a multisite trial consists of $J$ sites, indexed by $j=1,\ldots,J$, where identical treatments are administered within each site. The analytic setting considered in this paper is based on the Rubin (1981) model for parallel randomized experiments, also known as a random-effects meta-analysis (Raudenbush and Bryk,, 1985). For data inputs, this model requires only the observed or estimated effects $\widehat{\tau}_{j}$ from each of the $J$ sites, along with their corresponding squared standard errors $\widehat{se}_{j}^{2}$. The $\widehat{\tau}_{j}$ are estimating $\tau_{j}$, the true average impact in site $j$; under the potential outcomes framework, this would be the difference in the average outcomes of units if all were treated minus if all were control; we would make the usual assumptions about non-interference of treatment and a well-defined treatment (SUTVA) both within and across sites to make the estimand well defined (see, e.g., Imbens and Rubin, (2015) for further discussion). The $\widehat{\tau}_{j}$ and $\widehat{se}_{j}^{2}$ values are obtained through maximum likelihood (ML) estimation using data exclusively from site $j$. Our meta-analytic approach closely aligns with the fixed intercept, random coefficient (FIRC) estimator (Raudenbush and Bloom,, 2015), as it focuses on modeling only the true, randomly varying treatment effect $\tau_{j}$. This approach effectively circumvents some of the strong joint distributional assumptions typically made between $\tau_{j}$ and a randomly varying intercept or control group mean (Bloom et al.,, 2017).

The first stage of the hierarchical model describes the relationship between the observed statistics $\widehat{\tau}_{j}$ and the parameter $\tau_{j}$:

The second stage assumes that the site average impacts $\tau_{j}$ are independent and identically distributed ($i.i.d.$) with a specific prior distribution²²2We refer to $G$ as a ‘prior distribution’ to be consistent with standard practices in Bayesian hierarchical models (Gelman et al.,, 2013; Rabe-Hesketh and Skrondal,, 2022). In this framework, $G$ represents our prior knowledge or beliefs about the parameter $\tau_{j}$ before observing the data from site $j$. In contrast, within the frequentist framework, $G$ is interpreted as a distributional assumption regarding $\tau_{j}$. $G$,

The Bayesian hierarchical framework plays a vital role in separating the genuine heterogeneity in true effects $\tau_{j}$ across different sites ($\sigma^{2}$) from the sampling variation of each estimated effect $\widehat{\tau}_{j}$ around its true effect within sites ($\widehat{se}_{j}^{2}$). A key focus in our analytic setting is to assess the relative magnitudes of $\sigma^{2}$ and $\widehat{se}_{j}^{2}$, aiming to understand how this signal-to-noise ratio influences the efficacy of various estimation strategies. Particularly when $\widehat{se}_{j}^{2}$ is relatively large, the $\widehat{\tau}_{j}$ values will generally be overdispersed compared to true $\tau_{j}$ values. This can lead to a misrepresentation of the rank order of effects across different sites; e.g., when the $\widehat{se}_{j}^{2}$ vary due to factors such as site size, the smaller sample sizes might produce more extreme estimates. When estimating conditional effects, incorporating unit-level covariates into our analysis could concurrently reduce both $\sigma^{2}$ and $\widehat{se}_{j}^{2}$, although these reductions might occur at different rates, or in some cases, $\sigma^{2}$ might even increase. This introduces an additional layer of complexity to the comparative analysis of these two variances. Therefore, to ensure focused and clear exposition and formulas regarding the comparisons of $\sigma^{2}$ and $\widehat{se}j^{2}$, our study deliberately defers the inclusion of conditional effects and covariates interacted with treatment to future research; covariates for precision gain will simply decrease the standard errors and increase the information of the estimates regarding the estimands and all our results carry.

In the special case of the normality assumption for $\tau_{j}$ the conditional posterior distribution of $\tau_{j}$ given the hyperparameters $\tau$ and $\sigma^{2}$ is normal and the corresponding conditional posterior mean and variance of $\tau_{j}$ have simple closed-form expressions:

To assess the magnitude of $\sigma^{2}$ relative to the sampling error $\widehat{se}_{j}^{2}$ for each site, we use the shrinkage factor $S_{j}$ as a direct measure for site $j$. For a comprehensive evaluation across all $J$ sites, we calculate the geometric mean of $\widehat{se}_{j}^{2}$, yielding an average measure of reliability for the ML estimator $\hat{\tau}_{j}$, which we call the informativeness ($I$):

The Rubin model highlights two estimands of interest: the super population distribution of site-specific effects, $G$, and the finite-population distribution, $\left\{\tau_{j}\right\}_{j=1}^{J}$, for the $J$ sites in the experiment. These two estimands correspond to distinct research objectives. Super population inference aims to recover $G$, the effect distribution for the broader population of sites from which the $J$ experimental sites were sampled. This task is referred to as a Bayes deconvolution problem, wherein the observed, noisy sample estimates $\left\{\hat{\tau}_{j},\widehat{se}_{j}^{2}\right\}_{j=1}^{J}$ are used to compute an estimate $\widehat{G}$ of $G$ (Efron,, 2016). Performance evaluation studies typically focus on this research objective, as they often have sufficiently large $J$ (e.g., Gilraine et al.,, 2020). In contrast, finite-population inference examines the distribution of the true site-effect parameters for the specified $J$ sites, taking these sites as the entire finite population of interest.

In multisite trial settings, recovering the super population estimand $G$ can be challenging without assumptions about the form of $G$. For a normal model of $G$, only two hyperparameters, $\tau$ and $\sigma^{2}$, require estimation in the deconvolution process. However, when the prior $G$ is unknown, the performance of deconvolution estimators heavily relies on the shape of the underlying $G$ and the number of sites $J$ (McCulloch and Neuhaus,, 2011). Specifically, the mathematical guarantees of flexible estimators are typically asymptotic in nature, necessitating a large $J$ (Jiang and Zhang,, 2009). In the context of multisite trials in education, however, $J$ is often small to moderate (Weiss et al.,, 2017), rendering such asymptotics potentially less applicable. The core concern is that a random sample of 30 to 40 sites’ $\tau_{j}$, even if obtained without measurement error, may not give a sharp picture of the target distribution. Furthermore, it may be unreasonable to assume that the experimental sites in many multisite trials are adequately representative of the super population (Stuart et al.,, 2011). For example, Tipton et al., (2021) found that schools involved in grants from the Institute of Education Science were predominantly located in large urban districts, geographically close to one another and to the study’s principal investigators.

In this paper, we focus on finite-population estimands, represented by $\left\{\tau_{j}\right\}_{j=1}^{J}$. This approach helps us circumvent the risks associated with generalizing treatment effect estimates from a sample of sites to a super population when systematic discrepancies exist between them. Moreover, concentrating on a finite-population of sites typically results in more precise effect estimates compared to situations where researchers aim to generalize to a super population (Miratrix et al.,, 2021). This approach is statistically less demanding, allowing for more relevant recommendations in low-data settings for multisite trials.

Shen and Louis, (1998) identify three common inferential goals regarding a set of site-specific parameters: (1) estimating the individual site-specific effect parameters, $\tau_{j}$, (2) ranking the sites based on $\tau_{j}$, and (3) estimating the empirical distribution function (EDF) of the $\tau_{j}$’s. Each of these three goals can be associated with a different loss function, which may be minimized by a distinct estimator. The loss function typically associated with the first goal is the mean-squared error loss (MSEL):

For the second inferential goal, we can choose an estimator for the vector of ranks of $\tau_{j}$ that minimizes the mean squared error loss of the ranks (MSELR),

For the third inferential goal, Shen and Louis, (1998) suggest minimizing the integrated squared-error loss (ISEL):

The three inferential goals, their associated loss functions, and their optimal estimators reveal two primary challenges in achieving valid finite-population estimation. The first challenge is model misspecification, which arises when we assume an incorrect parametric form for the super population distribution $G$. If the true distribution $G$ is not Gaussian, assuming normality may result in insensitivity to skewness, long tails, multimodality, and other complexities in the EDF of the $\tau_{j}$’s (McCulloch and Neuhaus,, 2011). Consequently, estimators based on the standard Rubin model become unreliable, especially for the third inferential goal. The second challenge emerges when, for a given goal, an unsuitable estimator is chosen, even when the prior distribution $G$ is accurately specified in a model. The optimal estimator is contingent upon the chosen loss function or inferential goal. However, practitioners often use the same set of posterior mean effect estimates for all three goals, leading to suboptimal outcomes for at least some of them. For instance, the EDF of posterior mean effect estimates will tend to be underdispersed compared to the EDF of $\tau_{j}$ due to shrinkage towards the prior mean effect $\tau$. Conversely, the EDF of raw observed ML effect estimates $\hat{\tau}_{j}$ is overdispersed because of sampling error (Mislevy et al.,, 1992).

To address the threats outlined in the previous section, two strategies can be employed. The first strategy entails adopting flexible semiparametric or nonparametric specifications for the prior distribution $G$ to protect against model misspecification (Paddock et al.,, 2006). One prominent Bayesian nonparametric specification is the Dirichlet process (DP) prior, which has gained widespread use in the literature (Paganin et al.,, 2022). The second strategy focuses on employing posterior summary methods, such as the constrained Bayes or the triple-goal estimators. These estimators are designed to directly target the loss functions associated with specific inferential goals. The following sections discuss these two strategies and their implications for addressing the challenges in finite-population estimation.

Our primary target of interest is $\left\{\tau_{j}\right\}_{j=1}^{J}$, representing the EDF of the $\tau_{j}$’s. Both CB and GR generate estimates of the $\tau_{j}$ that are not as overly shrunk as the simple posterior means. The CB estimator directly modifies the posterior mean estimates $\tau_{j}^{*}$ to avoid under-dispersion, rescaling them to ensure that the estimated variance of site-specific effect estimates aligns with the estimated marginal variance of the latent parameters $\tau_{j}$. In equation (3), the posterior mean of $\tau_{j}$ is denoted as $\tau_{j}^{*}\equiv\mathrm{E}\left(\tau_{j}\mid\hat{\tau}_{j}\right)$, and the posterior variance of $\tau_{j}$ as $V_{j}\equiv\operatorname{Var}\left(\tau_{j}\mid\hat{\tau}_{j}\right)$. We introduce $\bar{\tau}\equiv J^{-1}\sum\tau_{j}^{*}$ to represent the finite-population mean of the posterior mean estimates and $v$ as the corresponding finite-population variance. The CB estimates are then:

In contrast, the GR estimator employs evenly spaced percentiles of the estimated $G$ to generate its estimates. The GR estimator is constructed via a three-stage process. In the first stage, an estimate of the distribution function that minimizes ISEL is constructed as

Instead of presupposing a known parametric form for the prior distribution $G$ in equation (2), a Dirichlet process (DP) can be employed as the prior for $G$ to accommodate uncertainty regarding its shape. Conceptually, the DP model can be viewed as an infinite mixture model in which both the priors and the data contribute to determining the number of observed mixture components.

The DP prior is characterized by two hyperparameters: a base distribution $G_{0}$ and a precision parameter $\alpha$ (Antoniak,, 1974). A two-stage hierarchical model, incorporating the DP prior, can be formulated as follows:

To finalize the model specification, we must define the base distribution $G_{0}$ and establish a hyperprior for the parameter $\alpha$. $G_{0}$ serves as the initial best guess for the shape of the prior distribution $G$. Typically, we assume $G_{0}$ follows a Gaussian distribution with mean $\tau$ and variance $\sigma^{2}$. For computational convenience, the base distribution $G_{0}$ is typically constructed as the product of conjugate distributions, such as a normal distribution for the mean parameter $\tau$ and an inverse-gamma distribution for the variance parameter $\sigma^{2}$ (Paganin et al.,, 2022).

The precision hyperparameter $\alpha$ controls the degree to which $G$ converges toward $G_{0}$ (see the Supplemental Material). The larger the $\alpha$, the larger the number of unique observed site impact estimates $K$. It is worth noting that $K$ does not always precisely correspond to the number of mixture components, $C$, representing latent subpopulations with substantive meaning, as would be defined in a finite mixture modeling approach. However, $K$ can be reasonably regarded as an upper limit for $C$ (Ishwaran and Zarepour,, 2000).

The hyperprior for $\alpha$ plays a pivotal role in determining $K$, which in turn shapes the posterior distribution over clusters. Using a Gamma distribution - with $a$ as the shape parameter and $b$ as the rate parameter - as a hyperprior for $\alpha$ is common practice, attributable to its computational advantages (Escobar and West,, 1995). While the choices of parameters $a$ and $b$ can significantly influence the posterior distribution of $\alpha$, thus impacting clustering behavior (Dorazio,, 2009; Murugiah and Sweeting,, 2012), other studies suggest that data are typically sufficient to produce a concentrated posterior, thereby reducing the influence of the hyperprior even when a high-variance prior for $\alpha$ is used (Leslie et al.,, 2007; Gelman et al.,, 2013).

Our objective is to assess the sensitivity under two distinct hyperprior options for $\alpha$, particularly in the context of varying levels of data informativeness. The first hyperprior, explored by Antonelli et al., (2016), is relatively diffuse, allowing for numerous clusters. This strategy, which we refer to as the DP-diffuse model, is applicable when prior knowledge about $K$ is absent. The DP-diffuse model chooses values of $a$ and $b$ such that $\alpha$ is centered between 1 and $J$, with a large variance to assign a priori mass to a wide range of possible $\alpha$ values. For instance, with $J=50$, we might set the mean of the $\alpha$ distribution to 25 and its variance to 250 , ten times the mean. With these initial values, we can use the moments of a Gamma distribution to derive the corresponding hyperparameters $a=2.5$ and $b=0.1$, based on the formulas $\mathrm{E}(\alpha\mid a,b)=$ $a/b$ and $\operatorname{Var}(\alpha\mid a,b)=a/b^{2}$.

This research proposes a novel approach, employing a $\chi^{2}$ distribution to formulate an informative hyperprior for $\alpha$. Assume we have $\operatorname{Pr}(K)$, a prior that encapsulates our knowledge of the distribution of $K$. Dorazio, (2009) provides an expression for $\operatorname{Pr}(K\mid J,a,b)$, the probability mass function for the prior distribution $K$ induced by a Gamma $(a,b)$ prior for $\alpha$ and $J$. We can then obtain optimal values for $a$ and $b$ by minimizing the discrepancy between the encoded prior $\operatorname{Pr}(K)$ and the $\alpha$-induced prior $\operatorname{Pr}(K\mid J,a,b)$ using the Kullback-Leibler divergence measure.

We propose setting $\operatorname{Pr}(K)\sim\chi^{2}(\mathrm{df}=u)$ to more intuitively encapsulate our prior knowledge on $K$ and its inherent uncertainty. This strategy is primarily motivated by the fact that the $\chi_{u}^{2}$ distribution has a mean and variance of $u$ and $2u$, respectively. For instance, if we regard five clusters ($K=5$) among 50 sites as a reasonable approximation of the underlying distribution, we could assume that $\operatorname{Pr}(K)$ follows a $\chi^{2}(5)$ distribution and specify a $\operatorname{Gamma}(a,b)$ that closely matches $\chi^{2}(5)$. In this case, the appropriate Gamma distribution would have parameters $(a,b)=(1.60,1.22)$ (see the Supplemental Material). This approach, which we term the DP-inform model, is particularly effective for crafting an informative prior for $K$, especially when aiming to impose near-zero probabilities beyond a certain threshold (e.g., $K=25$ as shown in Figure 1) and to clarify the prior mean and variance of $K$.

In Figure 1, we observe that the two DPM models represent contrasting beliefs regarding $K$. The DP-diffuse model favors a smoother $G$ with approximately 25-30 clusters for 50 sites, while the DP-inform model, based on $\chi^{2}(5)$, is more favorable to a more uneven $G$, assuming that about five clusters are needed for 50 sites. Further details about the choices of hyperpriors and their implementations can be found in the Supplemental Material.

To compare the performance of different models and posterior summary methods, we conduct a comprehensive Monte Carlo simulation. We generate data by systematically varying five factors: (a) the number of sites ($J$), (b) the average number of observations per site ($\bar{n}_{j}$), (c) the coefficient of variation in site sizes ($\mathrm{CV}$), (d) the true cross-site population distribution $(G)$ of $\tau_{j}$, and (e) the true cross-site population standard deviation of impacts ($\sigma$). We calibrate our simulated datasets to mirror typical data settings for multisite trials in education, as described by Weiss et al., (2017).

We consider five values for the number of sites in our study: $J=25,50,75,100$, and $300$. The number of random assignment blocks reported in the 16 studies from Weiss et al., (2017) ranged from 20 to 356, with a median of 78. Although multisite trials typically involve a small to moderate numbers of sites, we also wanted to investigate model performance for larger multisite trials, such as the Head Start Impact Study (e.g., Bloom and Weiland,, 2015).

We vary the average number of observations per site in our study as $\bar{n}_{j}=10,20,40,80$, and $160$. Weiss et al., (2017) reported a range of mean site sizes, from 11 in the Head Start Impact Study to 1,176 in the Welfare-to-Work Program, with the 25th and 75th percentiles at 75 and 163, respectively. Larger site sizes also capture contexts where we are using baseline covariates to increase precision, as both size and $R^{2}$ directly impact site-level uncertainty, which is what drives our overall performance (see equation (1) and the definition of informativeness in equation (4)). We control the variation in site sizes by setting the coefficient of variation of the site sizes to be $\mathrm{CV}=0.00,0.25,0.50,0.75$ by sampling $n_{j}$ from a gamma distribution with a mean of $\bar{n}_{j}$ and a standard deviation of $\bar{n}_{j}\cdot\mathrm{CV}$. We then truncate $n_{j}$ at the lower limit of $n_{j}=5$. Finally, we compute the within-site sampling variances as: $\widehat{se}_{j}^{2}=4/n_{j}$ (this is assuming we are in an effect size metric with a control-side standard deviation of 1 within site).

We vary the cross-site impact standard deviation as $\sigma=0.05,0.10,0.15,0.20,$ and $0.25$ in effect size units. Weiss et al., (2017) report that the cross-site standard deviation of program impact, in effect size units, ranges between 0.00 and 0.35, with most values falling between 0.10 and 0.25. They define 0.05, 0.15, and 0.25 as having a modest, moderate, and substantial impact variation, respectively. In conjunction with the range of $\widehat{se}_{j}^{2}$ determined by our selection of $\bar{n}_{j}$, $\sigma$ determines the magnitude of the average reliability of the ML estimates, denoted by $I$. Across our simulation conditions, $I$ varies between 0.01 and 0.71, with a mean of 0.25.

Finally, we consider three different population distributions $G$ for $\tau_{j}$: Gaussian, a mixture of two Gaussians, and asymmetric Laplace (AL). The Gaussian mixture represents a scenario where the actual $G$ is multimodal, while the AL distribution exemplifies a situation where the true $G$ is skewed and exhibits a long tail. We rescale the simulated $\tau_{j}$ so that the variance of $G$ equals $\sigma$. Additional details on the data-generation process can be found in the Supplemental Material.

Using the five design factors ($J$, $\bar{n}_{j}$, CV, $\sigma$, and $G$), we generate 1,500 simulation conditions, each with 100 replications. We fit each replication with three data-analytic models: Gaussian, DP-diffuse, and DP-inform, and, for each model, use Markov Chain Monte Carlo to obtain 4,000 posterior samples from the joint posterior distributions of the J site-specific effects. Finally, we apply all three posterior summary methods (PM, CB, and GR) to each model fit. This gives 9 estimators for each generated dataset.

To analyze the simulation results, we separate the three true population distributions $G$ into samples of 450,000 observations each, generated as a factorial combination of the four remaining design factors ($J$, $\bar{n}_{j}$, CV, and $\sigma$), two data-analytic factors (data-analytic model and posterior summary method), and 100 replications. We construct 50,000 cluster (i.e., dataset) identifiers based on the four design factors and 100 replications to account for the statistical dependence induced by analyzing the same simulated dataset.

We use meta-model regressions (Skrondal,, 2000) to examine the relationship between performance criteria and simulation factors. Unlike visual or descriptive methods, meta-models allow for a more precise detection of significant patterns while accounting for uncertainty from Monte Carlo error (Boomsma,, 2013). We focus on the three performance criteria outlined in Section 3.1 as our target outcome variables: RMSE, ISEL, and MSELP. To ease the interpretation of the meta-model regression coefficients, we apply a logarithmic transformation to the outcome variables. As for the explanatory variables, we create a set of dummy variables representing the six design and data-analytic factors, along with their two-way interaction terms. The reference groups are $J=25$, $\bar{n}_{j}=10$, $\sigma=0.05$, $\text{CV}=0.00$, the Gaussian model, and the PM summary method. Lastly, since we include each dataset nine times in each meta-model regression (corresponding to every model/posterior summary method combination), we cluster our standard errors on datasets.

We now examine a series of case studies to assess the relative performance of different model-method combinations in various real-world situations. Figure 5 displays the performance of each combination in a simulation scenario characterized by high informativeness, encompassing considerable between-site information ($J=300$, $\sigma=0.25$) and significant within-site information ($\bar{n}_{j}=80$). We note that Weiss et al., (2017) do not provide directly comparable cases, highlighting the challenge of acquiring highly informative data at both the within- and between-site levels in real-world multisite trials.

Panel A of Figure 5 plots histograms of the estimates made by the different model-method combinations against the super population density of true site-specific effects. Panel A aggregates estimates across all simulations, so it does not show how well these approaches estimate the true finite-population distribution of site-specific effects for any given simulation sample; nonetheless, it illustrates the general behaviors of each of these approaches. We see that in this highly informative setting, the DP models can consistently adapt to the sampled sites, while the Gaussian estimator tends to shrink estimates toward the overall average. Similarly, even in this highly informative setting, the PM estimator tends to shrink estimates much more than the CB or GR estimators.

To summarize these results, we visualize the log-outcomes predicted by the meta-model regressions in Panels B and C of Figure 5. In this highly informative setting, the DP models outperform the Gaussian model in terms of both RMSE and ISEL. For example, using the DP-inform model with the GR estimator significantly improves ISEL compared to the Gaussian model with the same estimator. We see no significant differences in performance between the two DP models, suggesting that with sufficient information, the DP models are quite data-adaptive and insensitive to the choice of prior on $\alpha$.

Figure 6 shows results in a scenario similar to the Head Start Impact Study (Puma et al.,, 2010), which has significant between-site information ($J=317$, $\sigma=0.30$) but limited within-site information ($\bar{n}_{j}=11.3$). In Panel A, we see how the lack of within-site information leads even the DP models to lean more heavily on the Gaussian base distribution, leading to roughly Gaussian histograms when aggregated across all simulations.

When there is significantly more between-site information than within-site information, shrinkage toward the overall mean effect can considerably influence model estimates. For instance, the DP-diffuse model is found to be unsuitable for both RMSE and ISEL due to undershrinkage, which can be attributed to the model’s assumption of a large number of clusters and the substantial between-site information. Consequently, this leads to excessive weight being placed on a broad spectrum of potential site-specific effects. Conversely, the Gaussian model exhibits excessive shrinkage as a result of its shape constraints. In this scenario, the DP-inform model demonstrates strong performance in terms of ISEL, as its flexible shape allows it to adapt appropriately to the considerable between-site information, while its limited number of assumed clusters prevents undershrinkage.

In general, the choice of posterior summary method has a more pronounced impact than model selection. Using the PM summary, RMSE improves by 16-42% compared to CB or GR methods. In contrast, only a 1% RMSE difference is observed between the Gaussian and DP-inform models with PM. For ISEL, CB or GR methods outperform PM by 1.85 to 2.14 times, a more significant gain than the 20-35% improvement when choosing the DP-inform over the Gaussian model with CB/GR methods.

Figure 7 presents simulation results for a setting akin to the Performance-Based Scholarships study, characterized by limited between-site information ($J=39$, $\sigma=0.05$) but abundant within-site information ($\bar{n}_{j}=177.9$). In this scenario, the choice of model does not substantially influence performance once a suitable posterior summary method is employed (i.e., PM for RMSE, CB/GR for ISEL). The data-adaptive nature of the DP models is evident; they exhibit less shrinkage than the Gaussian model in Case 2, where $\sigma$ is large,