Transfer Estimates for Causal Effects across Heterogeneous Sites

Transfer Estimates for Causal Effects across Heterogeneous Sites

\newtheorem{dfn}{Definition}[section] \newtheorem{rem}{Remark}[section] \newtheorem{cor}{Corollary}[section] \newtheorem{thm}{Theorem}[section] \newtheorem{lem}{Lemma}[section] \newtheorem{notn}{Notation}[section] \newtheorem{con}{Condition}[section] \newtheorem{prp}{Proposition}[section] \newtheorem{pty}{Property}[section] \newtheorem{ass}{Assumption}[section] \newtheorem{ex}{Example}[section] \newtheorem*{cst1}{Constraint S} \newtheorem*{cst2}{Constraint U} \newtheorem{qn}{Question}[section]

\onehalfspacing

Introduction

When scaling up an intervention or planning an implementation at a new location, it is often necessary to extrapolate experimental evidence to new sites or contexts. In such settings, average causal effects typically vary across contexts due to environmental factors, only some of which may be observed. We consider a problem in which cross-sectional information on outcomes and covariates is available for both experimental and target sites, and we formalize a process of predicting a causal response that uses disaggregated pre-intervention (baseline) outcome data from the target location to predict such a model shift reflecting site-specific heterogeneity. The underlying premise of such an approach is that the data-generating processes for potential outcomes for pre- and post-intervention outcomes are likely similar, and depend on the same unit- and site specific factors, so that baseline outcomes are in fact predictive for treatment effects. Such an assumption may be particularly plausible when the effect of the intervention is expected to be only incremental rather than fundamentally altering the relationship between unit or site characteristics and the outcome of interest.

It has widely been recognized that pre-intervention outcomes can be useful to predict or control for unobserved heterogeneity at the level of the individual unit (see e.g. DWa99). This paper proposes a strategy for doing so at the level of the entire subpopulation to account for shared unobserved heterogeneity at the level of the site rather than the individual. To that end, the relevant baseline information for a given site consists of the full conditional distribution of pre-intervention outcomes given unit covariates, that is we view the baseline as functional data. This choice is motivated by the observation that unobserved site-specific confounders may generally manifest themselves not only in average levels of outcomes, but also how these interact with observed unit-specific attributes. However, in most practically relevant settings, the number of observed sites is not large, forcing the researcher to make pragmatic decisions on how flexibly to model the observable data.\footnote{All15 considers a setting in which a policy was initially evaluated at 10 sites and eventually scaled up to 111 separate sites. DPS21 use 142 year/country samples from 61 different countries. The PROGRESA study of conditional cash transfers in Mexico was initially conducted in 506 rural communities across 7 states in Mexico (see TWo06). Mea22 aggregates across seven different RCTs for micro-credit interventions published in 2015.} The corresponding problem of predicting conditional average treatment effects from baseline outcome data can be viewed as functional regression where a realistic implementation can at best achieve a highly regularized solution. Moreover these data constraints also make it all the more important to choose a procedure that makes statistically efficient use of the available data.

Our approach corresponds to a finite-dimensional approximation to that problem, where we determine the optimal feature space in which to solve a linear version of the prediction problem. In our leading application, cross-validation recommends the use of as few as $K=2$ features for prediction, resulting in a highly regularized estimator. Compared to Ridge and other alternative regularization schemes, the resulting transfer estimate can always be interpreted as the best linear predictor given those constructed site-specific features regardless of the degree of regularization. We can furthermore assess whether there exist sites in the experimental population that are similar to a target location in terms of these site characteristics that were determined to be most predictive of conditional average treatment effects. Similar techniques could in principle be developed to predict conditional treatment effects for sites within the experimental sample when treatment assignment was randomized at the site level.

Conditioning on baseline data presumes a statistical framework that defines a joint distribution for pre- and post-intervention outcomes across sites. We choose a fixed-population construction that regards the combined (finite) population of experimental and target sites as fixed, but assumes that the number of cross-sectional units within each cluster is large. Statistical properties of extrapolation estimators are then evaluated under a randomization protocol that assigns experimental versus target status at random among those clusters. In analogy with the literature on conformal prediction, the constructed statistical experiment treats experimental and target locations as finitely exchangeable. We do not necessarily regard this assignment mechanism as factually accurate - e.g. the observed assignment may likely exhibit site-selection effects of the kind documented by All15. Rather, this data generating process can alternatively be viewed a device to define a pseudo-true treatment parameter that incorporates the available information on average effects and between cluster heterogeneity. A transfer estimate of this kind would summarize the relevant evidence from the available experimental data and could be subject to additional (qualitative or quantitative) sensitivity analysis with respect to potential violations of the exchangeability assumption.

Rather than imposing strong assumptions necessary for identification of counterfactuals in a target location, our focus is on prediction. Alternatively, we can impose conditions under which that predictor is asymptotically unbiased estimator for a version of the problem in which sites are drawn at random from an infinite superpopulation and consistent for average effects at the target site. Conditions under which the bias from linear interpolation vanishes are discussed in Appendix (ref).

The empirical application in this paper concerns the effect of conditional cash transfers (CCT) to households on children's school attendance. The effect of CCTs was first evaluated in a large multi-site trial of the PROGRESA/OPORTUNIDADES program in Mexico, which was followed by implementations and additional RCTs in many developing and middle-income countries, often modeled after the PROGRESA study. After applying selection criteria we construct a data set of 640 sites, combining data from five studies in Mexico, Morocco, Indonesia, Kenya, and Ecuador to illustrate our approach. One non-technical contribution of this paper is to exploit cross-site variation within and across studies for extrapolation across populations, where we find that site heterogeneity at baseline predicts cross-study differences in post-intervention responses and conditional average treatment effects.

The problem of adapting empirical findings to new contexts allowing for unobserved heterogeneity is certainly not limited to estimation of discrete treatment contrasts but is also relevant to make more model-based estimates generalizable or comparable across settings. A fully nonparametric approach appears to be well-suited for the particular problem of a binary policy intervention, but can be seen as a stand-in for a more pragmatic estimation approach based on a more explicit model for the outcome of interest. For more “structural" approaches, it may be preferable to choose low-dimensional models of site heterogeneity that can be directly incorporated into the model, possibly motivated by economic theory or empirical regularities.

Literature

A conceptual framework for the problem of extrapolation of estimated treatment effects across heterogeneous sites was developed in the seminal article by HIM05. Using their terminology, we assume unconfounded locations, but specifically want to allow for (site-specific) model shifts (“macro effects"), i.e. shared heterogeneity in potential outcomes and treatment effects within each context. We propose a mechanism to incorporate information on pre-treatment outcomes at the cluster/site level when no treated units are observed in the population of interest.

Extrapolation of treatment effects was considered by various studies, including DPS21, Gec23, Mea22, NIW21, ACh22, and CSO23. DPS21 considered the problem of predicting treatment effects at target sites based on observed site-specific covariates. Gec23, Man20, and NIW21 derive bounds that account for selection effects at the individual level, allowing individual heterogeneity to be distributed differently across sites. Our focus is on site-specific heterogeneity, in particular we do not require the support of unobservables $(U_{ig}',V_g')'$ to be shared across sites for the approach to be useful. ACh22 consider robust extrapolation of treatment rules when there is no separate data on the target site, but the distribution of potential outcomes is in a neighborhood of that for the experimental population.

A separate question concerns the transfer performance of extrapolation methods. GSDP19 use data from two conditional cash transfer programs to evaluate extrapolation of empirical treatment rules. KXCAL18 identify attributes that exhibit a stable predictive relationship to the outcome of interest across environments. AFLW22 analyze the problem of assessing transfer performance, where model estimates from data in one domain are transferred to another, whereas this paper optimizes cross-domain model performance within the experimental sample. While our analysis is formally design-based conditional on the experimental sample (rather than assuming i.i.d. draws of contexts from a meta-population), a sampling-based interpretation similar to theirs is also possible. GHLMMMRS23 discuss optimal selection of experimental locations for extrapolation to other sites.

The work closest to ours is CSO23 who propose to use the distribution of pre-intervention outcomes for the target site together with post-intervention outcomes from the experimental locations to predict outcomes under a synthetic transferability condition. Their approach predicts policy effects based on the assumption that the policy shift affects outcomes through an index where the supports of pre- and post-intervention index values overlap in the target population. We consider a setting in which the policy intervention is binary and not equivalent to a shift in other observed covariates. Under that scenario the supports for pre- and post-intervention values for that index are disjoint, so that no subpopulation of the target site can be directly matched to post-intervention outcomes in the experimental sample. Our approach predicts counterfactuals conditional on pre-treatment outcomes alone, and is therefore complementary to theirs.

The working assumption of exchangeability between experimental and target sites is shared by conformal prediction methods (see VGS05 and LGRTW18). The focus of the present paper is on a point estimate that is informed by the experimental sample rather than inference, however under an exchangeability assumption our approach could in principle be combined with classical or conformal methods for inference with either asymptotic or finite-sample guarantees. Sensitivity of conformal inference with respect to departures from exchangeability was characterized by BCRT23. We do not explore the problems of inference or sensitivity analysis in this paper but leave this for future research.

It is also worth comparing our approach to other conceptual frameworks for aggregation of causal estimates across different populations: PBa14 gave explicit conditions for transportability of causal estimates across populations in terms of selection diagrams. One interpretation of our approach is the construction of a site-level covariate from baseline outcome data conditional on which potential outcomes are, to an approximation, mean-independent of selection. This paper also differs in the interpretation of transfer estimates, where our focus is on cross-population prediction of causal effects, rather than assuming the idealized conditions that would guarantee transportability in the strict sense.

Conceptually, the extrapolation problem also has some resemblance with the method of synthetic controls (AGa03,ADH10). However our approach is developed with a setting in mind where we do not have (typically aggregate) time series information on a treated unit and the “donor pool" of potential controls. Rather we assume that each site/context provides rich cross-sectional information, where a fraction of units is treated in a study population of sites, and we then predict treatment effects for the (yet untreated) target cluster. For that problem, Gun23 is most similar to our approach in that he proposes to use cross-sectional variation in micro-data to calibrate synthetic weights, however our approach differs in that rather than optimizing weights to match the distribution of baseline outcomes as closely as possible, we construct factors that are optimized to predict post-treatment outcomes based on the conditional distribution given unit-specific attributes. Shi22 uses a k-means algorithm to model unobserved heterogeneity in a problem with cluster dependence in treatment assignment.

In order to model site-specific conditional mean functions as random objects, we use tools from functional data analysis (see RSi05 and also WCM16 for a more recent overview), where function-to-function regression was analyzed by HMW03 and HMWY10,YMS11, and BCF17. Our approach is also related to functional principal components approaches for function completion/reconstruction based on partially observed functional data, where our setting corresponds more closely to that of sparsely sampled functions analyzed in YMW05, rather than the dense case considered by Kra15 and KLi20, although we assume that the number of points sampled for each curve (site) grows large. Since our focus is on cases in which only a modest number of trajectories is observed, the basis functions for our approach is constructed in as to be optimal for prediction, using both covariate and outcome data rather than separate principal components for covariate and outcome trajectories.

Generally our problem differs from function reconstruction in that our objective is to predict the difference between two curves, corresponding to conditional mean functions for either potential value, rather than the trajectory of the partially observed curve, so that the functional principal components of the conditional mean functions themselves do not generally have the best basis property for this particular task. Our problem differs from that of covariate adaptive reconstruction (JWa10,Lie19) in that we consider unit-specific covariates which correspond to coordinates of the random trajectories, rather than site-specific covariates that shift the distribution of the random curve. Prediction of scalar outcomes based on functional principal components was analyzed by CHa06 and HHo07.

Our focus is on prediction of the conditional average treatment effect as a function of covariates, and we derive a choice of basis functions that is optimal for that prediction task in a sense to be made more specific below. We show that our solution bears some resemblance with, but is distinct from Hot36's classical problem of canonical correlation analysis. For functional data, functional canonical regression has first been proposed by HMWY10 whose approach differs from the present paper in terms of the approach to regularization. We derive our approach from optimality considerations and establish a (to our knowledge novel) formal optimality result.

Interpreting “locations" at which random trajectories are evaluated as covariates or causal variables also requires a few subtle adjustments relative to the classical literature on functional data analysis. In particular, the covariate distributions may differ across sites, so nonparametric estimation of moments of the distribution of the random function requires some local reweighting and support conditions.

Organization

The remainder of the paper is organized as follows: we first give a formal characterization of transfer estimation as a statistical problem. We then determine the optimal finite-dimensional subspace of features of the baseline data, and propose nonparametric estimators based on the experimental sample. Asymptotic properties of those estimators, assuming the number of experimental sites grows large, are given in Appendix (ref). The approach is then illustrated using an application to predicting the causal effect of conditional cash transfer programs to new locations.

Problem Description

The population of interest consists of $G$ sites (“clusters"/“contexts"), where the $g$th site consists of $N_g$ units. Our focus is on the case in which there is a single target site $g^*$ in addition to $G-1$ experimental sites $g\in\{1,\dots,G\}\setminus\{g^*\}$. We also use the dummy variable $R_g\in\{0,1\}$ to indicate whether $g$ is an experimental location ($R_g=1$), or a target site ($R_g=0$).

There is a binary policy variable (“treatment") $D_{gi}\in\{0,1\}$ which acts at the level of the unit $i$ at site $g$, where we assume that the outcome of interest is determined only by the unit's own treatment status (SUTVA). Specifically, the unit is associated with potential outcomes $Y_{gi}(0),Y_{gi}(1)$, where the realized outcome is given by $Y_{gi}:=Y_{gi}(D_{gi})$. Furthermore, each unit is associated with a finite-dimensional vector $X_{gi}$ of attributes whose distribution is given by the p.d.f. $f_g(x)$ for cluster $g$, where we assume that the support $\mathcal{X}$ of $X_{gi}$ is a compact subset of $\mathbb{R}^d$. For the purposes of this paper $N_g$ will be treated as infinite, but the researcher only observes a finite random sample of $n_g$ units for each cluster.

Adapting notation from NIW21, we can represent potential outcomes as

for some unspecified mapping $y(\cdot)$ and potentially multi-dimensional unobserved individual and site-specific heterogeneity $U_{gi}$ and $V_g$. We first define key objects in terms of a superpopulation model in which $V_g$ and $U_{gi}$ are random draws from an unspecified distribution. Our statistical approach will be conditional on a fixed population of $G$ sites with heterogeneity $V_1,\dots,V_G$ without additional restrictions on how those sites were selected.\footnote{Previous work by Gec23 and NIW21 proposed strategies to address cross-site differences in the conditional distribution of individual heterogeneity $U_{gi}$, whereas our focus is on site-specific heterogeneity $V_g$. While $V_g$ could be included with the vector $U_{gi}$ as a matter of notation, the approaches in Gec23 and NIW21 require $U_{gi}$ to have the same support across sites, which can't be satisfied by variables $V_g$ that are shared by all units at the site. We therefore prefer to keep site-specific heterogeneity explicit in our notation.}

Using this notation we can write the conditional expectation of post-intervention outcomes at site $g$ for $D_{gi}=d$ as \[\mu_g(x;d)\equiv\mu(x;1,V_g):=\mathbb{E}[Y_{gi}(d)|X_{gi}=x,V_g]\] The site-specific conditional average treatment effect is given by \[\tau_g(x)\equiv\tau(x;V_g):=\mathbb{E}[Y_{gi}(1)-Y_{gi}(0)|X_{gi}=x,V_g]\] In particular, $\mu_g(x;d)$ and $\tau_g(x)$ are functions of site-specific unobserved heterogeneity $V_g$ and therefore random objects whenever $V_g$ is regarded as stochastic. For a given superpopulation $V_g\sim F_V$, we can also define the cross-site averages $\mu(x;1):=\mathbb{E}_{F_V}\left[\mu(x;1,V_g)\right]$ and $\tau(x):=\mathbb{E}_{F_V}\left[\tau(x;V_g)\right]$ of the CATE.

Our goal is prediction of $\tau_g(x)$ rather than consistent estimation, although under a more restrictive superpopulation framework and a linearity assumption, the prediction problem can also be cast as estimation of $\tau_g(x)$, see Appendix (ref) for a dicussion. We aim to predict model shifts

using the site-specific distribution of pre-intervention outcomes, $Y_{gi}(0)|X_{gi},V_g$.

Prediction of site-specific CATE therefore seeks to account for model shifts $\Delta\tau_{g}(x)$. Our method aggregates information on the first two moments of the distribution of conditional expectation functions (pre- and post-intervention) across sites and does not require that we can estimate either conditional mean function consistently for any individual site. In particular, we also discuss a version of our aproach for the case in which treatment assignment was randomized at the site level. In principle, the arguments behind our method can therefore also be extended to imputation of site-specific CATE for experimental sites when treatment was randomized at the site level, or the researcher only observes a moderate number of units for each site.

Functional Predictors

Our approach aims to extract predictive information regarding the unobserved site-specific heterogeneity $V_g$ from baseline (pre-intervention) outcome data. Since $V_g$ is shared among all units at the same site, not only the baseline outcome $Y_{gi}(0)$ of a target unit is predictive of the post-intervention out come $Y_{gi}(1)$ of that same unit, but the conditional distribution of pre-intervention outcomes given covariates for that site, $Y_{gi}(0)|X_{gi},V_g$, contains additional information regarding the model shifts ((ref)) and ((ref)). This is particularly plausible when the effect of the intervention is expected to be only incremental so that pre- and post-intervention outcomes behave similarly and depend on the same unit- and site-specific factors. Under this view of the DGP, unobserved site-specific heterogeneity $V_g$ in expected outcomes is not necessarily separable, but site effects will often manifest themselves in interactions between attributes and outcome variables.

In practice, the researcher may in addition directly observe site-specific measures e.g. of price variables or the cost of attending school, and our approach could then be viewed as addressing residual site-specific heterogeneity after practically feasible adjustments for observable covariates. We discuss this further in Section (ref) below.

Our approach models the conditional distribution of $Y_{gi}(0)$ given $X_{gi},V_g$ as functional data, which is then used to extract site-specific factors $m_{g1},\dots,m_{gK}$ (say) to predict a model shift $\Delta\mu_g(x;1)$ or $\Delta\tau_g(x;1)$. While other moments of the conditional distribution of baseline outcomes may reveal additional information regarding $V_g$, in this paper we restrict our attention to the problem of using only the conditional first moment of baseline outcomes \[\mu_g(x;0)\equiv\mu(x;0,V_g):=\mathbb{E}[Y_{gi}(0)|X_{gi}=x,V_g]\] as a predictor of $\Delta\tau_g(x)$.\footnote{In our application, the outcome of interest $Y_{gi}$ is a binary indicator whether a school-age child attends school, so that all higher moments of potential outcomes are known functions of $\mu_g(x;0)$, but in general higher-order conditional moments of $Y_{gi}(0)$ given $X_{gi}$ may also be predictive of the CATE at the target site. IKi21 propose efficient aggregation of noisy site-specific estimates of unconditional ATEs. In our leading scenario, cluster size $n_g$ is large relative to $G$ so that error in estimating $\mu_g(x;0)$ is asymptotically negligible. In the sparsely sampled case in which the number of units per site is not large, estimation error in $\mu_g(x;0)$ gives rise to similar efficiency considerations which we do not address in this paper. We also do not consider the use potentially predictive information form the marginal distribution of covariates $f_{X_g|V_g}(x|V_g)$. It is also possible to incorporate observable site-specific covariates into our approach, as discussed in Section (ref) below.}

For a target site $g^*$ drawn from a superpopulation, $V_{g^*}\sim F_V$, the best (lowest variance) predictor of $\tau_{g^*}(x)$ given $\mu_{g^*}(\cdot;0)$ is the conditional expectation function,

Since $\mu(\cdot;0,V_{g^*})$ is generally infinite-dimensional (unless all attributes $X_{gi}$ are discrete), completely flexible interpolation between sites is generally not feasible as a practical matter. Instead, we restrict our attention to predictors that are linear in $\mu_{g^*}(x;0)$,

for a square integrable function $\beta(x_1,x_2)$. That is, we can view a linear predictor as a regression adjustment over the unconditional CATE $\tau(x)$. Finding the kernel $\beta(x_1,x_2)$ corresponding to the best linear predictor is the classical functional linear regression problem (see RSi05, HMWY10, and BCF17).

Estimation of ((ref)) from a modest number of experimental sites requires substantial regularization. The dimension of the function generally equals the cardinality of $\mathcal{X}$, and the researcher may choose to work with approximations in an $S$-dimensional sieve space, e.g. using functional principal components or using B-splines, as in our implementation below. We propose to substantially reduce the dimensionality of this problem by constructing a subspace of $K<<S$ predictive features from $\mu_g(x;0)$ in a way that is optimal for prediction in a sense to be made more precise below.

Our approach differs from a conventional application of functional regression techniques in that rather than aiming for consistent estimation, we regard regularization via the choice of $K$ as fixed and instead aim at constructing those predictive features optimally. In our application, cross-validation recommends an approximation using a subspace of dimension as low as $K=2$, a level at which other regularization approaches may be difficult to interpret. Our approach still yields a best linear predictor given those constructed predictive features, whose construction and distributions can be reported and analyzed explicitly in any empirical application.

Fixed-Population Approach

We are interested in solving the functional prediction problem ((ref)) for situations in which the researcher wishes to extrapolate from existing experimental data and therefore has limited control or knowledge on how those sites had been selected. In such a scenario, it is generally implausible to assume a well defined sampling mechanism from a particular superpopulation, however defined. Instead, we follow a fixed-population (design-based) approach to the problem of extrapolating from experimental to target sites, where our statistical theory will regard the combined population of experimental and target sites as fixed, but only the role of the target site $g^*$ is regarded as random.

To be specific, we analyze the statistical properties of a predictor under the distribution defined by the following hypothetical protocol: in a first stage, we select sites to the experimental arm by drawing $G_1$ sites at random from $\{1,\dots,G\}$ uniformly and without replacement. The remaining locations are assigned the role of a target site, and we take $R_g\in\{0,1\}$ to be an indicator variable that equals one if $g$ is an experimental site, and zero otherwise. In a second step, individualized treatments $D_{gi}\in\{0,1\}$ are assigned at random to units, and the intervention is implemented according to that assignment at each experimental site $g$. Our main results concern the case of unit-level randomization at each experimental site, but we also discuss the case of site-level randomization separately. Finally in a third step we sample units uniformly at random at all sites and use the resulting sample to construct extrapolation estimates for the CATE at each target site.

Under this fixed-population experiment, the cluster-specific conditional average treatment effects $\tau_1(x),\dots,\tau_G(x)$ are nonstochastic, however the assignment $R_g$ of sites to the experimental role as well as the selection $D_{gi}$ of treated units within each experimental cluster are random. In particular, the cross-site average and empirical covariance of the functions $\mu_g(x;d)$ can only be estimated with error since even for units included in the sample, only one of the two potential outcomes $Y_{gi}(0),Y_{gi}(1)$ is observed. For the remainder of the paper we focus on the case in which there is a single target cluster in addition to $G-1$ experimental clusters.

We consider a transfer estimate $\hat{\tau}_{g^*,1,\dots,G}(x)$ for extrapolating from the sites $\{1,\dots,G\}\backslash\{g^*\}$ to $g^*$. Such a transfer estimate combines information on covariates and outcomes from the $G$ experimental and target sites to predict the CATE for the target site, $g^*$. We evaluate the statistical performance of such a transfer estimate in terms of the integrated mean-squared error (IMSE) under the resulting statistical experiment,

with a weight function $f_0(x)$ that has the properties of a p.d.f. and is chosen by the researcher. That function could be e.g. the uniform distribution on a compact set, or an estimate of the covariate distribution across the $G$ sites.

The fixed population approach is therefore used as a way of formalizing the researcher's problem who aims to produce a forecast that performs as well as possible on average for prediction among this fixed population of sites. The best feasible prediction under those circumstances is a parameter that is specific to the set of observable experimental and target sites. The resulting transfer estimator represents a summary of site-specific unobserved model heterogeneity that can be quantified based on the experimental sample and used to predict the treatment effect at the target site. This is analogous to a situation that would arise when the sample average treatment effect (SATE) is used to predict the treatment effect for an individual participant in an experimental trial on subjects that were not sampled at random from a well-defined population.

We derive theoretical properties of the approach using fixed-population asymptotics (see AAIW17), where approximations are derived under a sequence of finite populations along which the number of sites $G$ grows large. We analyze scenarios at which sites are either sampled densely, where $n_g\rightarrow\infty$ for each site $g$, or sparsely, where $n_g$ remains bounded. While for any given application, $G$ is obviously fixed, embedding the fixed-population prediction problem into such a sequence of statistical experiments allows to establish stochastic orders of magnitude for estimation errors as long as $G$ is sufficiently large for those approximations to be close.

Optimal Predictors for Model Shift

This section concerns the optimal choice of basis functions (features) for estimation of the linear projection problem ((ref)). Our approach is based on a representation of the random processes $\mu_{g^*}(x;0)$ and $\tau_{g^*}(x)$ for the target site $g^*$ in terms of orthogonal bases. To be specific, for a given pair of orthogonal bases $\phi_1,\phi_2,\dots$ and $\psi_1,\psi_2,\dots$ of square integrable functions, respectively, we can write

for each $g=1,\dots,G$. We use a fixed-population framework in which the target site $g^*$ is a random draw from the deterministic population $\{1,\dots,G\}$, so that $\mu_{g^*}(x;0)$, $\tau_{g^*}(x)$, and the corresponding coefficients $\left\{m_{g^*k},t_{g^*k}\right\}_{k=1}^{\infty}$ are stochastic.

Our approach then estimates a truncated version of this expansion for $\tau_{g^*}(x)$,

to approximate the CATE at site $g$ at a low order $K<<G$. For the scenarios we are envisioning in this paper, the number of experimental clusters is not very large, so $K$ should be thought of as fairly small. In fact, for our empirical application, cross-validation (with respect to cross-site prediction) suggests a value of $K$ equal to 2 or 3, depending on the exact specification. So rather than aiming for consistent estimation of $\tau_g(x)$, we view the use of the first few leading factors in the expansion ((ref)) as a method of improving over the unconditional forecast $\tau(x)$ in order to account for site-specific heterogeneity.

It is therefore all the more important to have theoretical guidance on how to choose the basis of that expansion optimally so as to prioritize those features in the data that will be most predictive for $\tau_{g^*}(x)$. The need to truncate the expansion for purposes of estimation stems from ill-posedness in the problem of predicting $\tau_{g^*}(x)$ based on trajectories $\mu_{g^*}(x;0)$. While other continuous regularization methods are available (see CFR07), an advantage of this finite-dimensional approximation is that it can be interpreted as a linear prediction of the CATE based on the first $K$ factors in an analogous expansion of the function $\mu_{g^*}(x;0)$ for arbitrary fixed values of $K$.

Our approach requires nonparametric estimation of the mean functions \[\mu(x;d):=\frac1G\sum_{g=1}^G\mu_g(x;d),\hspace{0.5cm}d=0,1\] and \[\tau(x;d):=\frac1G\sum_{g=1}^G\tau_g(x;d),\hspace{0.5cm}d=0,1\] as well as the covariance kernels

These objects can be interpreted as expectations and covariances, respectively, with respect to a random draw of a site $g^*$ from the discrete uniform distribution over $\{1,\dots,G\}$.

A standard representation of the random processes $\mu_g(x;0)$ and $\tau_g(x)$ in ((ref)) is the Karhunen-Lo\`eve expansion, which chooses the basis functions $\phi_1,\phi_2,\dots$ and $\psi_1,\psi_2,\dots$ as eigenfunctions of the respective covariance operators $H_{\mu\mu}(\cdot),H_{\tau,\tau}(\cdot)$, see RSi05 and RWi06. These bases of eigenfunctions ordered by their associated eigenvalues are also known as the functional principal components (FPC) of the random functions $\mu_g(x;0)$ and $\tau_g(x)$. At any finite order, an reconstruction of the function by its leading $K$ FPC is known to be optimal with respect to the mean-square error of approximation. However, our goal is to extract those features of $\mu_g(x;0)$ that are “most predictive" for the average of $\tau_g(X_{gi})$, which generally do not coincide with the FPC. We show that instead, that optimal choice can be described in terms of a singular value decomposition of an operator characterizing the covariance between $\mu_g(x;0)$ and $\tau_g(x)$.

Optimal Basis Functions

Our main objective is to determine the optimal finite-dimensional feature space for the baseline data in which to solve the prediction problem ((ref)). We regard the conditional mean functions $\mu_{g}(x;d)$ and $\tau_{g}(x)$ as random elements of the Hilbert space $L_2(\mathcal{X},f_0)$ ($L_2(\mathcal{X})$ henceforth) of square integrable functions with norm induced by the scalar product \[\langle \phi,\psi\rangle = \int\phi(x)\psi(x)f_0(x)dx\] where $f_0(x)$ denotes the weighting function introduced in ((ref)).

We also define integral operators $T_{\mu\mu},T_{\mu\tau}$ associated with the covariance kernels

for any square integrable function $\varphi$. The operators $T_{\mu\mu},T_{\tau\tau}$ are self-adjoint, whereas the adjoint of $T_{\mu\tau}$ is given by \[(T_{\mu\tau}^*\varphi)(x):=\int H_{\mu\tau}(x,x_1)\varphi(x_1)f_0(x_1)dx_1.\]

We now turn to the construction of an optimal $K$-dimensional basis for predicting $\tau_g(x)$ based on $\mu_g(x;0)$. For a collection of $K$ functions $\phi_1,\dots,\phi_K\in L_2(\mathcal{X})$, we let $P_K:L_2(\mathcal{X})\rightarrow\mathcal{H}_K$ denote the operator associated with orthogonal projection onto the closed linear subspace \[\mathcal{H}_K:=\textnormal{span}\left(\phi_1,\dots,\phi_K\right):=\left\{\sum_{k=1}^Ka_k\phi_k:a_1,\dots,a_K\in\mathbb{R}\right\}\] By the classical projection theorem (Theorem 2 on p.51 in Lue69) that projection is well-defined.

We then consider the predictors $BP_K\mu_g$ for $\tau_g$ on that subspace corresponding to linear operators $B:L_2(\mathcal{X})\rightarrow L_2(\mathcal{X})$, where we define $B$ via \[(B h)(x):=\int h(x_1;0)\beta(x_1,x)f_0(x_1)dx_1\] for any function $h\in\mathcal{H}$. We then let

denote the integrated mean-square error of prediction, minimized over the set of linear predictors using those $K$ functions. We restrict our attention to basis functions in the closed linear subspace $\mathcal{N}^{\perp}$, the orthogonal complement of the null space of $T_{\mu\mu}$, $\mathcal{N}:=\ker(T_{\mu\mu})$. This restriction is of no practical consequence since for any function $h$ in the null space of $T_{\mu\mu}$, $\textnormal{Var}(\langle \mu_g,h\rangle)=\langle h,T_{\mu\mu}h\rangle=0$. Considering any possible choices of $\phi_1,\dots,\phi_K\in L_2(\mathcal{X})$, we first give a lower bound on $IMSE_K$

See the appendix for a proof. The operators $T_{\mu\mu},T_{\mu\tau}$ are known to be compact if the corresponding covariance kernels are square-integrable, that is if the integrals $\int H_{\mu\mu}(x_1,x_2)^2dx_1dx_2$ and $\int H_{\mu\tau}(x_1,x_2)^2dx_1dx_2$ are finite. Since the operator $T_{\mu\mu}$ in the constraint is compact, there is no guarantee that the infimum will be attained by square integrable functions $\phi_1,\dots,\phi_K$. Intuitively, this ill-posedness stems from the fact that there may be functionals of $\mu_g(x;0)$ that have small variance across sites but are highly predictive with respect to $\tau_{g^*}(x)$. This problem bears some resemblance with functional canonical analysis, where HMW03 propose high-level conditions on the cross-correlation operator which would also be sufficient to guarantee that the infimum in ((ref)) is in fact attained at elements in $L_2(\mathcal{X})$.

If such a solution exists, it can be easily seen from the expression for $IMSE_K^*$ that the optimal basis functions for linear prediction are given by the solutions to the generalized eigenvalue problem

where we select the eigenfunctions $\phi_1^*,\dots,\phi_K^*$ associated with the $K$ leading eigenvalues $|\lambda_1|\geq|\lambda_2|\geq\dots$.\footnote{Note that while the self-adjoint operators $T_{\mu\tau}T_{\mu\tau}^*$ and $T_{\mu\mu}$ are both nonnegative, the generalized eigenvalue problem may have solutions associated with a negative eigenvalue.}

Our results allow for multiplicities of eigenvalues rather than requiring the ordering of $\lambda_1,\lambda_2,\dots$ to be strict. In that case, ((ref)) holds equivalently for any re-ordering or linear combination of eigenfunctions associated with the same eigenvalue. However any such transformations also yield the same minimum in ((ref)) and are therefore equivalent for the purposes of minimizing the IMSE for prediction.

Rather than imposing conditions for existence, we focus instead on a regularized version of the problem, where we then demonstrate that the solution to that problem is approximately optimal in the sense that they achieve an IMSE that can be arbitrarily close to $IMSE_K^*$ when the regularization parameter is sufficiently small. We discuss conditions for existence of a non-regularized solution to that problem separately in Appendix (ref).

Specifically, we consider the following generalized eigenvalue problem

where $a>0$ is a regularization parameter. We then let $\phi_{1a}^*,\dots,\phi_{Ka}^*$ be the eigenvectors corresponding to the $K$ largest eigenvalues (in absolute value), that is $|\lambda_{1a}|\geq|\lambda_{2a}|\geq\dots|\lambda_{Ka}|\geq|\lambda_{K+sa}|$ for each $s\geq1$, where we impose the normalization $\langle\phi_{ka}^*,T_{\mu\mu}\phi_{ka}^*\rangle=1$ for each $k=1,\dots,K$. In what follows, we also denote the operator $T_{\mu\mu a}:=T_{\mu\mu}+a \textnormal{Id}$.

We denote the integrated mean-square error of prediction using the basis from the regularized problem ((ref)) with \[IMSE_K^*(a):=\int\min_{B\in\mathcal{H}_K^*\times L_2(\mathcal{X})}\mathbb{E}\left[(\Delta\tau_{g^*}(x)-BP_K^*\mu_{g^*}(x;0))^2\right]f_0(x)dx\] where $P_K^*$ is the orthogonal projector onto $\mathcal{H}_K^*:=\textnormal{span}\left(\phi_{1a}^*,\dots,\phi_{Ka}^*\right)$. We show that the solutions to ((ref)) corresponding to the $K$ largest eigenvalues are approximately optimal as $a\rightarrow0$:

See the appendix for a proof. We can interpret this result as establishing an optimal finite-dimensional feature space for $\mu_g(\cdot;0)$ for predicting the conditional average treatment effect, up to a regularization bias that can be made small in terms of its impact on the IMSE of prediction.

Prediction of Model Shift

Given the proposed choice of $\phi_1^*,\dots,\phi_K^*$, we also state the projection of $\tau_g$ onto the optimal basis:

See the appendix for a proof. In particular, given the operators $T_{\mu\mu},T_{\mu\tau}$ defined at the population level, the optimal projection depends on the site-specific mean function $\mu_{g^*}(x,0)$ only through $K$ scalar features $(t_{g^*1},\dots,t_{g^*K})$ that can be estimated consistently as the number $n_{g^*}$ of observations in the target cluster grows large.

Incidentally, we can also confirm that each of the functions $\psi_{1a}^*,\dots,\phi_{Ka}^*$ is an eigenfunction of the operator $T_{\mu\tau}^*T_{\mu\mu a}^{-1}T_{\mu\tau}$ at the eigenvalue $\lambda_{ka}$:

Hence one interpretation of the approach is as an approximation based on the $K$ leading components of a singular value decomposition of the operator $T_{\mu\mu a}^{-1/2}T_{\mu\tau}^*$ on a suitably chosen linear subspace of $L_2(\mathcal{X})$: Consider the eigensystem $\phi_{1a}^*,\phi_{2a}^*,\dots$ solving ((ref)) at any nonzero value for the generalized eigenvalue $\lambda_{ka}$, and the corresponding functions $\psi_{1a}^*,\psi_{2a}^*,\dots$. By standard properties of eigenfunctions, these systems form a basis for the orthogonal complements of the null spaces $\ker(T_{\mu\tau}T_{\mu\mu a}^{-1/2})$ and $\ker(T_{\mu\mu a}^{-1/2}T_{\mu\tau}^*)$, respectively. Hence, using these bases as test functions, we can confirm that $\{\phi_{1a}^*,\phi_{2a}^*,\dots\}$, $\{\psi_{1a}^*,\psi_{2a}^*,\dots\}$, and $\{\sqrt{|\lambda_{1a}|},\sqrt{|\lambda_{2a}|},\dots\}$ represent a singular value decomposition of the operator $T_{\mu\mu a}^{-1/2}T_{\mu\tau}^*$ where \[(T_{\mu\mu a}^{-1/2}T_{\mu\tau}^*h)(s) = \sum_{k=1}^{K^*}\sqrt{|\lambda_{ka}|}\phi_{ka}^*(s)\langle\psi_{ka}^*,h\rangle\] for any $h\in L_2(\mathcal{X})$.

Comparison to Existing Approaches for Functional Regression

We briefly discuss how this approach compares to existing methods in the literature on functional regression with a functional response.

While the basis functions $\phi_{1k},\phi_2,^*,\dots$ in our analysis are derived from optimality considerations, the procedure we arrive at has a close resemblance to canonical correlation analysis which has previously been proposed for functional regression problems by HMWY10, see also LMS93. Our results differ in that for one the basis $\phi_1^*,\dots,\phi_K^*$ is formally shown to be optimal for the linear prediction problem considered here. Moreover, the canonical variates need not be ordered according to the eigenvalues $\lambda_k$ which we show to be the relevant ordering for the IMSE-optimal choice among the eigenfunctions $\phi_1^*,\phi_2^*,\dots$.

To address the potential non-existence of an unregularized solution to ((ref)), HMW03 and HMWY10 impose high-level conditions on the cross-correlation operator to ensure existence (see Proposition 4.2 in HMW03). Since our focus is on prediction, we focus instead on the achievable IMSE, allowing for the possibility that unregularized canonical variates need not be well-defined. This approach parallels the analysis of CEGR08 who consider estimation of the largest canonical correlation between two $L_2$ processes and show that this scalar parameter can be approximated arbitrarily closely via regularized canonical correlation analysis.

YMS11 propose regression based on a singular value decomposition of the operator $T_{\mu\tau}$ rather than $T_{\mu\mu}^{-1/2}T_{\mu\tau}$, \[(T_{\mu\tau}h)(s):=\sum_{k=1}^K\sqrt{\nu}_k\langle\zeta_k^*,h\rangle\xi_k(s)\] Similarly, ROg07 propose functional partial least squares for functional regression. A finite-$K$ expansion based on spectral analysis of $T_{\mu\tau}$ has no known optimality properties but elegantly sidesteps the problem of inverting $T_{\mu\mu}$ and therefore works under weaker conditions and is numerically stable in the absence of regularization.

Another important approach proposed by BCF17 who directly minimize the mean-square error of prediction in a functional linear regression model, subject to a nuclear norm penalization of the projection operator $B$. The particular appeal of that approach is that it offers a “one-stop" approach towards regularization with a single tuning parameter, and directly optimizes the in-sample predictive performance subject to that penalty. Their approach assumes that $B$ is a Hilbert-Schmidt (kernel) operator which is not guaranteed under our assumptions. Their approach is also designed towards delivering a consistent estimator for $B$ in a setting where $G$ is large.

Our focus is instead on heavily regularized but interpretable solutions $B_{a,K}$ for moderate values of $G$, where the singular value representation delivers a sparse representation of the operator in terms of a functions of $x$. The estimated scores can then be used to assess whether the target site is comparable to the experimental sample in terms of the most predictive features identified by the method. The extrapolated CATE can be interpreted as a best linear predictor given the estimated basis functions, and regularization bias results in a potentially suboptimal (with respect to the IMSE), but ultimately valid construction of features from $\mu_g(x;0)$. As BCF17 point out, ridge regularization also yields more stable predictions in the presence of poorly separated eigenvalues than a truncation of the spectral expansion at a finite dimension, so if the eigenvalue $\lambda_{K}$ at the chosen cutoff is not well separated from $\lambda_{K+1}$, the resulting potential instability of predictions should be flagged when reporting estimation results.

Assumptions for Estimation

We next formalize the identifying conditions which are adapted from HIM05. We depart from their main framework in two substantial ways: for one our design-based approach treats experimental and target sites as random draws from a finite population of sites. Moreover, we also consider a version of the problem in which baseline data on pre-treatment outcomes for the target site are available and are to be used to predict site-specific “macro" effects. We highlight how this affects the interpretation of the assumptions on the assignment mechanism. While our derivation of optimal predictors in section (ref) is directly in terms of high-level properties of covariance operators, the following assumptions are maintained to establish asymptotic rates for estimates in Appendix (ref).

We assume throughout that for each cluster $g=1,\dots,G$ the researcher observes a sample of $n_g$ units that are drawn independently and uniformly at random from $\{1,\dots,N_g\}$, and also independently of potential values and unit attributes. For notational convenience our results will be stated for the case that the observed number of units is the same for each site, $n_g\equiv n$ for $g=1,\dots,G$. For each experimental site, we assume that selection of units into treatment is based only on observables $X_{gi}$,

This condition is met if $D_{gi}$ was assigned at random as part of a randomized controlled trial (RCT) at each experimental site, and it captures the idea of extrapolating from a collection of internally valid estimates of site-specific causal effects to a new site. In a practical application the set of confounders $X_{gi}$ may differ from the conditioning variables chosen by the researcher for define the relevant conditional average treatment effect, however for expositional clarity we only consider the case in which the conditioning variables are the same. It is also possible to adapt our approach to the case of randomization at the cluster level, $D_{gi}\equiv D_g$ for all $i=1,\dots,n_g$, see Appendix (ref) for a brief discussion.

Furthermore, we assume that among the $G$ sites, the $G-1$ experimental locations were selected independently of potential values, conditional on observable covariates:

This assumption is strengthened version of Assumption 2 in HIM05 and describes an idealized observational protocol that rules out systematic ex-ante site selection bias. It can be seen immediately that under this condition, for a randomly selected experimental site $\tilde{g}$ with $R_{\tilde{g}}=1$, $\left(Y_{\tilde{g}i}(0),Y_{\tilde{g}i}(1),X_{\tilde{g}i}\right)\stackrel{d}{=}\left(Y_{g^*i}(0),Y_{g^*i}(1),X_{g^*i}\right)$, where “$\stackrel{d}{=}$" denotes equality in marginal distributions. Therefore, Assumption (ref) implies that experimental and target sites are exchangeable, the fundamental assumption in the literature on conformal prediction (VGS05 and LGRTW18).

In practice, we do not expect that assumption to be an accurate description on how experimental (study) and target sites were selected. Rather, in the absence of additional knowledge regarding site selection, this auxiliary assumption defines a pseudo-true parameter, which aggregates estimates from experimental sites into a “best" prediction for the target population. The resulting transfer estimate should therefore be interpreted as a summary of the directly quantifiable relevant information from previous experiments, which could be subject to additional (qualitative or quantitative) sensitivity analysis with respect to suspected violations of that exchangeability condition (see e.g. BCRT23 for the problem of conformal prediction).

For the next assumption, we define the site-specific propensity score as \[p_g(x):=\mathbb{P}(D_{gi}=1|X_{gi}=x)\] We require that the supports of covariates overlap, both between treated and control units, as well as across the sites $g=1,\dots,G$.

The role of this assumption is to ensure that conditional moments of either potential value are identified and can be estimated consistently across sites. While Assumption (ref) does allow for experimental and target sites to differ in terms of the distribution of observables, we require that the site-specific supports overlap, potentially after trimming non-overlapping regions in the covariate space as suggested in HIM05. This assumption also does not cover site-specific aggregate covariates that may serve as additional predictors as analyzed by HIM05 and DPS21. Randomization at the level of the site would not satisfy the support condition on the site-specific propensity score and therefore requires a different approach which is discussed in Appendix (ref). Additional adjustments for site-specific variables may be possible, but would also be severely constrained by the small number of observable sites. While our focus is on the optimal use of cross-sectional information for extrapolation, we briefly discuss how to incorporate site-level covariates in Section (ref).

Nonparametric estimation of the first two conditional moments of potential values $Y_{gi}(d)$ given attributes $X_{gi}$ requires additional moment and smoothness conditions, where we specifically assume the following:

To avoid additional notation, we do not explicitly discuss the case in which some components of $X_{gi}$ may be discrete. With the exception of part (c), the conditions in Assumption (ref) are commonly assumed for nonparametric estimation of conditional moments, see e.g. Han08. Notice also that we effectively need to be able to estimate conditional moments separately for each site, and therefore require these conditions to hold uniformly over $g$. In the absence of covariate shifts, i.e. if the distribution of covariates $f_g(x)$ or propensity score $p_g(x)$ did not vary over $g$, this issue could be avoided (see YMW05), however we do not find such an assumption plausible for the problem considered here.

Implementation for Densely Sampled Clusters

The representation in Corollary (ref) motivates an estimator of the form \[\hat{\tau}_{g^*}(x):=\hat{\tau}(x) + \sum_{k=1}^{K}\hat{\tilde{t}}_{g^*k}\hat{\psi}_{ka}(x)\] where $\hat{\tau}(x):=\hat{\mu}(x;1)-\hat{\mu}(x;0)$, $\hat{\tilde{t}}_{g^*k}=\langle\hat{\mu}_{g^*},\hat{\phi}_{ka}\rangle$ for a nonparametric estimator $\hat{\mu}_g$ of $\mu_g(x;0)$, and the basis functions $\hat{\phi}_{1a},\dots,\hat{\phi}_{Ka}$ are obtained by solving an empirical analog of the generalized eigenvalue problem ((ref)).

Here we develop our approach for the case of densely sampled clusters, $n\rightarrow\infty$, separate results for the setting with sparse samples are given in Appendix (ref). In contrast to the densely sampled case, that approach requires that site-specific covariate distributions $f_g(x)$ are either known or can be estimated consistently, which does in general not allow those distributions to be fully nonparametric.

We estimate $\mu(x;d):=\mathbb{E}[\mu_{g^*}(x;d)]$ and $H(x_1,x_2;d_1,d_2):=\textnormal{Cov}(\mu_{g^*}(x_1;d_1),\mu_{g^*}(x_2;d_2))$ using nonparametric estimators $\hat{\mu}(x;d)$ and $\hat{H}(x_1,x_2;d_1,d_2)$. While our theory is not restricted to one particular choice of nonparametric estimators, following YMW05 we give results for local linear estimators: For each experimental cluster, let

with nonparametric weights \[w_{gi}(x;d):=1\hspace{-2.5pt}\textnormal{l}\{D_{gi}=d\}K\left(\frac{X_{gi}-x}{h}\right).\] Here, the notation “$\arg_{b_0}\min_{b_0,b_1}$" corresponds to the first component vector of the joint argmax of a function with respect to $b_0,b_1$.

Here, $K(u)$ is a kernel function with standard properties (see Assumption (ref) in the Appendix for formal conditions on $K(\cdot)$), and the bandwidth $h>0$ is chosen according to sample size $G,n$, the dimension of $X_{gi}$ and assumed smoothness of the estimands. We also let

where \[H_{gij}(x_1,x_2,\mathbf{b}):= \left(Y_{gi}Y_{gj}-b_0^{(g)}-b_{11}^{(g)}(X_{gi}-x_1) - b_{12}^{(g)}(X_{gj}-x_2)\right)^2.\] We then construct

In principle, the bandwidth could be chosen differently for estimation of $\hat{\mu}(x;d)$ and $\hat{H}(x_1,x_2;d_1,d_2)$, however in our theory in Appendix (ref), the optimal rate turns out to be the same for either estimator in the densely sampled case. Apart from kernel-based approaches, other possible methods include series estimators, random forests, or neural networks. The choice of nonparametric estimator will generally depend on the support of the covariates and other practical considerations.

This estimator is an average of separate local linear estimators for each of the $G-1$ experimental clusters, in a departure from the approach in YMW05 who propose a local linear estimator based on the pooled data from all $G-1$ clusters. There are two reasons for a different approach in the densely sampled case: for one we do not assume that attributes (“positions") are sampled from the same distribution in all clusters, but sites may differ in the distribution of $X_{gi}$. We furthermore assume “dense" samples from a small number of clusters, whereas they consider scenarios in which $n$ is small, but $G$ grows large. In our setup, cluster-specific moments can be estimated consistently, whereas between-cluster variation is the dominant source of estimation noise due to small $G$. That source of estimation error would be amplified in a nonparametric regression step, so our approach seeks to avoid that potential problem.

To describe the estimator for the basis functions $\hat{\phi}_1,\dots,\hat{\phi}_K$ let \[\hat{H}_{\mu\mu}(x_1,x_2):=\hat{H}(x_1,x_2;0,0)\hspace{0.5cm}\textnormal{and }\hat{H}_{\mu\tau}(x_1,x_2):=\hat{H}(x_1,x_2;1,0)-\hat{H}(x_1,x_2;0,0).\] In analogy to the definition for the operators $T_{\mu\mu}$ and $T_{\mu\tau}$, we can construct the estimators

for any square integrable function $h$, and let $\hat{T}_{\mu\mu a}:=\hat{T}_{\mu\mu} + a\textnormal{Id}$.

In order to estimate the eigenfunctions $\phi_{1a}^*,\phi_{2a}^*,\dots$, we solve the generalized eigenvalue problem ((ref)) after replacing the operators $T_{\mu\tau},T_{\mu\mu}$ with their estimates as defined above. Specifically, we can find the functions $\hat{\xi}_{1a},\dots,\hat{\xi}_{Ka}$ solving the eigenvalue problem

and that are associated with the $K$ largest eigenvalues in $\hat{\lambda}_1\geq\hat{\lambda}_2\geq\dots$. We then solve for

Since $\hat{T}_{\mu\mu}$ is a nonnegative (nonnegative definite) operator and $a>0$, the operator on the left-hand side of ((ref)) is Hermitian and compact, and the inverse problem ((ref)) is well-posed. To implement the procedure we use linear sieve approximations to the eigenfunctions, which converts ((ref)) into a finite-dimensional eigenvalue problem.\footnote{See e.g. RSi05, chapter 8.4.2.}

We then construct $\hat{\psi}_{ka}$ by applying the estimator of $T_{\mu\tau}^*$ to the estimated eigenfunction $\hat{\phi}_{ka}$, \[\hat{\psi}_{ka} (x):=\left(\hat{T}_{\mu\tau}^*\hat{\phi}_{ka}\right)(x)\equiv\int\hat{H}_{\mu\tau}(s,x)\hat{\phi}_{ka}(s)f_0(s)ds\] for $k=1,\dots,K$. Using these estimates, we then obtain \[\widehat{\tilde{t}}_{g^*k}:=\langle \hat{\mu}_{g^*},\hat{\phi}_{ka}\rangle\] Substituting this into the formula from Corollary (ref), our estimate of the conditional ATE $\tau_{g^*}(x)$ is \[\hat{\tau}_{g^*}(x)= \hat{\tau}(x) + \sum_{k=1}^{K}\widehat{\tilde{t}}_{g^*k}\hat{\psi}_{ka}(x)\]

Appendix (ref) gives convergence rates for these estimators both for densely and sparsely sampled sites. Specifically, assuming equal numbers of cross-sectional observations for each site, $n_g\equiv n$, Theorem (ref) gives the rate \[r_{Gn} = \frac1{G} + h^2 + \left(\frac{\log n}{Gnh^d}\right)^{1/2}\] for the preliminary nonparametric estimators of mean and covariance functions if sites are densely sampled ($n_G\rightarrow\infty$ as $G\rightarrow\infty$) and treatment is randomized among units in each site. If treatment is instead randomized at the site level, the approach to estimating the covariance function $H(x_1,x_2;d_1,d_2)$ has to be modified as discussed in the appendix, resulting in a rate \[r_{Gn} = \frac1{\sqrt{G}} + h^2 + \left(\frac{\log n}{Gnh^d}\right)^{1/2}\] for the densely sampled case. For sparsely sampled sites ($n_G$ bounded), $H(x_1,x_2;d_1,d_2)$ can still be estimated consistently as $G\rightarrow\infty$ by pooling observation pairs across sites. Convergence for eigenfunctions and the IMSE of prediction depends on the asymptotic rate of estimation of the covariance function $H(x_1,x_2;d_1,d_2)$, \[r_{Gn} = h^2 + \left(\frac{\log G}{Gh^{2d}}\right)^{1/2}\] whereas the rate for estimating the conditional mean function $\mu(x_1;d_1)$ is faster for reasonable bandwidth choices.

Given these preliminary rates, Theorem (ref) gives a rate \[\|\hat{\phi}_{k}-\phi_{ka}\|=O_p\left(a^{-3/2}r_{Gn}\right)\] for estimation of the eigenfunctions, and Corollary (ref) shows that the IMSE of prediction using the estimated basis function is

Appendix (ref) also provides comparable rates for nonparametric estimation of mean and covariance functions using B-splines instead of local linear regression. When the eigenvalues of the corresponding population problem ((ref)) are not simple, the eigenfunctions $\phi_k$ are only estimated up to a data-dependent rotation of each eigenspace associated with multiple eigenvalues. Since any such transformation yields the same value for the problem of minimizing the IMSE of prediction ((ref)), this does not affect the rate at which the IMSE is minimized in ((ref)).

Since the transfer estimate is always a linear projection on the constructed features $\hat{\phi}_{k}$, these rates illustrate how fast the quality of the prediction improves as we approximate the optimal basis functions $\phi_k^*$ more closely. In general, that approximation requires the number of sites $G$ to be not too small, especially if treatment was not randomized within each site. The difference in rates between the densely and sparsely sampled cases also illustrates how a larger number of cross-sectional observations $n_g$ within each site can be leveraged to retrieve the optimal predictors more accurately, although in practice typically the number of sites $G$ is the main limiting factor.

Site-Specific Covariates

A natural extension of the main framework concerns site-specific covariates $W_g$ which may be observed in addition to the unit-level attributes $X_{gi}$. In this section we sketch a conceptual extension to our approach under the assumption that these covariates satisfy unconfoundedness conditions analogous to those for $X_{gi}$. When the number of experimental sites is not very large, controlling nonparametrically for a significant number of site covariates is generally not feasible in practice, so we consider this extension to be primarily of theoretical interest. For the purposes of this section, we also regard the $G$ sites as random draws from a superpopulation in order to be able to define conditional expectations given the covariate $W_g$ in a meaningful way.

To be specific, we consider a version of the original problem, where Assumption (ref) is changed to \[D_{gi}\bot\hspace{-6pt}\bot(Y_{gi}(0),Y_{gi}(1))|X_{gi},W_g,R_g=1\] and Assumption (ref) is strengthened to assume that $g^*$ is drawn independently of $Y_{g^*i}(0),Y_{g^*i}(1),X_{g^*i}$, and $W_{g^*}$. Assuming that the $g$th cluster represents a random draw from a superpopulation, we can define the conditional expectation \[\mu(x,w;d):=\mathbb{E}[Y_{g^*i}(d)|X_{g^*i}=x,W_{g^*}=w]\] and covariance function

where expectations are with respect to the joint distribution of potential values, attributes, and $W_g$ in that superpopulation.

We can then apply the previous method conditional on $W_{g^*}=w_{g^*}$, where we replace the unconditional mean function $\mu(x;d)$ with an estimate of estimate $\mu(x;w_{g^*};d)$, and form the analogs of the covariance operators $T_{\mu\mu}$ and $T_{\mu\tau}$ from estimates of the conditional covariance function $H_{d_1d_2}(x_1,x_2;w_{g^*})$. The conditionally optimal basis functions $\phi_1^*,\dots,\phi_K^*$ are then obtained from an eigenanalysis of the conditional covariance operators given $W_{g^*}$. Such an approach would effectively amount to a regression adjustment for the mean and covariance functions for $\mu_{g^*}(x;d)$ with respect to $W_{g^*}$.

For modest values of $G$, the scope for fully nonparametric adjustments to site-specific covariates is fairly limited for practical purposes, in contrast to “micro" (unit-specific) covariates where our approach can leverage the size of the cross-sectional sample for each site to construct approximately optimal adjustments to estimates for the CATE. DPS21 used machine learning methods to adjust (unconditional) ATE estimates for site-specific covariates, however a fully nonparametric site-specific adjustment to the estimated CATE poses greater challenges given realistic sample sizes.

Empirical Application

We illustrate our approach with an empirical application to the estimation of the effect of conditional cash transfers on children's school attendance. In this literature, a conditional cash transfer is a recurring grant paid to an eligible household that is explicitly linked to a child attending school or other household decisions the policy maker wants to encourage, with transfer amounts of the order of 5-20 percent of average household consumption in the target population. In 1998-99 the government of Mexico conducted a large-scale randomized trial during the roll-out of the PROGRESA/OPORTUNIDADES program (Sch04 and TWo06), and similar programs have subsequently been implemented in over 50 other countries (see BHKO17 for a recent summary).

Data

We combine samples from PROGRESA with four additional randomized studies that were conducted in Indonesia (Program Keluarga Harapan (PKH), see Ala11 and CHOPSS20), Morrocco (Tayssir, see BDDDP16), Kenya (Kenya CT-OVC, CT_OVC12), and Ecuador (Bono de Desarrollo Humano (BDH), ESch12).\footnote{These studies were selected according to ease of access to the underlying microdata, where we excluded one additional study from Colombia (Subsidios Condicionados a la Asistencia Escolar, BBLP11) due to our inability to reconstruct baseline attendance data from the replication package.} Each of these field trials was a multi-site study conducted by the national government, where participants were recruited from a previously selected sample of clusters (schools, villages, or other comparable unit). In each study, clusters were drawn from a subset of the major administrative regions in each of these countries.

It should be noted that there were substantial differences in the exact design of the incentive between these five studies. In particular, Progresa and PKH explicitly make part of the transfer dependent on school attendance, whereas Tayssir experimented with a nudge rather than a strictly conditional transfer. For the remaining two studies in Kenya and Ecuador, cash transfers were unconditional. We deliberately pool the sites to replicate a realistic scenario for which a policy as been adapted to local circumstances, due to practical constraints and the policymaker's preferences.

Our main focus is on leveraging cross-site variation within each multi-site trial to extract predictive information on site-specific heterogeneity in the CATE. The five study populations in Mexico, Indonesia, Morocco, Ecuador, and Kenya are likely systematically different in terms of many factors that cannot be modeled explicitly, such as the local educational system, the chosen target population within the geographic reach of the study, the specific manner in which the transfer scheme was implemented, etc. Nevertheless, sites also vary substantially within each study, e.g. according to travel distance to urban centers or secondary school, or whether the language of instruction is widely spoken in the community. Hence, some communities in the heterogeneous pool of clusters in, say, Mexico, may still be sufficiently similar to a target location in Morocco or Indonesia in terms of the predictive attributes, as determined by our method. We will assess to what extent between-study variation can be predicted from between site variation on a more disaggregated level.

We retain all observations of households that met the eligibility criteria for the program, and for whom we can reconstruct measures of school attendance and per capita household expenditure at baseline and follow-up, along with children's age and gender, and the household head's level of education. For school attendance we use self-reports from baseline and follow up household surveys rather than data from school records or random checks which were only collected for some of the studies used in our analysis. After dropping households with incomplete data and locations with fewer than 15 school-aged children, we obtain a sample of 640 clusters (sites) with average cluster sizes ranging from 18 (PKH, Indonesia) to 47 (PROGRESA, Mexico) and 51 (BDH, Ecuador). PROGRESA and TAYSSIR (Morocco) contribute the largest number of clusters (297 and 238, respectively) compared to 50 for PKH (Indonesia), 31 for BDH (Ecuador), and 24 for CT-OVC (Kenya). For the purposes of this analysis we assign equal weight to each cluster. Of those clusters, 434 were treated, the remaining clusters were in the control group.

Implementation

We compare our approach across three different prediction tasks - as a benchmark, we report some results for the in-sample fit, with $\mu(\cdot)$ and $\textnormal{H}(\cdot)$ and resulting basis functions $\phi_k,\psi_k$ estimated from the full data set. We then consider cross-site prediction where for a given target site $g^*$, the basis functions are estimated from the remaining $G-1$ sites, and the transfer estimate is obtained by estimating the principal scores $m_{g^*1},\dots,m_{g^*K}$ from the baseline for the target site. Finally, we perform cross-study extrapolation, with the predictive model estimated from data excluding all other sites from the study that included the target site, for example predicting the outcome at a Progresa site using only data from sites in the remaining four studies.

Given the small to moderate cluster sizes, we choose an estimation approach suited to sparsely sampled functional data, see also Appendix (ref). The main difference to the densely sampled case is that the cluster-specific covariate distribution $f_g(x)$ for the weights in ((ref)) and ((ref)) cannot be estimated nonparametrically. We make the simplifying assumptions that gender and age are independent of location and household per capita expenditure, and per capita expenditure follows a log-normal distribution within each cluster, which we then estimate parametrically.

The setting also differs from the idealized setup discussed in the theoretical sections of the paper in that there is baseline data available for each experimental cluster. Furthermore, in each of theses studies, treatment was randomized at the cluster level. We therefore construct predictors from the observed baseline data for $\mu_{gt}(x;0):=\mathbb{E}[Y_{git}(0)|X_{git}=x]$ at $t=0$, which are then used to predict conditional expectations $\mu_{gt}(x;1):=\mathbb{E}[Y_{git}(d)|X_{git}=x]$ for $d\in\{0,1\}$ and $t=1$. The covariance operators between $\mu_{g0}(x;0)$ and $\mu_{g1}(x;1)$ or $\mu_{g1}(x;0)$ are then estimated using the treatment and control clusters in the experimental population. We first consider the problem of predicting post-treatment outcomes from baseline outcomes in the treated clusters, where we can validate predictions directly against the observed data at the site level. We then implement the algorithm for predicting conditional average treatment effects, which are not directly observed at the cluster level for any of the experimental sites.

Given the limited number of distinct sites, and also in order to apply consistent variable definitions across studies we restrict the unit-specific covariates $X_{gi}$ to four variables, the child's gender, the child's age in years, enrollment status at baseline, and log per-capita household expenditure. We also restrict the estimators for $\mu(\cdot)$ and $\textnormal{H}(\cdot)$ to be additively separable in covariates, where we flexibly dummy out gender and age in years, and use B-splines of degree 2 to model variation with respect to log expenditure. Tuning parameters are chosen using cross-validation across clusters, where we separately target the integrated mean-square error of estimating the mean and covariance functions to determine the bandwidths for local linear regression, and the mean-square error for cross-cluster prediction for the regularization parameter $a$ in ((ref)).

Results

We first report results for prediction of the model shift in post-intervention outcomes $\Delta\mu_{g}(x;1):=\mu_g(x;1)-\mu(x;1)$ using the estimated IMSE-optimal predictors from ((ref)), which were estimated using only the 434 treated sites. We assess their performance as predictors at the level of the individual site as well as after aggregating sites within each study. The number of knots for B-spline approximations was determined using (leave-one-site-out) cross-validation, targeting the mean function $\mu(x;d)$ and covariance function $H(x_1,x_2;d_1,d_2)$, respectively. The ridge parameter $a$ was chosen based on estimated cross-site predictive performance, and cross-validation also suggests that for this application the optimal number of basis functions is $K=2$.

Table (ref) reports the correlation coefficient between the predicted model shift for the average effect at site $g$, $\widehat{\Delta\mu_{g1K}}:=\frac1{n_g}\sum_{i=1}^{n_g}\sum_{k=1}^K\hat{t}_{gk}\hat{\psi}_k(x_{gi})$ with its post-hoc empirical counterpart, $\widehat{\Delta\mu_{g1}}:=\frac1{n_g}\sum_{i=1}^{n_g}(Y_{gi1}-\hat{\mu}_1(X_{gi}))$. A natural alternative strategy would be to predict post-intervention outcomes using separate regression estimates stratified by average pre-intervention outcomes. In the first column we therefore report correlation coefficients with the corresponding predictors as a benchmark, where sites were binned into three groups of equal size (terciles) according to average enrollment at baseline.

According to our results, optimal basis functions result in substantially more precise predictions relative to binned estimates and standard FPC, where gains are largest for the first two basis functions, and then plateau for 3 or more components. For example for cross-site prediction, we find a correlation coefficient of around $0.36$ (corresponding to an R-square of $0.13$) after using only the leading baseline function ($K=1$), which still gradually improves as additional terms are included. For $K$ larger than 5 or 6, terms are fairly noisily estimated and therefore do not lead to substantial additional improvements. As expected, the strength of correlation for cross-study extrapolation is lower than for cross-site prediction, but still substantial. Stratified estimation by pre-intervention levels of outcomes does not appear to extract much predictive information at all, suggesting that the gains observed for our estimator exploit information on how outcomes vary together with covariates at each site.

In Table (ref), we compare cross-site averages of predictions, where we let $\mathcal{G}_s$ denote the subset of $\{1,\dots,G\}$ corresponding to sites that were part of study $s=1,\dots,5$. For each study $s$ we then compare $\frac1{|\mathcal{G}_s|}\sum_{g\in\mathcal{G}_s}(\hat{\mu}_{g1K}-\hat{\mu}_1)$ to their “realized" empirical counterparts, $\frac1{|\mathcal{G}_s|}\sum_{g\in\mathcal{G}_s}(\hat{\mu}_{g1}-\hat{\mu}_1)$. We find that the predicted average outcomes reflect some of the systematic differences, although especially for BDH and CT-OVC, the numbers and sizes of clusters are smaller, so results are likely noisier than for the first three studies. It should also be noted that the baseline outcome $Y_{gi0}$ is already included as a control for post-intervention outcomes in the specification of $\mu_1$. Without controlling for state-dependence at the individual level (not reported here), the correlation between pre- and post-intervention outcomes at the site-level is substantially stronger, but the relative comparison between using baseline averages as the “naive" predictor and prediction using $K$ estimated basis functions is qualitatively similar.

We next repeat the same analysis using the respective functional PC for $\hat{\mu}_{g0}(x)$ and $\hat{\mu}_{g1}(x)$, see Table (ref). Since the general patterns of school attendance as a function of child and household attributes were unlikely to have shifted fundamentally between baseline and follow-up, and the effect of the intervention was sizeable but incremental, we should expect the functional PC for the baseline to be fairly closely aligned with those at follow-up, and therefore perform very well as predictors for post-intervention outcomes. This is confirmed by the quantitative results, where performance is very similar to the IMSE optimal predictors, likely within or close to the margin of error, although we do not formally quantify estimation error for these results.

In Figure (ref), we report estimates of the leading two leading optimal basis functions for predicting conditional post-intervention outcomes. These basis functions do not appear to vary much with income, so we plot $\phi_1,\phi_2$ only as functions of gender and age alone. Since post-intervention outcomes are also observed at all treated sites, we also plot the conditional mean square error for predicting the post-intervention response for those sites, where $IMSE_0$ corresponds to the case in which we use the unadjusted cross-site average as a predictor, and $IMSE_k$ for the prediction using the first $k$ basis functions $\phi_1,\dots,\phi_k$ as predictors. While the predictors appear to be responsive to differences in enrollments at young and old ages, most of the improvement in the forecast is for enrollment at ages 12 and above, where (within and across site) variation is generally highest. Most of the improvement in the conditional forecast results from including the first two factors, whereas additional predictors lead to a significant deterioration of the forecast at lower ages. This is in line with the number $K=2$ of factors selected by cross-validation.

One important question is whether these features constructed based on their predictive power capture systematic differences between the five different study countries (Mexico, Morocco, Indonesia, Kenya, Ecuador), but also whether there is substantial overlap between those populations. The latter is especially important since we use extrapolation that is linear in those features $\phi_k$. Figure (ref) plots the estimated scores corresponding to the leading two basis functions, $\hat{m}_{g1},\hat{m}_{g2}$, for each site. To visualize differences in the factor loadings between the five countries included in our analysis, we also plot study-specific variance ellipses corresponding to a 80 percent confidence set for jointly normal variates. We can see that while there is substantial overlap in the support, their distributions vary substantially across the five studies, with especially some sites in the BDH and CT-OVC differing quite substantially from those in the other three studies.

We next repeat the analysis for prediction of model shifts in site-specific treatment effects $\tau_g(x)-\tau(x)$, both using IMSE-optimal basis functions and functional PC as predictor (Tables (ref)). These results were obtained combining the data from the 434 treated and 206 control sites, and since the covariance operator of $\mu_g(x;0)$ with the baseline is estimated using only data from the substantially smaller control group, we should expect the resulting estimates to be less precise than for predicting post-intervention average outcomes.

In all five RCTs, treatment was randomized at the cluster level, so we can't directly assess the performance of either type of predictor at that level, but we can still aggregate actual and predicted effects at the level of the study. Here, average predictions based on the IMSE optimal basis functions match the sign and approximate magnitudes of post-hoc realized effects for all values of $K$, whereas at least 3 or 4 functional principal components appear to be necessary to match at least some qualitative aspects of study-level averages. We also report the estimated scores for predicting conditional ATEs plotted in Figure (ref).

For any of these comparisons, it should also be noted that both types of predictions are based on the unconfounded location assumption (Assumption (ref)), whereas realized conditional effects also reflect systematic differences between studies that can't be predicted by extrapolating intra-study variation among sites. Most importantly, the five studies differ in terms of the exact implementation of the incentive, and also country specific factors. Most importantly, cash transfer for CT-OVC in Kenya and BDH in Ecuador were unconditional, whereas transfers under Progresa, PKH, and Tayssir were conditioned on, or connected to, the child's enrollment in school. While cross-site average treatment effects for those two studies were indeed substantially lower than the cross-study average (first column in Table (ref)), our method appears to replicate most of that difference for the Kenyan sites, whereas it fails to reproduce the deviation from the cross-study average only for Ecuador.

Conclusion

We investigate how to exploit observed between-site variation within one or several studies to predict outcomes using baseline data for new “target" sites. The premise of our approach is that agent responses at the micro level follow some universal patterns across study populations. These responses are generally confounded by site-specific factors of an unknown structure, but cross-sectional patterns of attributes and outcome at baseline for each site typically contain useful information regarding those environmental factors in a target site, and may help identify “comparable" sites in the experimental sample. We chose to focus on a nonparametric, linear version of the problem primarily for clarity and ease of implementation, and nonseparable or structural models with sufficiently flexible specifications of site-specific heterogeneity may be another fruitful approach to this problem.

We give a finite-population formulation for the statistical problem of evaluating out of sample forecast performance. We define the target for the transfer estimate as a pseudo-true parameter which reflects the relevant information regarding likely outcomes at the target site that may be learned from previously observed contexts. The corresponding prediction problem is equivalent to functional regression, but given the limited number of sites can only estimate heavily regularized version of the problem. We therefore choose a regularization approach that targets a small number of “most predictive" features of the distribution of outcomes in the baseline. Those optimal predictors are solutions to a generalized eigenvalue problem in terms of the covariance operators of $\mu_g$ and $\tau_g$. The approach can be adapted to sparsely or densely sampled sites, as well as randomization within or between clusters, resulting in different convergence rates.