# Transfer Estimates for Causal Effects across Heterogeneous Sites

Konrad Menzel
New York University

(Date: April 2023. Preliminary - any comments and suggestions will be
appreciated. The author thanks Michael Gechter for sharing data used for this
paper.)

###### Abstract.

We consider the problem of extrapolating treatment effects across
heterogeneous populations (“sites”/“contexts”). We consider an idealized
scenario in which the researcher observes cross-sectional data for a large
number of units across several “experimental” sites in which an intervention
has already been implemented to a new “target” site for which a baseline
survey of unit-specific, pre-treatment outcomes and relevant attributes is
available. We propose a transfer estimator that exploits cross-sectional
variation between individuals and sites to predict treatment outcomes using
baseline outcome data for the target location. We consider the problem of
obtaining a predictor of conditional average treatment effects at the target
site that is MSE optimal within a certain class and subject to data
constraints. Our approach is design-based in the sense that the performance of
the predictor is evaluated given the specific, finite selection of
experimental and target sites. Our approach is nonparametric, and our formal
results concern the construction of an optimal basis of predictors as well as
convergence rates for the estimated conditional average treatment effect
relative to the constrained-optimal population predictor for the target site.
We illustrate our approach using a combined data set of five multi-site
randomized controlled trials (RCTs) to evaluate the effect of conditional cash
transfers on school attendance.

JEL Classification: C14, C21, C53
Keywords: Treatment Effects, Transfer Learning, Nonparametric Estimation,
Functional Regression, Design-Based Estimation

##  1\. Introduction

When scaling up an intervention, it is often necessary to extrapolate
experimental evidence to new sites or contexts even though average causal
effects are likely heterogeneous across sites. We consider a problem in which
cross-sectional information on outcomes and covariates are available for both
experimental and target sites, and we formalize a process of predicting a
causal response that uses pre-intervention outcome data from the target
location to obtain a forecast that is responsive to site-specific
heterogeneity.

In observational studies or randomized controlled trials (RCTs) a researcher
may on occasion observe populations at several distinct sites, either as part
of a larger multi-site study, or by aggregating information from various
studies that were conducted separately in different locations. 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.111Allcott (2015) considers a setting in which a policy
was initially evaluated at 10 sites and eventually scaled up to 111 separate
sites. Dehejia, Pop-Eleches, and Samii (2021) 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 Todd and Wolpin (2006)). Meager (2022) aggregates across seven
different RCTs for Microcredit interventions published in 2015. We therefore
do not aim to provide consistent estimators of treatment parameters for a
theoretical superpopulation. Instead we view the researcher’s task as that of
providing the most accurate feasible prediction given the available
information, and we therefore pursue a more limited goal of characterizing
transfer estimates as pseudo-true parameters that are specific to the finitely
many observed experimental and target sites. That is, we try to find a
procedure that makes efficient use of the data given the specific, finite
population of selected sites.

We cast this problem as a design-based approach 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 an
randomization protocol that assigns experimental versus target status at
random among those clusters. We do necessarily not regard this assignment
mechanism as factually accurate - e.g. the observed assignment may likely
exhibit site-selection effects of the kind documented by Allcott (2015).
Rather, this “as if” data generating process serves a device to define a
pseudo-true treatment parameter that incorporates information on average
effects and between cluster heterogeneity but does not account for site
selection bias. A transfer estimate of this kind should therefore always be
interpreted in light of other sources of extrapolation error that cannot be
addressed without additional data or auxiliary assumptions.

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 656 sites, combining data from five
studies in Mexico, Morocco, Indonesia, Kenya, and Ecuador to illustrate our
approach.

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. 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.

###  1.1. Literature

A conceptual framework for the problem of extrapolation of estimated treatment
effects across heterogeneous sites was developed in the seminal article
byHotz, Imbens, and Mortimer (2005). Using their terminology, we assume
unconfounded locations, but specifically want to allow for (site-specific)
“macro effects,” that is shared heterogeneity in potential outcomes and
treatment effects within each cluster. 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.

Using the same conceptual framework Dehejia, Pop-Eleches, and Samii (2021)
considered the problem of predicting treatment effects at target sites based
on site-specific covariates. Meager (2022) proposes a hierarchical Bayesian
framework for distributional treatment effects. Gechter, Samii, Dehejia, and
Pop-Eleches (2019) use data from two conditional cash transfer programs to
evaluate extrapolation of empirical treatment rules. Adjaho and Christensen
(2022) consider robust extrapolation of treatment rules in the absence of
separate data on the target site. Recent work by Canen and Song (2022) and
Canen and Song (2023) predicts policy effects based on an invariance condition
on an underlying unit-specific response function.

Conceptually, this problem has some resemblance with the method of synthetic
controls (Abadie and Gardeazabal (2003),Abadie, Diamond, and Hainmueller
(2010)). 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. Shin (2022) uses a k-means
algorithm to model unobserved heterogeneity in a problem with cluster
dependence in treatment assignment.

Our approach is 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 Yao, Müller, and Wang (2005a), rather than the dense
case considered by Kraus (2015) and Kneip and Liebl (2020), 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 the classical case 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 a
single 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 (Jiang and Wang (2010),Liebl (2019)) 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 Cai and Hall (2006) and Hall and Horowitz
(2007).

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
Hotelling (1936)’s classical problem of canonical correlation analysis. For
functional data, functional canonical regression has first been proposed by
He, Müller, Wang, and Yang (2010) 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. Our proposal also differs from the elegant approach towards
functional regression with a single regularization parameter in Benatia,
Carrasco, and Florens (2017). Since we do not anticipate to estimate the
regression kernel consistently given typical data constraints in our setting,
we argue that a regularization scheme in terms of finitely many predictive
features (“basis functions”) may be more transparent regarding the basis of
extrapolation for a heavily regularized predictive model.

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.

The remainder of the paper is organized as follows: we first give a formal
characterization of transfer estimation as a statistical problem. We then
derive the mean-square optimal predictors and propose nonparametric estimators
based on the experimental sample. We then give asymptotic properties of those
estimators, assuming the number of experimental sites grows large. The
approach is then illustrated using an application to predicting the causal
effect of conditional cash transfer programs to new locations.

##  2\. Problem Definition

The population of interest consists of $G$ clusters (“sites”/“contexts”),
where the $g$th cluster consists of $N_{g}$ units. There is a binary policy
variable (“treatment”) $D_{gi}\in\\{0,1\\}$ which is assigned at the level of
the unit $i$ in cluster $g$, where we assume that the outcome of interest is
determined only by the unit’s own treatment status. 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 units for
each cluster.

The object of interest is the conditional average treatment effect (CATE)
among units in the $g$th cluster,

$\tau_{g}(x):=\mathbb{E}[Y_{gi}(1)-Y_{gi}(0)|X_{gi}=x]\equiv\mu_{g}(x;1)-\mu_{g}(x;0)$

where we denote the conditional mean functions of potential outcomes with

$\mu_{g}(x;d):=\mathbb{E}[Y_{gi}(d)|X_{gi}=x]$

In Assumption 2.4 below we impose sufficient conditions to ensure that
$\mu_{g}(x;d)$ and $\tau_{g}(x)$ have bounded first and second moments.

The challenge consists in predicting average counterfactuals in the target
cluster when not a single treated unit in that cluster is observed. The
transfer estimate will use information from other sites in combination with
baseline outcomes for the target clusters, and the prediction error is
evaluated unconditionally with respect to the identity of the target cluster.

We assume that the assignment mechanism for $D_{gi}$ has two stages: in a
first stage, clusters are assigned to the experimental population at random,
where we let $R_{g}\in\\{0,1\\}$ denote an indicator that equals one if the
$g$th cluster is assigned to the experimental group, and zero otherwise. For
the experimental clusters $g$ with $R_{g}=1$ each individual unit $i$ is
assigned a binary treatment status $D_{gi}\in\\{0,1\\}$ independently with
probability $p_{g}(x):=\mathbb{P}(D_{gi}=1|X_{gi}=x)$ where $X_{i}$ are the
observable exogenous attributes of the $i$th unit in cluster $g$. All units in
the remaining “target” clusters $g$ with $R_{g}=0$, are assigned treatment
$D_{gi}=0$. In the remainder of the paper we only consider the case of a
single target cluster, which we denote with $g^{*}$.

###  2.1. Transfer Estimates

We consider a transfer estimate $\hat{\tau}_{g^{*},1,\dots,G}(x)$ for
extrapolating from the sites $\\{1,\dots,G\\}\backslash\\{g^{*}\\}$ to
$g^{*}$. That transfer estimate combines information on covariates and
outcomes from the $G$ experimental and target sites to predict the CATE for
the target site, $g^{*}$. Specifically, the site-specific conditional average
treatment effect

$\tau_{g}(x):=\mu_{g}(x;1)-\mu_{g}(x;0)$

is identified for each experimental site
$g\in\\{1,\dots,G\\}\backslash\\{g^{*}\\}$ under standard conditions on the
assignment mechanism, whereas for the target site, only the conditional
distribution of $Y_{g^{*}i}(0)$ given $X_{g^{*}i}$ is identified from a pre-
intervention baseline sample. While other features of the conditional
distribution may also have predictive power, we restrict our focus to linear
predictors for $\tau_{g^{*}}(x)$ based on the conditional mean function
$\mu_{g^{*}}(x)$. This is a functional regression problem with a functional
response, where our focus will be on low-MSE prediction when $G$ may not be
very large.

######  Example 2.1.

(Extrapolation for Conditional Cash Transfer Program) To illustrate the
potential usefulness of using functional data $\mu_{g}(x;0)$ as a predictor
for $\tau_{g}(x;0)$ consider a hypothetical RCT on the effectiveness of a
conditional cash transfer program to incentivize households to invest in a
child’s educational or health outcome $Y_{gi}$, where we also observe a
nominal measure $X_{gi}$ of household income. If household behavior did not
vary fundamentally across locations, the effect of the same nominal transfer
given $X_{gi}$ may still vary across locations if our measure of income was
only imperfectly adjusted for local purchasing power, e.g. using nominal
exchange rates of PPP at the national level.

Specifically suppose that, unbeknownst to the researcher,
$\mu_{g}(x;d)\equiv\mu(x/\pi_{g};d)$, where $\pi_{g}$ is the local price
index, and the function $\mu(x;d)$ is the same for each site $g=1,\dots,G$.
Then if baseline outcomes and/or the response to the cash incentive vary with
household income in local purchasing power, the conditional distribution of
$Y_{gi}(0)$ given measured income $X_{gi}$ reveals information about $\pi_{g}$
and would therefore be predictive of $\tau_{g}(x)$ beyond an adjustment of the
estimated CATE by the average baseline outcome or nominal incomes.

We follow a 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 clusters as fixed. In
particular, the cluster-specific conditional average treatment effects
$\tau_{1}(x),\dots,\tau_{G}(x)$ are treated as deterministic, however the
assignment $R_{g}$ of sites to the experimental population as well as the
selection $D_{gi}$ of treated units within each experimental cluster is
random. In particular, the cross-cluster 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 only
consider the case in which there is a single target cluster in addition to
$G-1$ experimental clusters.

The statistical properties of the estimator are evaluated in terms of the
distribution resulting from the above sampling protocol. For simplicity we
will focus on optimality regarding the integrated mean-squared error (IMSE),

$IMSE:=\frac{1}{G}\sum_{g=1}^{G}\mathbb{E}\left[(1-R_{g})\int(\hat{\tau}_{g,1\dots,G}(x)-\tau_{g}(x))^{2}f_{0}(x)dx\right]$
(2.1)

with a weight function $f_{0}(x)$ that has the properties of a p.d.f. and is
chosen by the researcher. The finite population approach is used as a way of
formalizing the researcher’s problem who aims to produce a forecast based on
data from a limited number of sites. The best feasible prediction under those
circumstances is a parameter that is specific to the set of observable
experimental and target sites.

We show how to use baseline information on observed potential outcomes
$Y_{gi}(0)$ in the target clusters to improve over an unconditional forecast
that doesn’t account for site heterogeneity. Specifically, we consider an
expansion

$\tau_{g}(x)\equiv\tau(x)+\sum_{k=1}^{\infty}t_{gk}\psi_{k}(x)$ (2.2)

where

$\tau(x):=\frac{1}{G}\sum_{g=1}^{G}\tau_{g}(x)$

denotes the expectation of the CATE across clusters, and
$\left\\{\psi_{k}\right\\}_{k\geq 1}$ are a sequence of basis functions to be
discussed in greater detail below. Classical functional principal components
analysis (FPCA) usually works with the Karhunen-Loève representation of a
random function $\tau_{g}(x)$ in terms of the eigenfunctions of the covariance
operator of $\tau_{g}(x)$. However we will choose a different sequence of
basis functions that is then shown to be optimal for this prediction problem
in a sense to be made more precise below.

Those optimal basis functions depend on the unknown covariance operator of the
trajectories $\mu_{g}(x;d)$ across clusters and therefore need to be
estimated. The loadings $t_{g^{*}k}$ for the target cluster are then recovered
based on an estimator for the trajectory $\mu_{g^{*}}(x;0)$. Given the
estimated basis functions $\left\\{\hat{\psi}_{k}\right\\}_{k\geq 1}$ and
loadings $\hat{t}_{g^{*}k}$, the predicted conditional treatment effect is

$\hat{\tau}_{g^{*}}(x):=\hat{\tau}(x)+\sum_{k=1}^{K}\hat{t}_{g^{*}k}\hat{\psi}_{k}(x)$
(2.3)

where we truncate the expansion at an order $K$.

For the scenarios we are envisioning here, the number of experimental clusters
$\sum_{g=1}^{G}R_{g}$ is not very large, so $K$ should be thought of as fairly
small. So rather than aiming for consistent estimation of $\tau_{g}(x)$, we
view the use of the first few leading factors in the expansion (2.2) as a
method of improving over an unconditional forecast that doesn’t account for
site heterogeneity. We could alternatively describe our approach as consistent
estimation of a pseudo-true parameter

$\tau_{g^{*}}^{K}(x):=\tau(x)+\sum_{k=1}^{K}t_{g^{*}k}\psi_{k}(x)$ (2.4)

which approximates the CATE at cluster $g$ at a low order. 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, and may be
preferable to spectral cutoff in some cases (see Carrasco, Florens, and
Renault (2007)), an advantage of our approach 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 trajectory $\mu_{g^{*}}(x;0)$.

###  2.2. Assumptions

We next formalize the identifying conditions which are adapted from Hotz,
Imbens, and Mortimer (2005). 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.

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}$,

######  Assumption 2.1.

(Unconfounded Assignment) For all $g$ with $R_{g}=1$,

$D_{gi}\bot\bot(Y_{gi}(0),Y_{gi}(1))|X_{gi},R_{g}=1$

where $D_{g1},\dots,D_{gN_{g}}$ are also conditionally independent across
units and clusters given attributes and $R_{1},\dots,R_{G}$.

This condition can be satisfied if a randomized controlled trial (RCTs) is
conducted at each location, unconfoundedness is also a standard assumption for
observational studies. Assumption 2.1 therefore 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 A 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:

######  Assumption 2.2.

(Unconfounded Location) $g^{*}$ is drawn uniformly at random from
$\\{1,\dots,G\\}$ independently of
$\left\\{Y_{gi}(0),Y_{gi}(1),X_{gi}:g=1,\dots G,i=1,\dots,N_{g}\right\\}$.

This assumption is an adaption of Assumption 2 in Hotz, Imbens, and Mortimer
(2005) and describes an idealized observational protocol that rules out
systematic ex-ante site selection bias. Our approach specifically aims to
predict ex-post “macro” effects using baseline outcome data after sites have
been selected. In many applications, especially those involving randomized
field trials, there are good reasons to expect this assumption to be
misspecified. In the absence of a systematic approach towards correcting for
site selection bias, our preferred interpretation of the extrapolation
estimate proposed here is therefore that of a pseudo-true parameter, which
aggregates estimates from experimental sites into a constrained “best”
prediction for the target population.

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$.

######  Assumption 2.3.

(Support Conditions) There exists $\delta$, $0<\delta<1$ such that

$\delta<p_{g}(x)<1-\delta\hskip 8.5359pt\textnormal{ and
}\delta<f_{g}(x)/f_{0}(x)<1-\delta$

for all $g\in\\{1,\dots,G\\}\backslash\\{g^{*}\\}$ and $x$ in the support of
$f(x)$.

While Assumption 2.2 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 Hotz, Imbens, and Mortimer (2005). Our
asymptotic rate results imply that nonparametric estimators tend to perform
substantially worse if covariate distributions vary widely across sites.

This assumption also does not cover site-specific aggregate covariates that
may serve as additional predictors as analyzed by Hotz, Imbens, and Mortimer
(2005) and Dehejia, Pop-Eleches, and Samii (2021). 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 3.4.

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:

######  Assumption 2.4.

(Distribution and Moments) For $g=1,\dots,G$, (a) $X_{gi}$ is continuously
distributed with marginal p.d.f. $\sup_{x}f_{g}(x)\leq B_{0}<\infty$. (b) The
site-specific density $f_{g}(x)$, propensity score $p_{g}(x)$, and conditional
mean functions $\mu_{g}(x;d)$ are twice continuously differentiable in $x$
with uniformly bounded derivatives, and (c) the log-derivative $\nabla_{x}\log
f_{g}(x)$ is bounded uniformly over $g=1,\dots,G$. We also assume that (d)
potential outcomes have bounded moments $\mathbb{E}|Y_{gi}(d)|^{s}<\infty$ for
$d=0,1$ and some $s>3$, (e) there exist finite constants $B_{0},B_{1}$ such
that $\sup_{x}f(x)\leq B_{0}$ and
$\sup_{x}\mathbb{E}[|Y_{gi}(d)|^{s}|X_{gi}=x]f_{g}(x)\leq B_{1}$ for all $g$.

With the exception of part (c), these conditions are commonly assumed for
nonparametric estimation of conditional moments, see e.g. Hansen (2008), and
part (c) could be omitted if the Nadaraya-Watson (kernel regression) estimator
was used instead of the local linear estimator. Another feature to note is
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$. 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 Yao,
Müller, and Wang (2005a)), however we do not find such an assumption plausible
for the problem considered here.

##  3\. Estimation Approach

The best (lowest variance) predictor for the treatment effect in the target
cluster $g^{*}$ given $\mu_{g^{*}}(x;0)$ is given by the conditional
expectation function,

$\mathbb{E}[\tau_{g^{*}}^{*}|\mu_{g^{*}}(:,0)]=\mathbb{E}\left[\mu_{g^{*}}(X_{gi};1)-\mu_{g^{*}}(X_{gi};0)|\mu_{g^{*}}(\cdot;0)\right]$

As a practical matter, this predictor is infeasible because the site-specific
conditional mean function $\mu_{g^{*}}(:;0)$ has to be estimated, and the
dependence structure between $\mu_{g}(\cdot;1)$ and $\mu_{g}(\cdot;0)$ has to
be inferred from a fairly small number of sites.

Our more limited objective is therefore to form a predictor that makes
efficient use of information on $\mu_{g^{*}}(\cdot;0)$ given the practical
limitations given the available data. In order to operationalize this idea, we
restrict our attention to predictors that are linear in $\mu_{g}(x)$, i.e.
that can be expressed in terms of a linear operator $B$ acting on $\mu_{g}$,

$\tilde{\tau}_{g}(x):=\int\mu_{g}(x_{1};0)\beta(x_{1},x)f_{0}(x_{1})dx_{1}=:(B\mu_{g})(x)$
(3.1)

where $\beta(x_{1},x_{2})$ is a square integrable function. That is, we can
view the predictor as a Hilbert-Schmidt integral operator acting on
$\mu_{g}(0)$ with kernel function $\beta(x_{1},x_{2})$. We then estimate a
finite-dimensional approximation to $\beta(x_{1},x_{2})$ in terms of $K$ basis
functions $\phi_{k}(x_{1}),\psi_{k}(x_{2})$. One of our main results
establishes the optimal choice of those basis functions at any order of
approximation.

Our general approach is based on a representation of the random processes
$\mu_{g}(x;0)$ and $\tau_{g}(x)$ in terms of orthogonal bases,
$\phi_{1},\phi_{2},\dots$ and $\psi_{1},\psi_{2},\dots$,

$\displaystyle\mu_{g}(x;0)$ $\displaystyle=$
$\displaystyle\mu(x;0)+\sum_{k=1}^{\infty}m_{gk}\phi_{k}(x)$
$\displaystyle\tau_{g}(x)$ $\displaystyle=$
$\displaystyle\tau(x)+\sum_{k=1}^{\infty}t_{gk}\psi_{k}(x)$

A standard representation of the random processes $\mu_{g}(x;0)$ and
$\tau_{g}(x)$ is the Karhunen-Loève 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)$. 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 approximation 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)$.

We then use an appropriate (small) number $K$ of basis functions for
prediction, where $K<<G$. This approach requires nonparametric estimation of
the mean functions

$\mu(x;d):=\mathbb{E}[\mu_{g^{*}}(x;d)]\equiv\frac{1}{G}\sum_{g=1}^{G}\mu_{g}(x;d),\hskip
14.22636ptd=0,1$

and

$\tau(x;d):=\mathbb{E}[\mu_{g^{*}}(x;d)]\equiv\frac{1}{G}\sum_{g=1}^{G}\tau_{g}(x;d),\hskip
14.22636ptd=0,1$

as well as the covariance kernels

$\displaystyle H_{\mu\mu}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\textnormal{Cov}(\mu_{g}(x_{1};0)-\mu(x_{1};0),\mu_{g}(x_{2};0)-\mu(x_{2};0))$
$\displaystyle=$
$\displaystyle\frac{1}{G}\sum_{g=1}^{G}(\mu_{g}(x_{1};0)-\mu(x_{1};0))(\mu_{g}(x_{2};0)-\mu(x_{2};0))$
$\displaystyle H_{\mu\tau}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\textnormal{Cov}(\mu_{g}(x_{1};0)-\mu(x_{1};0),\tau_{g}(x_{2})-\tau(x_{2}))$
$\displaystyle=$
$\displaystyle\frac{1}{G}\sum_{g=1}^{G}(\mu_{g}(x_{1};0)-\mu(x_{1};0))(\tau_{g}(x_{2})-\tau(x_{2}))$
$\displaystyle H_{\tau\tau}(x_{1},x_{2})$ $\displaystyle:=$
$\displaystyle\textnormal{Cov}(\tau_{g}(x_{1};0)-\tau(x_{1};0),\tau_{g}(x_{2})-\tau(x_{2}))$
$\displaystyle=$
$\displaystyle\frac{1}{G}\sum_{g=1}^{G}(\tau_{g}(x_{1};0)-\tau(x_{1};0))(\tau_{g}(x_{2})-\tau(x_{2}))$

where covariances are with respect to a random draw of a site $g$ from the
discrete uniform distribution over $\\{1,\dots,G\\}$.

###  3.1. Optimal Basis Functions

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 (2.1).

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

$\displaystyle(T_{\mu\mu}\varphi)(x)$ $\displaystyle:=$ $\displaystyle\int
H_{\mu\mu}(x_{1},x)\varphi(x_{1})f_{0}(x_{1})dx_{1}$
$\displaystyle(T_{\mu\tau}\varphi)(x)$ $\displaystyle:=$ $\displaystyle\int
H_{\mu\tau}(x_{1},x)\varphi(x_{1})f_{0}(x_{1})dx_{1}$
$\displaystyle(T_{\tau\tau}\varphi)(x)$ $\displaystyle:=$ $\displaystyle\int
H_{\tau\tau}(x_{1},x)\varphi(x_{1})f_{0}(x_{1})dx_{1}$

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
$\phi_{1}(x),\dots,\phi_{K}(x)$ of $K$ functions 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}^{K}a_{k}\phi_{k}:a_{1},\dots,a_{K}\in\mathbb{R}\right\\}$

By the classical projection theorem (Theorem 2 on p.51 in Luenberger (1969))
that projection is well-defined. We also let

$\displaystyle IMSE_{K}$ $\displaystyle\equiv$ $\displaystyle
IMSE_{K}[\phi_{1},\dots,\phi_{K}]$ $\displaystyle:=$
$\displaystyle\int\min_{B\in\mathcal{H}_{K}\times
L_{2}(\mathcal{X})}\mathbb{E}\left[(\tau_{g}(x)-BP_{K}\mu_{g}(x))^{2}\right]f_{0}(x)dx$

denote the integrated mean-square error of prediction 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\ker(T_{\mu\mu})$, 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}$

######  Lemma 3.1.

Suppose that $T_{\mu\mu},T_{\mu\tau}$ are compact operators, and define

$IMSE_{K}^{*}:=\inf_{\phi_{1},\dots,\phi_{K}}\sum_{k=1}^{K}\langle\phi_{k},T_{\mu\tau}T_{\mu\tau}^{*}\phi_{k}\rangle$
(3.2)

where the infimum is taken over $\phi_{1},\dots,\phi_{K}\in
L_{2}(\mathcal{X})$ such that
$\langle\phi_{k},T_{\mu\mu}\phi_{l}\rangle=\delta_{kl}$ for all
$k,l=1,\dots,K$. Then for an arbitrary choice of
$\phi_{1}(x),\dots,\phi_{K}(x)$ we have that $IMSE_{K}\geq IMSE_{K}^{*}$.

See the appendix for a proof. 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}$. This problem bears
some resemblance with functional canonical analysis, where He, Müller, and
Wang (2003) propose high-level conditions on the cross-correlation operator
which would also be sufficient to guarantee that the infimum in (3.2) 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 a regularized version of the generalized eigenvalue
problem

$T_{\mu\tau}T_{\mu\tau}^{*}\phi_{k}^{*}=\lambda_{k}T_{\mu\mu}\phi_{k}^{*}\hskip
14.22636pt\textnormal{for each }k=1,\dots,K$ (3.3)

where we select the eigenfunctions $\phi_{1}^{*},\dots,\phi_{K}^{*}$
associated with the $K$ leading eigenvalues
$|\lambda_{1}|\geq|\lambda_{2}|\geq\dots$. 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.

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. Specifically, we consider the
following generalized eigenvalue problem

$T_{\mu\tau}T_{\mu\tau}^{*}\phi_{ka}^{*}=\lambda_{ka}(T_{\mu\mu}+a\textnormal{Id})\phi_{ka}^{*}\hskip
14.22636pt\textnormal{for each }k=1,\dots,K$ (3.4)

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\geq 1$, 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 (3.4) with

$IMSE_{K}^{*}(a):=\int\min_{B\in\mathcal{H}_{K}^{*}\times
L_{2}(\mathcal{X})}\mathbb{E}\left[(\tau_{g}(x)-BP_{K}^{*}\mu_{g}(x))^{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 (3.4) corresponding to the $K$ largest
eigenvalues are approximately optimal as $a\rightarrow 0$:

######  Theorem 3.1.

(Optimal Basis for $\mu_{g}(x;0)$) Suppose that $T_{\mu\mu}$ and $T_{\mu\tau}$
are compact operators. Then for any $a>0$ and fixed $K$ there exists a
solution $\phi_{1a}^{*},\dots,\phi_{Ka}^{*}$ the functions solving the
generalized eigenvalue problem (3.4), and the resulting IMSE satisfies

$IMSE_{K}^{*}(a)\leq IMSE_{K}^{*}+o(1)$

as $a\rightarrow 0$.

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.

We next give results on the behavior of $IMSE_{K}(a)$ under conditions for
which a solution to the non-regularized problem (3.3) exists. Following He,
Müller, Wang, and Yang (2010), we state those sufficient conditions in terms
of separate Karhunen-Loève expansions for $\mu_{g}(x;0)$ and $\tau_{g}(x;0)$,

$\displaystyle\mu_{g}(x;0)$ $\displaystyle=$
$\displaystyle\mu(x;0)+\sum_{k=1}^{\infty}\alpha_{gk}\xi_{k}(x)$
$\displaystyle\tau_{g}(x)$ $\displaystyle=$
$\displaystyle\tau(x)+\sum_{k=1}^{\infty}\beta_{gk}\zeta_{k}(x)$ (3.5)

where $\xi_{1},\xi_{2},\dots$ and $\zeta_{1},\zeta_{2},\dots$ are orthonormal
systems corresponding to the eigenfunctions of the covariance operators for
$\mu_{g}(\cdot;0)$ and $\tau_{g}(\cdot)$ associated with eigenvalues
$\mathbb{E}[\alpha_{gk}^{2}]$ and $\mathbb{E}[\beta_{gk}^{2}]$, respectively.
By the usual properties of the Karhunen-Loève representation, the coefficients
satisfy $\mathbb{E}[\alpha_{gk}]=\mathbb{E}[\beta_{gk}]=0$ and
$\mathbb{E}[\alpha_{gk}\alpha_{gl}]=\mathbb{E}[\beta_{gk}\beta_{gl}]=0$ for
all $k=1,2,\dots$, and $l\neq k$.

We then analyze regularization bias under the following condition:

$\sum_{k=1}^{\infty}\sum_{l=1}^{\infty}\frac{\mathbb{E}[\alpha_{gk}\beta_{gl}]^{2}}{\mathbb{E}[\alpha_{gk}^{2}]^{3/2}}<\infty$
(3.6)

Note that by compactness of $T_{\mu\mu}$,
$\mathbb{E}[\alpha_{gk}^{2}]\rightarrow 0$ as $k$ grows, so that this
condition is slightly weaker than the sufficient condition given by He,
Müller, Wang, and Yang (2010) for the existence of a functional linear
regression model representing the linear projection of $\tau_{g}$ on $\mu_{g}$
(see their Proposition 2.4).

While the squared correlation coefficient
$r_{kl}^{2}:=\frac{\mathbb{E}[\alpha_{gk}\beta_{gl}]^{2}}{\mathbb{E}[\alpha_{gk}^{2}]\mathbb{E}[\beta_{gl}^{2}]}$
is bounded by 1 for any fixed $k,l$ and $\mathbb{E}[\beta_{gl}^{2}]$ are
square summable in $l$, (3.6) may fail when there are many features of the
conditional mean of baseline outcomes $Y_{gi}(0)$ given $X_{gi}$ that have low
variability in $\alpha_{gk}$ but are highly predictive of $\tau_{g}$.

######  Proposition 3.1.

Suppose that $\mu_{g}(x;0)$ and $\tau_{g}(x)$ have a Karhunen-Loève
representation (3.5) with coefficients satisfying (3.6). Then the solutions to
the problem (3.3) are well-defined, and the conclusion of Theorem 3.1 holds
with

$IMSE_{K}^{*}(a)\leq IMSE_{K}^{*}+O(a)$

See the appendix for a proof. This result establishes not only that the
optimal IMSE can be achieved at a linear rate in $a$, but also that the
condition (3.6) is sufficient to guarantee that the optimal predictors defined
in (3.3) exist independently of a particular regularization scheme, and can
also be estimated consistently under conditions given in the next section. We
are not aware of any clear correspondence of the summability condition (3.6)
to conventional assumptions on the potential outcomes model, and we therefore
also want to expressly allow for the case in which it fails.

Given the proposed choice of $\phi_{1}^{*},\dots,\phi_{K}^{*}$, we also state
the projection of $\tau_{g}$ onto the optimal basis:

######  Corollary 3.1.

For any $K=1,\dots,K^{*}$, the projection based on the solution of (3.4) takes
the form

$\tau_{gK}^{*}(x):=\sum_{k=1}^{K}t_{gk}\psi_{ka}^{*}(x)$

where $t_{gk}:=\frac{1+a}{1-a}\langle\mu_{g},\phi_{ka}^{*}\rangle$ and

$\psi_{ka}^{*}(x):=\left(T_{\mu\tau}^{*}\phi_{ka}^{*}\right)(x)$ (3.7)

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_{g1},\dots,t_{gK})$ 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}$:

$\displaystyle T_{\mu\tau}^{*}T_{\mu\mu a}^{-1}T_{\mu\tau}\psi_{ka}^{*}$
$\displaystyle=$ $\displaystyle T_{\mu\tau}^{*}T_{\mu\mu
a}^{-1}T_{\mu\tau}T_{\mu\tau}^{*}\phi_{ka}^{*}$ $\displaystyle=$
$\displaystyle\lambda_{k}T_{\mu\tau}^{*}\phi_{ka}^{*}=\lambda_{k}\psi_{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 (3.4) 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})$.

###  3.2. Related Approaches

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

The procedure we arrive at starting from optimality considerations has a close
resemblance to canonical correlation analysis which has previously been
proposed for functional regression problems by He, Müller, Wang, and Yang
(2010). 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 (3.3),
He, Müller, and Wang (2003) and He, Müller, Wang, and Yang (2010) impose high-
level conditions on the cross-correlation operator to ensure existence (see
Proposition 4.2 in He, Müller, and Wang (2003)). 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 Cupidon, Eubank, Gilliam, and
Ruymgaart (2008) 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.

Yang, Müller, and Stadtmüller (2011) 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)$

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 Benatia, Carrasco, and Florens (2017)
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)$.

###  3.3. Estimation

The representation in Corollary 3.1 motivates an estimator of the form

$\hat{\tau}_{g}(x):=\sum_{k=1}^{K}\hat{\tilde{m}}_{gk}\hat{\psi}_{k}(x)$

where $\hat{\tilde{m}}_{gk}=\langle\hat{\mu}_{g},\hat{\phi}_{k}\rangle$ for a
nonparametric estimator $\hat{\mu}_{g}$ of $\mu_{g}(x;0)$ and basis functions
$\hat{\phi}_{1},\dots,\hat{\phi}_{K}$ solving the generalized eigenvalue
problem (3.3) where $T_{\mu\mu}$ and $T_{\mu\mu}^{-1/2}T_{\mu\tau}$ are
replaced with estimates.

Here we develop our approach for the case of densely sampled clusters,
$n\rightarrow\infty$, for a discussion of the setting with sparse samples, see
Appendix A. 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 Yao, Müller, and Wang
(2005a) we give results for local linear estimators: For each experimental
cluster, let

$\hat{\mu}_{g}(x;d):=\arg_{b_{0}}\min_{b_{0},b_{1}}\sum_{i=1}^{n_{g}}w_{gi}(x;d)(Y_{gi}-b_{0}-b_{1}(x-X_{gi}))^{2}$
(3.8)

where we let the nonparametric weights

$w_{gi}(x;d):=1\textnormal{l}\\{D_{gi}=d\\}K\left(\frac{X_{gi}-x}{h}\right).$

Here, $K(u)$ is a kernel function, 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

$\hat{M}_{g}(x_{1},x_{2};d_{1},d_{2}):=\arg_{b_{0}^{(g)}}\min_{b_{0}^{(g)},b_{11}^{(g)},b_{12}^{(g)}}\sum_{j\neq
i}H_{gij}(x_{1},x_{2},\mathbf{b})w_{gi}(x_{1};d)w_{gj}(x_{2};d)$ (3.9)

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

$\displaystyle\hat{\mu}(x;d)$ $\displaystyle:=$
$\displaystyle\frac{1}{G-1}\sum_{g=1}^{G}R_{g}\hat{\mu}_{g}(x;d)$
$\displaystyle\hat{H}(x_{1},x_{2};d_{1},d_{2})$ $\displaystyle:=$
$\displaystyle\frac{1}{G-1}\sum_{g=1}^{G}R_{g}\hat{M}_{g}(x_{1},x_{2};d_{1},d_{2})-\hat{\mu}(x_{1};d_{1})\hat{\mu}(x_{2};d_{2})$

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 Section 4, the optimal rate turns out to be the same for either estimator.
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 Yao,
Müller, and Wang (2005a) 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{\tau}(x):=\hat{\mu}(x;1)-\hat{\mu}(x;0)$

and

$\hat{H}_{\mu\mu}(x_{1},x_{2}):=\hat{H}(x_{1},x_{2};0,0)\hskip
14.22636pt\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

$\displaystyle(\hat{T}_{\mu\mu}h)(x)$ $\displaystyle=$
$\displaystyle\int\hat{H}_{\mu\mu}(x,s)h(x)f_{0}(s)ds$
$\displaystyle(\hat{T}_{\mu\tau}h)(x)$ $\displaystyle=$
$\displaystyle\int\hat{H}_{\mu\tau}(x,s)h(x)f_{0}(s)ds$
$\displaystyle(\hat{T}_{\mu\tau}^{*}h)(x)$ $\displaystyle=$
$\displaystyle\int\hat{H}_{\mu\tau}(s,x)h(x)f_{0}(s)ds$

for any square integrable function $h$.

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

$\hat{T}_{\mu\mu
a}^{-1/2}\hat{T}_{\mu\tau}\hat{T}_{\mu\tau}^{*}\hat{T}_{\mu\mu
a}^{-1/2}\hat{\xi}_{k}=\hat{\lambda}_{k}\hat{\xi}_{k}$ (3.11)

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

$\hat{\phi}_{k}:=\hat{T}_{\mu\mu a}^{-1/2}\hat{\xi}_{k}.$ (3.12)

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

We then construct $\hat{\psi}_{k}$ by applying the estimator of
$T_{\mu\tau}^{*}$ to the estimated eigenfunction $\hat{\phi}_{k}$,

$\hat{\psi}_{k}(x):=\left(\hat{T}_{\mu\tau}^{*}\hat{\phi}_{k}\right)(x)\equiv\int\hat{H}_{\mu\tau}(s,x)\hat{\phi}_{k}(s)f_{0}(s)ds$

for $k=1,\dots,K$. Using these estimates, we then obtain

$\widehat{\tilde{m}}_{gk}:=\langle\hat{\mu}_{g},\hat{\phi}_{k}\rangle$

Substituting this into the formula from Corollary 3.1, our estimate of the
conditional ATE $\tau_{g}(x)$ is

$\hat{\tau}_{g^{*}}(x)=\sum_{k=1}^{K}\widehat{\tilde{m}}_{g^{*}k}\hat{\psi}_{k}(x)$

###  3.4. 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 $G-1$ of experimental sites
is not very large, the scope for adjustments for $W_{g}$ will be fairly
limited 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 2.1 is changed to

$D_{gi}\bot\bot(Y_{gi}(0),Y_{gi}(1))|X_{gi},W_{g},R_{g}=1$

and Assumption 2.2 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_{gi}(d)|X_{gi}=x,W_{g}=w]$

and covariance function

$\textnormal{
}H(x_{1},x_{2},w;d_{1},d_{2}):=\mathbb{E}\left[\left.\frac{}{}(Y_{g1}(d_{1})-\mu(x_{1},w;d_{1}))(Y_{g2}(d_{2})-\mu(x_{2},w;d_{2}))\right|X_{g1}=x_{1},X_{g2}=x_{2},W_{g}=w\right]$

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_{1}d_{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 potential 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. Dehejia, Pop-Eleches, and Samii (2021)
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.

##  4\. Asymptotic Theory

This section gives convergence rates for the estimated eigenfunctions
$\hat{\phi}_{k},\hat{\psi}_{k}$ and the resulting integrated mean square error
relative to the optimal benchmark $IMSE_{K}^{*}$. We consider rates as both
$G$ and $n_{g}\equiv n$ grow to infinity, possibly at different rates. In this
section we only consider the case of densely sampled clusters,
$n\rightarrow\infty$, estimators and rates for the sparsely sampled case are
discussed in Appendix A.

###  4.1. Nonparametric Estimation of Covariance Functions

We first derive convergence rates for the local linear estimator of the
covariance function $H(x_{1},x_{2};d_{1},d_{2})$. Following Hansen (2008), we
assume the following regarding the kernel function

######  Assumption 4.1.

(Kernel Function) The multivariate kernel function
$K:\mathbb{R}^{d}\rightarrow\mathbb{R}$ satisfies (a)
$|K(u)|\leq\bar{K}<\infty$ and $\int|K(u)|du\leq\mu<\infty$. Furthermore, (b)
$K(u)$ is differentiable and there exist $\Lambda_{1},L<\infty$ and $\nu>1$
such that
$\left\|\nabla_{u}K(u)\right\|\leq\Lambda_{1}\min\left\\{1,\|u\|^{-\nu}\right\\}$
for $\|u\|>L$. (c) The first two moments of the kernel satisfy the conditions
$\int K(u)du=1$, $\int uK(u)du=0$, and $\int uu^{\prime}K(u)du=\Omega$, a
finite, positive definite matrix.

Parts (a) and (c) are fairly standard in the literature, the bound in part (b)
is important for the uniformity arguments. These assumptions are satisfied by
commonly used kernel functions such as the Gaussian or Epanechnikov kernel.
Hansen (2008) and Graham, Niu, and Powell (2021) consider an alternative set
of conditions for kernel functions with bounded support that need not be
differentiable, which could be used to replace part (b).

Following Hansen (2008) we consider uniform convergence over a potentially
expanding subset of the support of $X$, $\\{x:\|c\|\leq c_{n}\\}$. We also
denote

$\delta_{Gn}:=\min_{d\in\\{0,1\\}}\min_{g=1,\dots,G}\inf_{\|x\|\leq
c_{n}}f_{g}(x;\mathbf{d})>0$

where the minimum is attained for each finite $G$. The sequence $\delta_{Gn}$
reflects that the density $f_{g}(x;\mathbf{d})$ is not necessarily bounded
away from zero as we expand the domain on which the conditional mean function
is estimated and add new clusters under the asymptotic sequence for which the
distribution of the covariate $X_{gi}$ may differ substantially.

We can now state our main result regarding the rate of consistency for the
local linear estimator for the conditional mean functions and their covariance
kernel:

######  Theorem 4.1.

Suppose that Assumptions 2.1-2.4 and 4.1 hold. Furthermore, the bandwidth $h$
satisfies $\left(\frac{n}{\log n}\right)^{1/3}h^{d}\rightarrow\infty$ and
$\delta_{Gn}^{-1}\left(\frac{\log n}{Gnh^{d}}\right)^{1/2}\rightarrow 0.$ Then
the local linear estimators for $\mu(x;\mathbf{d})$ and the covariance
operator $H_{d_{1}d_{2}}(x_{1},x_{2})$ are consistent at the rate

$\displaystyle\sup_{d,x_{1}}\left|\hat{\mu}(x_{1},d)-\mu(x_{1},d)\right|$
$\displaystyle=$ $\displaystyle O_{p}\left(r_{Gn}\right)$
$\displaystyle\sup_{d_{1},d_{2},x_{1},x_{2}}\left|\hat{H}(x_{1},x_{2};d_{1},d_{2})-H(x_{1},x_{2};d_{1},d_{2})\right|$
$\displaystyle=$ $\displaystyle O_{p}\left(r_{Gn}\right)$

where $r_{Gn}=\frac{1}{G}+\delta_{Gn}^{-1}\left(h^{2}+\left(\frac{\log
n}{Gnh^{d}}\right)^{1/2}\right)$, and the suprema are for
$\|x_{1}\|,\|x_{2}\|\leq c_{n}$. The rate optimal bandwidth minimizing the
second term of either error is of the order $h^{*}=O\left(\left(\frac{\log
n}{Gn}\right)^{\frac{1}{4+d}}\right)$, resulting in a rate
$\frac{1}{G}+\delta_{Gn}^{-1}\left(\frac{\log n}{Gn}\right)^{\frac{2}{4+d}}$.

The proof uses uniform consistency results by Hansen (2008) and Graham, Niu,
and Powell (2021) for nonparametric regression using cross-sectional and
dyadic data and is provided in the appendix. The convergence rate consists of
two separate contributions: the first component represents the error from
estimating moments from the $G-1$ experimental clusters, which excludes the
single target cluster. The rate for this component is of the order
$O_{p}(G^{-1})$ rather than its square root, owing to the fact that the
estimands are the mean and covariance function for the sample of $G$ sites,
$G-1$ out of which are used for estimation. If instead the population of
interest consisted of a greater number additional target clusters growing at
least at the order of $G$ or an infinite meta-population, the rate on that
leading term would instead be $O(G^{-1/2})$.

The second component represents the sampling error in nonparametric estimation
of conditional mean functions in each cluster, where the optimal bandwidth
sequence balances the respective rates for the smoothing bias and sampling
error. It is also instructive to compare the rate to that in Theorem 1 of Yao,
Müller, and Wang (2005a) who consider the case of sparse (finite-$n$) rather
than dense samples from each site. For one the effective dimension for
nonparametric estimation of the covariance kernel is only $d$ in our problem
rather than $2d$, reflecting the fact that the contribution for each cluster
is a U-statistic averaging over $\binom{n}{2}$ terms, so that allowing $n$ to
grow results in a more favorable rate. If $G$ is small relative to $n$, the
first term may dominate the second contribution of nonparametric estimation
error in which case slower rates due to cross-site variation in the covariate
distribution (as captured by $\delta_{Gn}$) need not be a first-order concern.

Moreover, in the rate for the local linear estimator derived by Hansen (2008),
$\delta^{-1}$ appears as a larger power than for the Nadaraya Watson (local
constant) estimator, suggesting a greater fragility of the local linear
estimator on unbounded domains when the density $f_{g}(x;d)$ can be close to
zero. We impose mild restrictions on the bandwidth in Theorem 4.1, as well as
a bound on the log derivative of $f_{g}(x;d)$ in Assumption 2.4 part (c),
eliminating the source of potential instability of the local linear estimator,
thereby matching the rate for the Nadaraya-Watson estimator.

###  4.2. Estimation of Basis Functions $\phi_{k},\psi_{k}$

We next consider convergence rates for the estimated eigenfunctions given a
preliminary estimator for mean and covariance function for the conditional
average treatment effect function $\mu_{g}(x)$. We distinguish between
settings where the optimal predictors $\phi_{k}^{*}$ are well-defined in the
absence of regularization according to (3.3) and the general case in which we
consider estimation of $\phi_{ka}^{*}$ for the regularized version of the
problem (3.4), and $\psi_{ka}^{*}$ given by Corollary 3.1.

We first state result in terms of a generic nonparametric estimator and its
convergence rate, both under the inner product norm $\|\cdot\|$ as well as the
sup norm. Asymptotic rates based on the local linear estimator are given
further below.

######  Theorem 4.2.

Suppose the estimator $\hat{H}(x_{1},x_{2};d_{1},d_{2})$ is consistent with
rate

$\sup_{x_{1},x_{2},d_{1},d_{2}}\left|\hat{H}(x_{1},x_{2};d_{1},d_{2})-H(x_{1},x_{2};d_{1},d_{2})\right|=O_{p}(r_{Gn}).$

Then the estimators for the generalized eigenvalue problem (3.4) with
regularization parameter $a>0$ are consistent at respective rates
$|\hat{\lambda}_{k}-\lambda_{ka}|=O_{p}\left(a^{-3/2}r_{Gn}\right)$,
$\|\hat{\psi}_{k}-\psi_{ka}\|=O_{p}\left(a^{-3/2}r_{Gn}\right)$,
$\|\hat{\phi}_{k}-\phi_{ka}\|=O_{p}\left(a^{-3/2}r_{Gn}\right)$ for each
$k=1,\dots,K$. Furthermore,
$\sup_{x}|\hat{\psi}_{k}(x)-\psi_{ka}(x)|=O_{p}(a^{-3/2}r_{Gn})$ and
$\sup_{x}|\hat{\phi}_{k}(x)-\phi_{ka}(x)|=O_{p}(a^{-3/2}r_{Gn})$ for each
$k=1,\dots,K$.

The proof has a similar overall structure as that for similar results in
Gobet, Hoffmann, and Reiß (2004) and Yao, Müller, and Wang (2005a)
(Proposition 4.2 and Theorem 2, respectively), but requires some adjustments
for the present site. For the problem considered here, estimation of the
operator itself requires regularization, and furthermore the rank of all
covariance operators is less than or equal to $G$. We also allow the dimension
of $x$ to be greater than one and some of the relevant eigenvalues need not be
unique. While each of these adjustments is incremental and leaves the general
structure of the argument unchanged, we provide a self-contained proof in the
appendix.

For the local linear estimator we can immediately obtain the following from
Theorems 4.1 and 4.2:

######  Corollary 4.1.

For each $k=1,\dots,K$, let $\hat{\psi}_{k},\hat{\phi}_{k}$ and
$\hat{\nu}_{k},\hat{\lambda}_{k}$ the estimators for eigenfunctions and
eigenvalues using local linear estimators for $\mu(x)$ and
$H_{d_{1}d_{2}}(x_{1},x_{2})$ with a bandwidth sequence $h$ satisfying
$\left(\frac{n}{\log n}\right)^{1/3}h^{d}\rightarrow\infty$ and
$\delta_{Gn}^{-1}\left(\frac{\log n}{nh^{d}}\right)^{1/2}.$ Under the
assumptions for Theorem 4.1, the conclusions of Theorem 4.2 hold for these
estimators with the rate

$r_{Gn}=\frac{1}{G}+\delta_{Gn}^{-1}\left(h^{2}+\left(\frac{\log
n}{Gnh^{d}}\right)^{1/2}\right)$

Similarly, we can give the asymptotic rate for estimating the best linear
predictor for the conditional average treatment effect as defined in (2.3),
where the rates of individual components follow immediately from Theorem 4.2.

######  Corollary 4.2.

Suppose the estimators $\hat{\mu}(x;d)$ and
$\hat{H}_{d_{1}d_{2}}(x_{1},x_{2})$ are consistent with rates
$\sup_{x_{1},d_{1}}\left|\hat{\mu}(x_{1};d_{1})-\mu(x_{1};d_{1})\right|=O_{p}(r_{Gn})$
and
$\sup_{x_{1},x_{2},d_{1},d_{2}}\left|\hat{H}(x_{1},x_{2};d_{1},d_{2})-H(x_{1},x_{2};d_{1},d_{2})\right|=O_{p}(r_{Gn})$.
Then for any fixed value of $K$, the estimator $\hat{\tau}_{g}^{K}(x)$ based
on (2.3) is consistent at the rate

$\sup_{x\in\mathcal{X}}\left|\hat{\tau}_{g}(x)-\tau_{g}^{K}(x)\right|=O_{p}\left(a^{-3/2}r_{Gn}\right)$

At present we do not derive an asymptotic distribution for functionals of
$\hat{\tau}_{g}(x)$. Asymptotic normality of certain functionals of estimated
eigenfunctions were derived in a different site by Christensen (2014). Whether
his strategy of proof can be adapted to derive distributional results for our
setup is left for future research.

We finally consider the case in which condition (3.6) holds and the solution
to the unregularized problem (3.3) is well defined. It can be seen from the
proof of Theorem 4.2 that regularization is generally still necessary for
estimation of the optimal functions $\phi_{1}^{*},\dots,\phi_{K}^{*}$, but we
can give a rate with which the IMSE given the estimated functions converges to
the lower bound in Lemma 3.2. Specifically, for any $K$ square-integrable
functions $\phi_{1},\dots,\phi_{K}$, we define

$IMSE_{K}[\phi_{1},\dots,\phi_{K}]:=\int\min_{B\in\hat{\mathcal{H}}_{K}\times
L_{2}(\mathcal{X})}\mathbb{E}\left[(\tau_{g}(x)-BP_{K}\mu_{g}(x))^{2}\right]f_{0}(x)dx$

where $P_{K}:L_{2}(\mathcal{X})\rightarrow\mathcal{H}_{K}$ denotes the
operator associated with orthogonal projection onto the closed linear subspace
$\mathcal{H}_{K}:=\textnormal{span}\left(\phi_{1},\dots,\phi_{K}\right)$.
Combining the rate results in Proposition 3.1 and Theorem 4.2, we can then
give the following

######  Corollary 4.3.

Suppose that the Assumptions of Theorems 4.1 and 4.2 hold along with Condition
(3.6). For the local linear estimator with bandwidth
$h^{*}=O\left(\left(\frac{\log n}{Gn}\right)^{\frac{1}{4+d}}\right)$ and
regularization parameter $a>0$,

$|IMSE_{K}[\hat{\phi}_{1},\dots,\hat{\phi}_{K}]-IMSE_{K}^{*}|=O_{P}\left(a+a^{-3/2}r_{Gn}\right)$

where $r_{Gn}:=\frac{1}{G}+\delta_{Gn}^{-1}\left(\frac{\log
n}{Gn}\right)^{\frac{2}{4+d}}$.

The rate for the IMSE in $a$ and $r_{Gn}$ suggests that this upper bound is
minimized at a rate $a=r_{nG}^{2/5}$ for the regularization parameter. For
methods of functional data analysis, it is common to determine regularization
parameters by cross-validation (see e.g. He, Müller, Wang, and Yang (2010)),
however we do not formally establish consistency of such a method.

##  5\. Empirical Application [preliminary]

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
(Schultz (2004) and Todd and Wolpin (2006)), and similar programs have
subsequently been implemented in over 50 other countries (see Banerjee, Hanna,
Kreindler, and Olken (2017) for a recent summary).

###  5.1. Data

We combine samples from _PROGRESA_ with four additional randomized studies
that were conducted in Indonesia (_Program Keluarga Harapan (PKH)_ , see
Alatas (2011) and Cahyadi, Hanna, Olken, Prima, Satriawan, and Syamsulhakim
(2020)), Morrocco (_Tayssir_ , see Benhassine, Devoto, Duflo, Dupas, and
Pouliquen (2015)), Kenya (_Kenya CT-OVC_ , Team (2012)), and Ecuador (_Bono de
Desarrollo Humano (BDH)_ , Edmonds and Schady (2012)).333These 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_ , Barrera-Osorio, Bertrand, Linden, and Perez-Calle
(2011)) 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.

Our main focus is on leveraging unobserved heterogeneity across locations for
cross-study comparisons. 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, 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 try to assess to what extent between-study variation can be predicted
from between cluster variation on a more disaggregated level.

We retain all observations of eligible households 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. After dropping households with incomplete data and
locations with fewer than 15 school-aged children, we obtain a sample of 656
clusters 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 (302 and 238, respectively) compared
to 45 for PKH (Indonesia), 58 for BDH (Ecuador), and 13 for CT-OVC (Kenya).
For the purposes of this analysis we assign equal weight to each cluster.

###  5.2. Implementation

Given the small to moderate cluster sizes, we choose an estimation approach
suited to sparsely sampled functional data, see also Appendix A. The main
difference to the densely sampled case is that the cluster-specific covariate
distribution $f_{g}(x)$ for the weights in (3.8) and (3.9) 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 differs from the idealized setup discussed in the theoretical
sections of the paper in that (a) 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 from the data at the cluster 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.

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 (3.4).

###  5.3. Results

We first report results for post-intervention outcomes using the estimated
IMSE-optimal predictors from (3.4). We assess their performance as predictors
at the level of the individual site as well as after aggregating sites within
each study.

|  | Prediction using IMSE-Optimal Basis Functions
---|---|---
baseline mean |  | K=1 | K=2 | K=3 | K=4 | K=5 | K=6
0.0915 |  | 0.2954 | 0.3251 | 0.3290 | 0.3562 | 0.3590 | 0.3713
|  | Prediction using Functional PC
baseline mean |  | K=1 | K=2 | K=3 | K=4 | K=5 | K=6
0.0915 |  | 0.3010 | 0.3173 | 0.3443 | 0.3319 | 0.3332 | 0.3197
Table 1. Prediction of Post-Intervention Outcomes $\mu_{g1}(x;1)-\mu(x;1)$
using the leading $K$ IMSE-optimal basis functions (top panel) and functional
PC (bottom panel) as predictors. Entries correspond to correlation
coefficients between actual mean at the site level and the baseline average
(first column) and the prediction based on the leading $K$ basis functions
(remaining columns).

Table 1 reports the correlation coefficient between the predicted average
effect at site $g$,
$(\hat{\mu}_{g1K}-\hat{\mu}_{1}):=\frac{1}{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,
$(\hat{\mu}_{g1}-\hat{\mu}_{1}):=\frac{1}{n_{g}}\sum_{i=1}^{n_{g}}(Y_{gi1}-\hat{\mu}_{1}(X_{gi}))$.
Since a natural alternative strategy would be to control for cluster-specific
average outcomes at baseline, we also report the correlation coefficient of
$(\hat{\mu}_{g1}-\hat{\mu}_{1})$ with the baseline mean
$(\hat{\mu}_{g0}-\hat{\mu}_{0}):=\frac{1}{n_{g}}\sum_{i=1}^{n_{g}}(Y_{gi0}-\hat{\mu}_{0}(X_{gi}))$
as a benchmark. We find a correlation coefficient of around $0.3$
(corresponding to an R-square of $0.09$) 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.

|  | Prediction using IMSE-Optimal Basis Functions
---|---|---
Study |  | actual |  | K=1 | K=2 | K=3 | K=4 | K=5 | K=6
MEX |  | -0.0067 |  | 0.0003 | -0.0006 | -0.0029 | -0.0058 | -0.0056 | -0.0067
MOR |  | 0.0213 |  | -0.0012 | 0.0029 | 0.0088 | 0.0124 | 0.0121 | 0.0135
IDN |  | 0.0008 |  | -0.0005 | 0.0020 | -0.0015 | 0.0078 | 0.0073 | 0.0119
KEN |  | 0.0219 |  | 0.0526 | 0.0521 | 0.0548 | 0.0492 | 0.0481 | 0.0456
ECU |  | -0.0456 |  | -0.0106 | -0.0237 | -0.0286 | -0.0194 | -0.0188 | -0.0163
|  | Prediction using Functional PC
Study |  | actual |  | K=1 | K=2 | K=3 | K=4 | K=5 | K=6
MEX |  | -0.0067 |  | 0.0000 | -0.0029 | -0.0054 | -0.0052 | -0.0060 | -0.0060
MOR |  | 0.0213 |  | -0.0005 | 0.0060 | 0.0120 | 0.0120 | 0.0128 | 0.0128
IDN |  | 0.0008 |  | -0.0017 | -0.0027 | 0.0084 | 0.0069 | 0.0105 | 0.0111
KEN |  | 0.0219 |  | 0.0514 | 0.0531 | 0.0525 | 0.0525 | 0.0524 | 0.0472
ECU |  | -0.0456 |  | -0.0098 | -0.0123 | -0.0222 | -0.0242 | -0.0213 | -0.0206
Table 2. Predictions of post-intervention outcomes $\mu_{g1}(x;1)-\mu(x;1)$
using the leading $K$ IMSE optimal basis functions (top panel) and functional
PC (bottom panel) as predictors. Rows correspond to averages across sites for
each of the five studies (Progresa, Tayssir, PKH, CT-OVC, BDH), “actual”
corresponds to the empirical mean of $\mu_{g1}(x;1)-\mu(x;1)$.

In Table 2, 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
$\frac{1}{|\mathcal{G}_{s}|}\sum_{g\in\mathcal{G}_{s}}(\hat{\mu}_{g1K}-\hat{\mu}_{1})$
to their “realized” empirical counterparts,
$\frac{1}{|\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)$. 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.

Figure 1. Estimated scores for leading $K=2$ IMSE-optimal predictors of
$\mu_{g1}(x)$.

Finally, Figure 1 plots the estimated scores corresponding to the leading two
basis functions, $\hat{m}_{g1},\hat{m}_{g2}$ together with 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.

|  | Prediction using IMSE-Optimal Basis Functions
---|---|---
Study |  | Actual |  | K=1 | K=2 | K=3 | K=4 | K=5 | K=6
MEX |  | 0.0046 |  | -0.0007 | -0.0008 | 0.0009 | 0.0009 | 0.0022 | 0.0020
MOR |  | -0.0153 |  | 0.0076 | 0.0083 | 0.0059 | 0.0059 | 0.0073 | 0.0091
IDN |  | -0.0423 |  | -0.0034 | -0.0073 | -0.0101 | -0.0095 | -0.0182 | -0.0257
KEN |  | -0.0257 |  | -0.0614 | -0.0605 | -0.0461 | -0.0459 | -0.0563 | -0.0582
ECU |  | -0.0240 |  | -0.0082 | -0.0076 | -0.0162 | -0.0157 | -0.0260 | -0.0265
|  | Prediction using Functional PC
Study |  | actual |  | K=1 | K=2 | K=3 | K=4 | K=5 | K=6
MEX |  | 0.0046 |  | 0.0002 | 0.0047 | 0.0026 | 0.0019 | 0.0047 | 0.0030
MOR |  | -0.0153 |  | -0.0010 | -0.0138 | -0.0079 | -0.0069 | -0.0088 | -0.0062
IDN |  | -0.0423 |  | 0.0006 | -0.0012 | 0.0051 | 0.0086 | 0.0009 | -0.0002
KEN |  | -0.0257 |  | 0.0205 | 0.0086 | -0.0322 | -0.0320 | -0.0415 | -0.0492
ECU |  | -0.0240 |  | -0.0021 | 0.0034 | 0.0056 | 0.0058 | -0.0069 | 0.0011

Table 3. Prediction of Conditional ATE $\tau_{g}(x)-\tau(x)$ averaging across
sites for each of the five studies (Progresa, Tayssir, PKH, CT-OVC, BDH),
using the leading $K$ IMSE optimal basis functions (top panel) and functional
PC (bottom panel) as predictors. “Actual” corresponds to the empirical mean of
$\tau_{g}(x)-\tau(x)$.

We next repeat the analysis for prediction of site-specific treatment effects
$\tau_{g}(x)-\tau(x)$, both using IMSE-optimal basis functions and functional
PC as predictors (Table 3). Since in all five RCTs, treatment was randomized
at the cluster level, 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. For any of these comparisons, it should also
be noted that both types of predictions are based on the unconfounded location
assumption (Assumption 2.2), whereas realized conditional effects also reflect
systematic differences between studies that can’t be predicted by
extrapolating intra-study variation among sites. Interestingly, the estimated
scores for predicting conditional ATEs plotted in Figure 2 appear to be more
closely aligned across the five studies than for predicting post-intervention
outcomes in Figure 1.

Figure 2. Estimated scores for leading $K=2$ IMSE-optimal predictors of
$\tau_{g}(x)$.

##  6\. 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.

##  Appendix A Extensions

The asymptotic results in the main text assume that we observe a dense sample
from each cluster, i.e. $n_{g}\rightarrow\infty$ for each $g$, and that
treatment was assigned at random within each cluster. These assumptions were
only used in the derivation of convergence rates for the nonparametric
estimators for $\mu(x;d)$ and $H(x_{1},x_{2};d_{1},d_{2})$. In this appendix
we discuss strategies for estimation and asymptotic rates for the case of
sparse samples as well as cluster-based randomization.

###  A.1. Sparsely Sampled Clusters

If the number of observed units in each cluster $n_{g}$ is bounded or grows at
a slow rate relative to $G$, we have to adapt the strategy for estimating the
mean and covariance functions of $\mu_{g}(x;d)$ along the lines proposed by
Yao, Müller, and Wang (2005a) for the problem of functional principal
components.

A challenge relative to the densely sampled case is that the site-specific
distribution of attributes $f_{g}(x)$ can no longer be estimated consistently.
Instead, we assume that the cluster-specific distribution of attributes
$X_{gi}$ is either known or can be estimated parametrically with sufficient
precision from the $n_{g}$ observations in each cluster, so that individual
observations can be reweighted accordingly to match the target distribution
$f_{0}(x)$ in each cluster.444For the combined studies of conditional cash
transfers in the empirical application, this is unproblematic for some of the
demographic variables, including the child’s age and gender. On the other
hand, the means and variances of log per capita household expenditure vary
between sites, but separate histograms for each study suggest that the log-
normal distribution provides a plausible parametric approximation to the
marginal distribution.

Specifically, let

$\hat{\mu}(x;d):=\arg_{b_{0}}\min_{b_{0},b_{1}}\sum_{g=1}^{G}\sum_{i=1}^{n_{g}}w_{gi}(x;d,h_{\mu})(Y_{gi}-b_{0}-b_{1}(x-X_{gi}))^{2}$

for $\mu(x;d)$. As before, $K(u)$ is a kernel function satisfying Assumption
4.1, and the bandwidth $h>0$ is chosen according to sample size $G,n_{g}$, the
dimension of $X_{gi}$, and assumed smoothness of the estimands. The weights
$w_{gi}(x;d,h):=1\textnormal{l}\\{D_{gi}=d\\}\frac{f_{0}(x)}{f_{g}(x)}K\left(\frac{X_{gi}-x}{h}\right)$
are assumed to be known, and we define

$\delta_{Gn}:=\min_{d\in\\{0,1\\}}\min_{g=1,\dots,G}\inf_{\|x\|\leq
c_{n}}\frac{f_{g}(x;\mathbf{d})}{f_{0}(x;\mathbf{d})}>0$

We also let

$\hat{H}(x_{1},x_{2};d_{1},d_{2}):=\arg_{b_{0}^{(g)}}\min_{b_{0}^{(g)},b_{11}^{(g)},b_{12}^{(g)}}\sum_{g=1}^{G}\sum_{j\neq
i}H_{gij}(x_{1},x_{2},\mathbf{b})w_{gi}(x_{1};d,h_{H})w_{gj}(x_{2};d,h_{H})$

for an appropriately chosen bandwidth $h_{H}$, 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}$

Then using arguments parallel to the proof of Theorem 4.1, the convergence
rates of the local linear estimators for $\mu(x;\mathbf{d})$ and the
covariance operator $H_{d_{1}d_{2}}(x_{1},x_{2})$ are given by

$\displaystyle\sup_{d,x_{1}}\left|\hat{\mu}(x_{1},d)-\mu(x_{1},d)\right|$
$\displaystyle=$ $\displaystyle O_{p}\left(r_{G\mu}\right)$
$\displaystyle\sup_{d_{1},d_{2},x_{1},x_{2}}\left|\hat{H}(x_{1},x_{2};d_{1},d_{2})-H(x_{1},x_{2};d_{1},d_{2})\right|$
$\displaystyle=$ $\displaystyle O_{p}\left(r_{GH}\right)$

where $r_{G\mu}=\delta_{Gn}^{-1}\left(h_{\mu}^{2}+\left(\frac{\log
G}{Gh_{\mu}^{d}}\right)^{1/2}\right)$ and
$r_{GH}=\delta_{Gn}^{-1}\left(h_{H}^{2}+\left(\frac{\log
G}{Gh_{H}^{2d}}\right)^{1/2}\right)$. A similar result was proven under
slightly different conditions in Theorem 1 by Yao, Müller, and Wang (2005a).

Comparing these rates to those for the densely sampled case in Theorem 4.1,
the contributions of order $\frac{1}{G}$ are now unambiguously dominated by
the remaining errors. Furthermore, in the rate for the covariance kernel, the
bandwidth $h_{H}$ now appears at the power $h_{H}^{2d}$ (compared to
$h_{H}^{d}$ in the dense case). This is a consequence of the fact that the
number of unit pairs in each cluster no longer increases to infinity under the
sparsely sampled case. As a result, the optimal bandwidths for estimating
$\mu(x;d)$ and $H(x_{1},x_{2};d_{1},d_{2})$ are generally at different rates
and should therefore also be chosen separately in this case. As in the densely
sampled case, we can then use Theorem 4.2 to determine rates for
eigenfunctions and eigenvalues. Unbiased estimators for the scores
$m_{gk},t_{gk}$ are also available, however consistent estimation requires
$n_{g}$ growing large for the sites of interest.

###  A.2. Cluster-Based Randomization

We can also adapt the approach to the scenario in which treatment assignment
is randomized at the cluster level, but a baseline survey of pre-intervention
outcomes is available. Specifically, we let $D_{gi}\equiv D_{g}$ be the
assigned treatment for all units in cluster $g$, and $Y_{git}(d)$ denote the
potential value for the policy variable $D_{git}=d$ unit $i$ at site $g$ in
period $t=0,1$. We then consider the problem of predicting post-intervention
conditional average treatment effects

$\tau_{g1}(x):=\mathbb{E}\left[Y_{gi1}(1)-Y_{gi1}(0)|X_{gi}=x\right]$

based on $\mu_{g0}(x):=\mathbb{E}\left[Y_{gi0}(0)|X_{gi}=x\right]$, where the
covariates $X_{gi}$ are also assumed to be measured at baseline.

With cluster-level randomization, we only observe potential values
$Y_{gi1}(1)$ for all units in a cluster $g$ assigned to treatment, $D_{g}=1$,
or only $Y_{gi1}(0)$ for all units in a control cluster with $D_{g}=0$.
However under ignorable assignment,

$D_{g}\bot\bot Y_{gi0}(0),Y_{gi1}(0),Y_{gi1}(1)|X_{gi}$

and independent assignment across clusters, we can estimate the covariance
function$H_{\mu\tau}(x_{1},x_{2}):=\textnormal{Cov}(\mu_{g0}(x_{1}),\tau_{g}(x_{2}))$
consistently as $G\rightarrow\infty$ via

$\hat{H}_{\mu\tau}(x_{1},x_{2}):=\hat{H}(x_{1},x_{2};0,1)-\hat{H}(x_{1},x_{2};0,0).$

where $\hat{H}(x_{1},x_{2};0,1)$ and $\hat{H}(x_{1},x_{2};0,0)$ are
nonparametric estimators for
$\textnormal{Cov}(\mu_{g0}(x_{1}),\mu_{g1}(x_{2};1))$
and$\textnormal{Cov}(\mu_{g0}(x_{1}),\mu_{g1}(x_{2};0))$, respectively,
obtained separately from the treatment and control clusters.

With minor modifications of the proofs of Theorem 4.1 and Theorem 1 in Yao,
Müller, and Wang (2005a) we can then obtain the convergence rates
$r_{Gn}=\frac{1}{\sqrt{G}}+\delta_{Gn}^{-1}\left(h^{2}+\left(\frac{\log
n}{Gnh^{d}}\right)^{1/2}\right)$ for the densely sampled case where
$n_{g}\equiv n\rightarrow\infty$ for each cluster. For the sparsely sampled
case where $n_{g}$ remains fixed, we obtain
$r_{G\mu}=\delta_{Gn}^{-1}\left(h_{\mu}^{2}+\left(\frac{\log
G}{Gh_{\mu}^{d}}\right)^{1/2}\right)$ and
$r_{GH}=\delta_{Gn}^{-1}\left(h_{H}^{2}+\left(\frac{\log
G}{Gh_{H}^{2d}}\right)^{1/2}\right)$, so that the rates coincide with the case
of within cluster-randomization. The eigenfunctions in (3.4) can then be
obtained from Theorem 4.2 as before.

##  Appendix B Proofs for Section 3

Before proving the main result, we give a characterization of the integrated
mean square error of projection into closed linear subspaces of
$L_{2}(\mathcal{X})$. Specifically, we consider the mean square error of
projection of $\tau_{g}$ with respect to a general linear transformation
$\mu_{g}$.

We consider the problem where the domain of the operator $B$ is restricted to
$\mathcal{H}$, a closed linear subspace $\mathcal{H}$ of
$\mathcal{N}^{\perp})$, which was defined as the orthogonal complement of the
null space of $T_{\mu\mu}$, $\mathcal{N}:=\ker(T_{\mu\mu})$. It is known that
for any linear operator $T:L_{2}(\mathcal{X})\rightarrow L_{2}(\mathcal{X})$,
$\textnormal{im}(T)$ and $\ker(T)$ are linear subspaces of
$L_{2}(\mathcal{X})$, and that $\ker(T)$ is closed if the linear operator $T$
is bounded. Since $\mathcal{H}$ is the orthogonal complement of the null space
$\ker(T_{\mu\mu})$, it is a closed linear subspace of $L_{2}(\mathcal{X})$, so
this allows for the choice $\mathcal{H}\equiv\mathcal{N}^{\perp}$. We also let
the operator $P:L_{2}(\mathcal{X})\rightarrow\mathcal{H}$ to denote the
orthogonal projection on $\mathcal{H}$. Since $\mathcal{H}$ is closed by
assumption, that projection is well-defined by the Classical Projection
Theorem (Theorem 2 on p.51 in Luenberger (1969)).

We can then characterize the integrated mean square error of projection as
follows:

######  Theorem B.1.

(Integrated MSE of Projection) There exists a best linear predictor for
$\tau_{g}$ of the form (3.1) based on $\mu_{g}$ that is of the form

$\tau_{g}^{*}(x):=\int\mu_{g}(s)\beta(s,x)f_{0}(s)ds$

If furthermore the operator
$PT_{\mu\mu}A^{*}P:\mathcal{H}\rightarrow\mathcal{H}$ possesses an inverse,
the minimized mean square error of prediction satisfies

$\int\mathbb{E}\left[(\tau_{g}(x)-\tau_{g}^{*}(x))^{2}\right]f_{0}(x)dx=\|\tau_{g}\|-\textnormal{tr}\left(T_{\mu\tau}^{*}P(PT_{\mu\mu}P)^{-1}PT_{\mu\tau}\right)$
(B.1)

For the definition of the inverse in (B.1), note that the operator and
therefore its inverse are understood to be restricted to the subspace
$\mathcal{H}$.

Proof: We first consider the case $P=\textnormal{Id}$. The set $\mathcal{M}$
of linear predictors of the form (3.1) can be identified with the Hilbert
space $L_{2}(\mathcal{X}\times\mathcal{X})$ endowed with the scalar product
$\langle T_{1},T_{2}\rangle=tr(T_{1}^{*}T_{2})=\int\int
g(x_{2},x_{1})h(x_{1},x_{2})f_{0}(x_{1})f_{0}(x_{2})dx_{1}dx_{2}$ and the
trace norm induced by that scalar product.

To obtain a representation of the minimized projection error, we derive a
projection analog for the regression model developed in He, Müller, and Wang
(2000). We first show that $\mathcal{M}$ is a closed linear subspace of
$L_{2}(\mathcal{X})$: Consider a random element $\phi_{g}\in
L_{2}(\mathcal{X})$ and define
$H_{\phi_{f},\mu_{g}}(x_{1},x_{2}):=\mathbb{E}[\phi_{g}(x_{1})\mu_{g}(x_{2};0)]$
and let $T_{\phi\mu}$ denote the corresponding Hilbert-Schmidt integral
operator. We then have that $\phi_{g}$ is orthogonal to $B\mu_{g}$ if

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\langle\phi_{g},B\mu_{g}\rangle$ $\displaystyle=$
$\displaystyle\mathbb{E}\left[\int\phi_{g}(x_{1})\mu_{g}(x_{2})\beta(x_{1},x_{2})f_{0}(x_{1})f_{0}(x_{2})dx_{1}dx_{2}\right]$
$\displaystyle=$
$\displaystyle\int\mathbb{E}\left[\phi_{g}(x_{1})\mu_{g}(x_{2})\right]\beta(x_{1},x_{2})f_{0}(x_{1})f_{0}(x_{2})dx_{1}dx_{2}$
$\displaystyle=$ $\displaystyle\int
H_{\phi_{g},\mu_{g}}(x_{1},x_{2})\beta(x_{1},x_{2})f_{0}(x_{1})f_{0}(s)dx_{1}dx_{2}=tr\left(T_{\phi\mu}B\right)$

Since $B$ may in particular include the identity, $\phi_{g}$ is orthogonal to
$\mathcal{M}$ for any $B\in L_{2}(\mathcal{X}\times\mathcal{X})$ if and only
if $H_{\phi_{g},\mu_{g}}(x_{1},x_{2})\equiv 0$. Since $\mathcal{M}$ is the
orthogonal complement of the set of all such vectors $\phi_{g}\in
L_{2}(\mathcal{X})$, it is in particular a closed linear subspace of
$L_{2}(\mathcal{X})$.

By the classical projection theorem (Theorem 2 in Luenberger (1969), p.51), it
then follows that there exists a unique minimizer
$\tau_{0}^{*}\in\mathcal{M}$. Furthermore, that minimizer also satisfies the
orthogonality condition $\langle m,\tau_{g}-\tau_{0}^{*}\rangle$ for any
$m\in\mathcal{M}$. It follows that

$\displaystyle\min_{B\in
L_{2}(\mathcal{X}\times\mathcal{X})}\left\|\tau_{g}-B\mu_{g}\right\|^{2}$
$\displaystyle=$
$\displaystyle\|\tau_{g}\|^{2}-2\langle\tau_{g},\tau_{0}^{*}\rangle+\|\tau_{0}^{*}\|^{2}$
(B.2) $\displaystyle=$
$\displaystyle\|\tau_{g}\|^{2}-2(\|\tau_{g}^{*}\|^{2}+\langle\tau_{g}-\tau_{g}^{*},\tau_{g}^{*}\rangle)+\|\tau_{g}^{*}\|^{2}$
$\displaystyle=$
$\displaystyle\left\|\tau_{g}\right\|^{2}-\left\|\tau_{g}^{*}\right\|^{2}$

We next characterize the optimal solution $\tau_{g}^{*}$ in terms of the
operators $T_{\mu\mu}$ and $T_{\mu\tau}$. Suppose that $B_{0}$ with kernel
$\beta_{0}^{*}(x_{1},x_{2})$ is such that $\tau_{g}^{*}=B_{0}\mu_{g}$. Then,
we have for any $h\in L_{2}(\mathcal{X})$ that

$\displaystyle(T_{\mu\tau}h)(x_{1})$ $\displaystyle=$
$\displaystyle\int\mathbb{E}[\mu_{g}(x_{1})\tau_{g}(x_{2})]h(x_{2})f_{0}(x_{2})dx_{2}$
(B.3) $\displaystyle=$
$\displaystyle\int\mathbb{E}\left[\mu_{g}(x_{1})\left[(B_{0}\mu_{g})(x_{2})+\left\\{\tau_{g}(x_{2})-(B_{0}\mu)(x_{2})\right\\}\right]\right]h(x_{2})f_{0}(x_{2})dx_{2}$
$\displaystyle=$
$\displaystyle\int\mathbb{E}[\mu_{g}(x_{1})\mu_{g}(x_{2})]\beta_{0}^{*}(x_{2},x_{3})f_{0}(x_{2})dx_{2}h(x_{3})f_{0}(x_{3})dx_{3}$
$\displaystyle+\int\mathbb{E}\left[\mu_{g}(x_{1})\left\\{\tau_{g}(x_{3})-(B_{0}\mu)(x_{3})\right\\}\right]h(x_{3})f_{0}(x_{3})dx_{3}$
$\displaystyle=$ $\displaystyle\left(T_{\mu\mu}B_{0}h\right)(x_{1})+0$

where the last equality follows from orthogonality of the projection error,
noting that $B\mu_{g}:=\textnormal{Id}\mu_{g}$ is in $\mathcal{M}$.

Now suppose that the minimum is attained at both $B_{0}\mu_{g}\in\mathcal{M}$
and $B_{1}\mu_{g}\in\mathcal{M}$. Then by orthogonality of the projection
error,

$\displaystyle\|\tau_{g}-B_{1}\mu_{g}\|^{2}$ $\displaystyle=$
$\displaystyle\|\tau_{g}-B_{0}\mu_{g}+(B_{0}-B_{1})\mu_{g}\|^{2}$
$\displaystyle=$
$\displaystyle\|\tau_{g}-B_{0}\mu_{g}\|^{2}+2\langle\tau_{g}-B_{0}\mu_{g},(B_{0}-B_{1})\mu_{g}\rangle+\|(B_{0}-B_{1})\mu_{g}\|^{2}$
$\displaystyle=$
$\displaystyle\|\tau_{g}-B_{0}\mu_{g}\|^{2}+\|(B_{0}-B_{1})\mu_{g}\|^{2}$

Hence the minimum can be attained at both $B_{1}$ and $B_{0}$ iff

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\|(B_{0}-B_{1})\mu_{g}\|^{2}$
$\displaystyle=$
$\displaystyle\int(\beta_{0}(x_{1},x_{2})-\beta_{1}(x_{1},x_{2}))\mathbb{E}[\mu_{g}(x_{1})\mu_{g}(x_{3})](\beta_{0}(x_{3},x_{2})-\beta_{1}(x_{3},x_{2}))f_{0}(x_{1})f_{0}(x_{3})f_{0}(x_{2})dx_{1}dx_{2}dx_{3}$
$\displaystyle=$
$\displaystyle\textnormal{tr}\left((B_{0}-B_{1})T_{\mu\mu}(B_{0}^{*}-B_{1}^{*})\right)$

or equivalently iff $(B_{0}-B_{1})\in\ker(T_{\mu\mu})$. In particular, the
orthogonal projection $B_{0}^{*}$ of any solution $B\in
L_{2}(\mathcal{X}\times\mathcal{X})$ onto the closed subspace
$\mathcal{N}^{\perp}\times L_{2}(\mathcal{X})$ exists and is a solution to the
same minimization problem. So without loss of generality we can restrict our
attention to the minimization problem

$\min_{B\in\mathcal{N}^{\perp}\times
L_{2}(\mathcal{X})}\left\|\tau_{g}-B\mu_{g}\right\|$ (B.4)

Noticing that restricted to its range $\mathcal{N}^{\perp}$, the operator
$T_{\mu\mu}$ possesses an inverse, we can solve the normal equations (B.3) to
obtain

$B_{0}^{*}h=T_{\mu\tau}^{*}T_{\mu\mu}^{-1}h$

for any $h\in\mathcal{N}^{\perp}$. Substituting this expression into (B.2), we
therefore obtain

$\min_{B\in
L_{2}(\mathcal{X}\times\mathcal{X})}\left\|\tau_{g}-B\mu_{g}\right\|^{2}=\|\tau_{g}\|^{2}-\textnormal{tr}\left(T_{\mu\tau}^{*}T_{\mu\mu}^{-1}T_{\mu\tau}\right)$

establishing the claim for $P_{1}=\textnormal{Id}$. For the general case,
notice that $\mathcal{H}$ was a closed linear subspace of
$\mathcal{N}^{\perp}$ and $A$ was assumed to be invertible on $\mathcal{H}$,
so that the argument continues to apply when restricted to the subspace of
linear predictors based on $P\mu_{g}$∎

###  B.1. Proof of Lemma 3.1

We use the formula from (B.1) with
$\tilde{\mathcal{H}}_{K}:=\textnormal{span}\left(\tilde{\phi}_{1},\dots,\tilde{\phi}_{K}\right)$
as the linear subspace of $L_{1}(\mathcal{X})$ spanned by the $K$ orthonormal
functions $\tilde{\phi}_{1}(x),\dots,\tilde{\phi}_{K}(x)$.

We first consider the operator
$\tilde{S}_{K}:=\tilde{P}_{K}T_{\mu\mu}\tilde{P}_{K}$. Since
$\tilde{\mathcal{H}}_{K}$ is finite-dimensional, we can identify
$\tilde{S}_{K}$ with a $K\times K$ matrix of coefficients with respect to the
basis $\tilde{\phi}_{1},\dots,\tilde{\phi}_{K}$ of $\tilde{\mathcal{H}}_{K}$,

$\tilde{S}_{K}:=\left(\langle\tilde{\phi}_{k},T_{\mu\mu}\tilde{\phi}_{l}\rangle\right)_{k,l}$

Noting that $T_{\mu\mu}$ is injective on $\mathcal{N}^{\perp}$ we will impose
the normalization

$\langle\tilde{\phi}_{k},T_{\mu\mu}\tilde{\phi}_{l}\rangle=\delta_{k,l}$ (B.5)

so that $\tilde{S}_{K}=I_{K}$, the $K$-dimensional identity matrix.

Furthermore, evaluating the trace using an arbitrarily chosen orthonormal
system for $L_{2}(\mathcal{X})$, we can also verify that

$\displaystyle\textnormal{tr}(T_{\mu\tau}^{*}\tilde{P}_{K}\tilde{P}_{K}T_{\mu\tau})$
$\displaystyle=$ $\displaystyle\int\sum_{k=1}^{K}\left(\int
H_{\mu\tau}(x_{1},x_{2})H_{\mu\tau}(x_{3},x_{1})\tilde{\phi}_{k}(x_{2})\tilde{\phi}_{k}(x_{3})f_{0}(x_{2})f_{0}(x_{3})dx_{2}dx_{3}\right)f_{0}(x_{1})dx_{1}$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{K}\langle\tilde{\phi}_{k},T_{\mu\tau}T_{\mu\tau}^{*}\tilde{\phi}_{k}\rangle$

Hence, using the formula from Theorem B.1,

$\int\mathbb{E}\left[(\tau_{g}(x)-\tau_{g}^{*}(x))^{2}\right]f_{0}(x)dx=\sum_{k=1}^{K}\langle\tilde{\phi}_{k},T_{\mu\tau}T_{\mu\tau}^{*}\tilde{\phi}_{k}\rangle$

for any collection of functions $\tilde{\phi}_{1},\dots,\tilde{\phi}_{K}$
satisfying the constraint (B.5).

Hence, the mean-square error optimal basis functions are determined by the
quadratic program

$\inf_{\phi_{1},\dots,\phi_{K}}\sum_{k=1}^{K}\langle\phi_{k},T_{\mu\tau}T_{\mu\tau}^{*}\phi_{k}\rangle\hskip
22.76228pt\textnormal{s.t.
}\langle\phi_{k},T_{\mu\mu}\phi_{l}\rangle=\delta_{kl}$ (B.6)

which establishes the claim∎

###  B.2. Proof of Theorem 3.1

From the definition of $IMSE_{K}^{*}$, there exist
$\tilde{\phi}_{1\varepsilon},\dots,\tilde{\phi}_{K\varepsilon}\in
L_{2}(\mathcal{X})$, such that
$\langle\tilde{\phi}_{k\varepsilon},T_{\mu\mu}\phi_{l\varepsilon}\rangle=\delta_{kl}$
and

$IMSE_{K}\left[\tilde{\phi}_{1\varepsilon},\dots,\tilde{\phi}_{K\varepsilon}\right]\leq
IMSE_{K}^{*}+\frac{\varepsilon}{2}$

Since $\tilde{\phi}_{1\varepsilon},\dots,\tilde{\phi}_{K\varepsilon}\in
L_{2}(\mathcal{X})$, we can find $L_{\varepsilon}<\infty$ such that
$\|\tilde{\phi}_{k\varepsilon}\|\leq L_{\varepsilon}$.

These $K$ functions do not necessarily satisfy the regularized orthogonality
constraint, rather we find “close” alternative functions
$\bar{\phi}_{1\varepsilon},\dots,\bar{\phi}_{K\varepsilon}$ such that
$\langle\check{\phi}_{k\varepsilon},T_{\mu\mu
a}\check{\phi}_{l\varepsilon}\rangle=\delta_{kl}$. Since the operator
$T_{\mu\mu a}$ is Hermitian and positive, we can define a scalar product
$\langle u,v\rangle_{a}:=\langle u,T_{\mu\mu a}v\rangle$. We first obtain $K$
functions $\check{\phi}_{1\varepsilon},\dots,\check{\phi}_{K\varepsilon}$ that
satisfy the orthogonality constraints
$\langle\check{\phi}_{k\varepsilon},\check{\phi}_{l\varepsilon}\rangle=0$ for
all $k\neq l$. To that end, we use the Gram-Schmidt procedure with respect to
the scalar product $\langle\cdot,\cdot\rangle_{a}$, where
$\check{\phi}_{1\varepsilon}:=\tilde{\phi}_{1\varepsilon}$, and

$\check{\phi}_{k\varepsilon}:=\tilde{\phi}_{k\varepsilon}-\sum_{l=1}^{k-1}\frac{\langle\tilde{\phi}_{k\varepsilon},\check{\phi}_{l\varepsilon}\rangle_{a}}{\langle\check{\phi}_{l\varepsilon},\check{\phi}_{l\varepsilon}\rangle_{a}}\check{\phi}_{l\varepsilon}$

Defining $c_{kk}=1$ and
$c_{lk}:=-\frac{\langle\tilde{\phi}_{k\varepsilon},\check{\phi}_{l\varepsilon}\rangle_{a}}{\langle\check{\phi}_{l\varepsilon},\check{\phi}_{l\varepsilon}\rangle_{a}}$
for $k\neq l$, we can confirm by induction that

$\check{\phi}_{k\varepsilon}=\check{\phi}_{k\varepsilon}-\sum_{l=1}^{k-1}d_{lk}\tilde{\phi}_{l\varepsilon}$

where $d_{lk}=\sum_{\mathcal{J}(l,k)}\prod_{s=2}^{|k-l|}c_{j_{s}j_{s-1}}$, and
$\mathcal{J}(l,k)$ is the set of all $(j_{1},\dots,j_{s})$ such that
$s\leq|k-l|$ and $j_{1}<j_{2}\dots<j_{s}$. For notational convenience, we also
define $d_{kk}=1$.

We now recursively bound $c_{kl}$ and $d_{kl}$ for $k\neq l$. Due to the
constraint on $\tilde{\phi}_{1\varepsilon},\dots,\tilde{\phi}_{K\varepsilon}$
it follows for any $k\neq l$,

$\langle\tilde{\phi}_{k\varepsilon},\tilde{\phi}_{l\varepsilon}\rangle_{a}=\langle\tilde{\phi}_{k\varepsilon},T_{\mu\mu}\tilde{\phi}_{l\varepsilon}\rangle+a\langle\tilde{\phi}_{k\varepsilon},\tilde{\phi}_{l\varepsilon}\rangle=a\langle\tilde{\phi}_{k\varepsilon},\tilde{\phi}_{l\varepsilon}\rangle$

Hence we can calculate the scalar products

$\displaystyle\langle\tilde{\phi}_{k\varepsilon},\check{\phi}_{l\varepsilon}\rangle_{a}$
$\displaystyle=$ $\displaystyle
a\sum_{j=1}^{l}d_{lj}\langle\tilde{\phi}_{k\varepsilon},\tilde{\phi}_{j\varepsilon}\rangle$
$\displaystyle\langle\check{\phi}_{k\varepsilon},\check{\phi}_{k\varepsilon}\rangle_{a}$
$\displaystyle=$ $\displaystyle
1+\sum_{l=1}^{k-1}d_{kl}^{2}+a\sum_{j=1}^{k-1}\sum_{l=1}^{k-1}d_{kj}d_{kl}\langle\tilde{\phi}_{j\varepsilon},\tilde{\phi}_{l\varepsilon}\rangle$

for any $k\neq l$. In particular, denoting
$\bar{d}_{k}:=\max\\{|d_{kl}|:k>l\\}$, we can bound

$\displaystyle\left|\langle\tilde{\phi}_{k\varepsilon},\check{\phi}_{l\varepsilon}\rangle_{a}\right|$
$\displaystyle\leq$ $\displaystyle a(1+(k-1)\bar{d}_{k})L_{\varepsilon}^{2}$
$\displaystyle\left|\langle\check{\phi}_{k\varepsilon},\check{\phi}_{k\varepsilon}\rangle_{a}-1-\sum_{l=1}^{k-1}d_{kl}^{2}\right|$
$\displaystyle\leq$ $\displaystyle
aL_{\varepsilon}^{2}(1+(k-1)^{2}\bar{d}_{k}^{2})$

Hence for $a$ satisfying $2a^{2}\leq(L_{\varepsilon}^{2}K^{2})^{-1}$, we can
bound $\bar{d}_{2}:=|c_{21}|\leq 2aL_{\varepsilon}^{2}$. It then follows by
induction over $k$ that $|c_{kl}|\leq 4aL_{\varepsilon}^{2}$ and

$\bar{d}_{k}\leq 8aL_{\varepsilon}^{2}$

for each $k=1,\dots,K$.

To obtain functions
$\bar{\phi}_{1\varepsilon},\dots,\bar{\phi}_{K\varepsilon}$ with the desired
property, we furthermore need to impose the scale normalization
$\langle\bar{\phi}_{k\varepsilon},\bar{\phi}_{k\varepsilon}\rangle_{a}=1$.
Since
$\check{\phi}_{k\varepsilon}=\sum_{j=1}^{k}d_{kj}\tilde{\phi}_{k\varepsilon}$,

$\displaystyle\langle\check{\phi}_{k\varepsilon},\check{\phi}_{k\varepsilon}\rangle_{a}$
$\displaystyle=$
$\displaystyle\langle\sum_{j=1}^{k}\tilde{\phi}_{j\varepsilon},T_{\mu\mu}\langle\check{\phi}_{k\varepsilon},\check{\phi}_{k\varepsilon}\rangle\tilde{\phi}_{j\varepsilon}\rangle+a\langle\check{\phi}_{k\varepsilon},\check{\phi}_{k\varepsilon}\rangle$

Since $d_{kk}=1$, we can then bound

$|\langle\check{\phi}_{k\varepsilon},\check{\phi}_{k\varepsilon}\rangle_{a}-1|\leq
K\bar{d}_{k}^{2}+aL_{\varepsilon}^{2}\leq(K+8)aL_{\varepsilon}^{2}$

We can therefore form

$\bar{\phi}_{k}:=(\langle\check{\phi}_{k\varepsilon},\check{\phi}_{k\varepsilon}\rangle_{a})^{-1/2}\check{\phi}_{k\varepsilon}$

By construction, $\bar{\phi}_{1\varepsilon},\dots,\bar{\phi}_{K\varepsilon}$
satisfy the constraints
$\langle\bar{\phi}_{k\varepsilon},\bar{\phi}_{l\varepsilon}\rangle_{a}=\delta_{kl}$
and, using standard bounding arguments

$\|\bar{\phi}_{k\varepsilon}-\tilde{\phi}_{k\varepsilon}\|\leq
10KaL_{\varepsilon}^{3}$

for any sufficiently small value of $a>0$.

It remains to be shown by a continuity argument that the IMSE achieved by
$\bar{\phi}_{1\varepsilon},\dots,\bar{\phi}_{K\varepsilon}$ is greater than
$IMSE_{K}^{*}+\varepsilon$ for $a$ sufficiently close to zero. Using the
formula from (B.1), we can write the IMSE for a given set of basis functions
as

$IMSE_{K}[\phi_{1},\dots,\phi_{K}]=tr(C_{K}^{-1}A_{K})$

where $C_{K}$ and $A_{K}$ can be identified with $K\times K$ matrices
$A_{K}=(A_{kl})$ and $C_{K}=(C_{kl})$, and

$C_{kl}:=\langle\phi_{k},T_{\mu\mu}\phi_{l}\rangle$

and

$A_{kl}:=\langle\phi_{k},T_{\mu\tau}T_{\mu\tau}^{*}\phi_{l}\rangle$

By assumption, for the basis functions
$\tilde{\phi}_{1\varepsilon},\dots,\tilde{\phi}_{K\varepsilon}$ we have
$C_{kl}=\delta_{kl}$ so that the corresponding matrix is the $K$-dimensional
identity matrix. Furthermore, the operators $T_{\mu\mu}$ and
$T_{\mu\tau}T_{\mu\tau}^{*}$ are compact with largest eigenvalue equal to
$\nu_{1}<\infty$ and $\lambda_{1}<\infty$, respectively. Since
$\|\bar{\phi}_{k\varepsilon}-\tilde{\phi}_{k\varepsilon}\|\leq
aB_{9}K{\varepsilon}^{3}$, there exists a constant
$\kappa_{\varepsilon}<\infty$ such that

$|IMSE_{K}[\bar{\phi}_{1\varepsilon},\dots,\bar{\phi}_{K\varepsilon}]-IMSE_{K}[\tilde{\phi}_{1\varepsilon},\dots,\tilde{\phi}_{K\varepsilon}]|\leq\kappa_{\varepsilon}a$
(B.7)

for $a$ sufficiently small, which establishes the desired conclusion∎

We also note that the convergence rate for $|IMSE_{K}(a)-IMSE_{K}^{*}|$ as a
function of $a$ generally depends on the eigenvalues of $T_{\mu\mu}$ and
$T_{\mu\tau}$, where $L_{\varepsilon}$ in the preceding proof may diverge as
$\varepsilon\rightarrow 0$. We leave this question for future research.

### Proof of Proposition 3.1

Since $T_{\mu\mu}$ is injective on $\ker(T_{\mu\mu})^{\perp}$, any function
$\phi_{k}^{*}$ solving the generalized eigenvalue problem (3.3) at eigenvalue
$\lambda_{k}$ can be equivalently characterized by
$\phi_{k}^{*}:=T_{\mu\mu}^{-1/2}\chi_{k}^{*}$, where $\chi_{k}^{*}$ is an
eigenfunction

$T_{\mu\mu}^{-1/2}T_{\mu\tau}T_{\mu\tau}^{*}T_{\mu\mu}^{-1/2}\chi_{k}^{*}=\lambda_{k}\chi_{k}^{*}$

at the same eigenvalue $\lambda_{k}$. In what follows we also write
$S:=T_{\mu\mu}^{-1/2}T_{\mu\tau}T_{\mu\tau}^{*}T_{\mu\mu}^{-1/2}$. Without
loss of generality we also assume that eigenvalues are ordered according to
$\lambda_{1}\geq\lambda_{2}\geq\dots$, and $\lambda_{K}>0$.

We next argue that under the condition (3.6), $S$ is compact: it follows
immediately from the Karhunen-Loève expansion (3.5) that we can write

$\displaystyle(T_{\mu\mu}h)(x)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\infty}\mathbb{E}[\alpha_{gk}^{2}]\langle\xi_{k},h\rangle\xi_{k}(x)$
$\displaystyle(T_{\mu\tau}h)(x)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{\infty}\sum_{l=1}^{\infty}\mathbb{E}[\alpha_{gk}\beta_{gl}]\langle\zeta_{l},h\rangle\xi_{k}(x),$

with the analogous expression for the adjoint $T_{\mu\tau}^{*}$.

Since $\xi_{1},\xi_{2},\dots$ is a basis of $\ker(T_{\mu\mu})^{\perp}$, we can
therefore evaluate the trace of $S$ in terms of that basis to obtain

$\textnormal{tr}(S)=\sum_{k=1}^{\infty}\sum_{l=1}^{\infty}\frac{\mathbb{E}[\alpha_{gk}\beta_{gl}]^{2}}{\mathbb{E}[\alpha_{gk}^{2}]}$

which is finite by (3.6). In particular, the sequence of compact operators
$S_{K}$ defined as the analog of $S$ after truncating the respective Karhunen-
Loève expansions for $\mu_{g}$ and $\tau_{g}$ after the $K$ leading principal
components converges to $S$ under the trace norm, so that $S$ is also compact.
Since $S$ is self-adjoint and nonnegative and therefore has a discrete
spectrum of nonnegative eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\dots\geq
0$ with associated eigenfunctions $\chi_{1}^{*},\chi_{2}^{*},\dots\in
L_{2}(\mathcal{X})$.

It remains to be shown that the norm of
$\phi_{k}^{*}:=T_{\mu\mu}^{-1/2}\chi_{k}^{*}$ is also bounded. Since
$\chi_{k}^{*}$ is an eigenfunction of $S$ at eigenvalue $\lambda_{k}>0$, we
can write

$\phi_{k}^{*}:=\lambda_{k}^{-1}T_{\mu\mu}^{-1/2}S\chi_{k}^{*}.$

Hence, the norm of $\phi_{k}^{*}$ is given by

$\|\phi_{k}^{*}\|^{2}=\lambda_{k}^{-2}\langle
T_{\mu\mu}^{-1/2}S\chi_{k}^{*},T_{\mu\mu}^{-1/2}S\chi_{k}^{*}\rangle=\lambda_{k}^{-2}\langle\chi_{k}^{*},ST_{\mu\mu}^{-1}S\chi_{k}^{*}\rangle$

noting that $S$ and $T_{\mu\mu}$ are self-adjoint. Moreover, expressing
$T_{\mu\mu}$ and $T_{\mu\tau}$ in terms of the Karhunen-Loève representations
(3.5), we can use the Cauchy-Schwarz inequality to bound

$\textnormal{tr}(ST_{\mu\mu}^{-1}S)\leq\left(\sum_{k=1}^{\infty}\sum_{l=1}^{\infty}\frac{\mathbb{E}[\alpha_{gk}\beta_{gl}]^{2}}{\mathbb{E}[\alpha_{gk}^{2}]^{3/2}}\right)^{2}$

which is finite under (3.6). Noting that $S^{r}T_{\mu\mu}^{-1}S^{r}$ is self-
adjoint, this establishes that the operator is also trace-class and therefore
compact with operator norm bounded by $\textnormal{tr}(ST_{\mu\mu}^{-1}S)$.

In particular, since $\lambda_{1}\geq\dots\lambda_{K}>0$ by assumption, we can
bound

$\max_{k=1,\dots,K}\|\phi_{k}^{*}\|^{2}\leq|\lambda_{K}|^{-2}tr(ST_{\mu\mu}^{-1}S)\max_{k}\|\chi_{k}^{*}\|^{2}$

where $\|\chi_{k}^{*}\|$ is finite for each $k$ since $\chi_{k}^{*}\in
L_{2}(\mathcal{X})$.

We can then apply the argument from the proof of Theorem 3.1 where we choose
$\tilde{\phi}_{k\varepsilon}(x)\equiv\phi_{k}^{*}(x)$, the solutions to (3.3)
corresponding to the $k$th largest eigenvalue for $k=1,\dots,K$. Since those
functions are chosen independently of $\varepsilon$, the bound
$B_{\varepsilon}$ is fixed at some finite value $B_{0}$, so the claim of the
proposition follows immediately from (B.7), noting that
$\kappa_{\varepsilon}\equiv\kappa_{0}$ is constant∎

###  B.3. Proof of Corollary 3.1

By construction, the optimal basis for the regularized problem 3.4 satisfies
the constraint
$\langle\phi_{k}^{*},(T_{\mu\mu}+a\textnormal{Id})\phi_{l}^{*}\rangle=\delta_{kl}$.
Hence, we can rewrite the $K\times K$ matrix
$P_{K}^{*}T_{\mu\mu}P_{K}^{*}=I_{K}-aP_{K}^{*}$, and use the Neumann series to
obtain its inverse,
$(P_{K}^{*}T_{\mu\mu}P_{K}^{*})^{-1}=I_{K}+\frac{a}{1-a}P_{K}^{*}$. From the
proof of Theorem B.1, the projection of $\tau_{g}$ onto the optimal basis is
therefore given by

$\tau_{gK}^{*}(x):=\sum_{k=1}^{K}\langle\mu_{g},\phi_{k}^{*}\rangle\frac{1+a}{1-a}\left(T_{\mu\tau}^{*}\phi_{k}^{*}\right)(x)$

which establishes the formula given in the Corollary∎

##  Appendix C Proofs for Section 4

###  C.1. Proof of Theorem 4.1

We use the main result in Graham, Niu, and Powell (2021) which adapts Theorems
2 and 10 in Hansen (2008) to nonparametric regression for dyadic data, rather
than conventional sample averages. We first apply their results separately for
each clusters $g=1,\dots,G$, where we strengthen the rate conditions to ensure
uniformity across clusters. We then aggregate the cluster-specific estimates
to obtain the first and second conditional moments across clusters.

####  C.1.1. Convergence Rate for Kernel Averages

The local linear estimator can be expressed in terms of weighted averages of
products of $X_{gi},d_{gi},Y_{gi}$. For a general notation, let
$W_{sgi}=(a_{0s}+a_{1s}X_{gi})(b_{0s}+b_{1s}Y_{gi})$ be a function of
$X_{gi},Y_{gi}$ that is affine in $Y_{gi}$ given fixed coefficients
$a_{0s},b_{0s},a_{1s},b_{1s}$. We consider uniform convergence of conventional
and dyadic kernel averages

$\displaystyle\hat{\Psi}_{g}(x_{1};d_{1})$ $\displaystyle:=$
$\displaystyle\frac{1}{nh^{d}}\sum_{i=1}^{n}W_{1gi}K\left(\frac{x_{1}-x_{gi}}{h}\right)1\textnormal{l}\\{d_{gi}=d_{1}\\}$
$\displaystyle\hat{\Omega}_{g}(x_{1},x_{2};d_{1},d_{2})$ $\displaystyle:=$
$\displaystyle\binom{n}{2}^{-1}h^{-2d}\frac{1}{2}\sum_{i\neq
j}W_{1gi}W_{2gj}K\left(\frac{x_{1}-x_{gi}}{h}\right)K\left(\frac{x_{2}-x_{gj}}{h}\right)1\textnormal{l}\\{d_{gi}=d_{1},d_{gj}=d_{2}\\}$

for each $g=1,\dots,G$ such that $R_{g}=1$.

######  Lemma C.1.

Suppose that the bandwidth $g$ satisfies $\left(\frac{n}{\log
n}\right)^{1/3}h^{d}\rightarrow\infty$ and $\delta_{Gn}^{-1}\left(\frac{\log
n}{nh^{d}}\right)^{1/2}$. Under Assumptions 2.4 and 4.1, the kernel averages
$\hat{\Psi}_{g}$ and $\hat{\Omega}_{g}$ converge uniformly to their respective
expectations at the rate

$\displaystyle\max_{g=1,\dots,G}R_{g}\sup_{\|x_{1}\|\leq
c_{n}}\left|\hat{\Psi}_{g}(x_{1};d_{1})-\mathbb{E}[\hat{\Psi}_{g}(x_{1};d_{1})]\right|$
$\displaystyle=$ $\displaystyle O_{P}\left(\delta_{Gn}^{-1}\left(\frac{\log
n}{nh^{d}}\right)^{1/2}\right)$
$\displaystyle\max_{g=1,\dots,G}R_{g}\sup_{\|x_{1}\|,\|x_{2}\|\leq
c_{n}}\left|\hat{\Omega}_{g}(x_{1},x_{2};d_{1},d_{2})-\mathbb{E}[\hat{\Omega}_{g}(x_{1},x_{2};d_{1},d_{2})]\right|$
$\displaystyle=$ $\displaystyle O_{P}\left(\delta_{Gn}^{-1}\left(\frac{\log
n}{nh^{d}}\right)^{1/2}\right)$

for any $d_{1},d_{2}\in\\{0,1\\}$.

Proof: We prove this result using Theorem 3.2 in Graham, Niu, and Powell
(2021). Note first that, since $W_{gi}$ are conditionally i.i.d., Assumption
2.4 implies Assumptions 3.1 and 3.3 (a) in Graham, Niu, and Powell (2021) for
$Z_{gij}:=W_{gi}W_{gj}$. Also, Assumption 4.1 subsumes Assumptions 3.2 and 3.3
(b) in Graham, Niu, and Powell (2021). Moreover, for $s\geq 3$, lengthy but
elementary rate calculations confirm that any bandwidth sequence with
$\left(\frac{n}{\log n}\right)^{1/3}h^{d}\rightarrow\infty$ satisfies the
additional bandwidth conditions required for their theorem.

Their argument can then be adapted to achieve uniform convergence of
$\hat{\Psi}_{gs}(x_{1};d_{1})$ and
$\hat{\Omega}_{gs}(x_{1},x_{2};d_{1},d_{2})$ with respect to $x_{1},x_{2}$ and
$g$. To that end, the grid $\mathcal{X}^{*}_{n}\\{w_{n1},\dots,w_{nL_{n}}\\}$
is chosen in a way such that the set $\left\\{x:\|x\|\leq c_{n}\right\\}$ is
covered by the collection of balls of radius $a_{n}h$. We then replace the
approximating grid introduced on p.19 of Graham, Niu, and Powell (2021) with
$\mathcal{X}_{n}^{*}\times\\{1,\dots,G\\}$, resulting in
$L_{n}:=G\left(c_{n}h^{-1}\left(\frac{\log n}{nh^{d}}\right)^{1/2}\right)^{d}$
partition elements. We can therefore conclude that

$\max_{g=1,\dots,G}R_{g}\sup_{\|x_{1}\|,\|x_{2}\|\leq
c_{n}}\left|\hat{\Omega}_{gs}(x_{1},x_{2};d_{1},d_{2})-\mathbb{E}[\hat{\Omega}_{g1}(x_{1},x_{2};d_{1},d_{2})]\right|=O_{p}\left(\frac{\log
n}{nh^{d}}\right)^{1/2}$ (C.2)

for $s=1,2$. The claim regarding $\hat{\Psi}(x_{1};d_{1})$ is proven in an
analogous manner using Theorem 2 in Hansen (2008), whose conditions are
subsumed under those for Theorem 3.2 in Graham, Niu, and Powell (2021)∎

####  C.1.2. Proof of Theorem 4.1

We now complete the proof of Theorem 4.1. We consider the general case of
estimating the conditional expectation

$\Psi_{g}(x;d)\equiv\Psi_{g}(x_{1},x_{2};d_{1},d_{2}):=\mathbb{E}\left[W_{1gi},W_{2gj}|X_{gi}=x,X_{gj}=x_{2},D_{gi}=d_{1},D_{gj}=d_{2}\right]$

and
$\Lambda(x_{1},x_{2};d_{1},d_{2}):=\frac{1}{G}\sum_{g=1}^{G}\Lambda_{g}(x_{1},x_{2};d_{1},d_{2}$.
for general $W_{1gi},W_{2gj}$. The corresponding cluster-specific local linear
estimator

$\displaystyle\hat{\Psi}_{g}(x_{1},x_{2};d_{1},d_{2})$ $\displaystyle:=$
$\displaystyle\arg_{b_{0}^{(g)}}\min_{b_{0}^{(g)},b_{11}^{(g)},b_{12}^{(g)}}\sum_{i=1}^{n_{g}}\sum_{j\neq
i}\Psi_{gij}(x_{1},x_{2};\mathbf{b})w_{gi}(x_{1};d_{1})w_{gj}(x_{2};d_{2})$

where we denote
$\Psi_{gij}(x_{1},x_{2};\mathbf{b}):=\left(W_{1gi}W_{2gj}-b_{0}^{(g)}-b_{11}^{(g)}(X_{gi}-x_{1})-b_{12}^{(g)}(X_{gj}-x_{2})\right)$,
and as before,
$w_{gi}(x_{1};d_{1}):=K\left(\frac{X_{gi}-x_{1}}{h}\right)1\textnormal{l}\\{d_{gi}=d_{1}\\}$.

As in the proof of Theorem 10 in Hansen (2008), we can write each local linear
estimator as

$\hat{\Lambda}_{g}(x;d)=\frac{\hat{m}_{g}+\hat{S}_{g}^{\prime}\hat{M}_{g}^{-1}\hat{N}_{g}}{\hat{f}_{g}-\hat{S}_{g}^{\prime}\hat{M}_{g}^{-1}\hat{S}_{g}}$

where for greater legibility, we write $(x,d):=(x_{1},x_{2};d_{1},d_{2})$ and
suppress dependence on $d_{1},d_{2}$ wherever possible, and terms on the
right-hand side are defined as follows:

$\displaystyle\hat{m}_{g}$ $\displaystyle=$
$\displaystyle\hat{m}_{g}(x;d):=\binom{n}{2}^{-1}h^{-2d}\frac{1}{2}\sum_{i\neq
j}W_{gi}W_{gj}w_{gi}(x_{1};d_{1})w_{gj}(x_{2};d_{2})$
$\displaystyle\hat{f}_{g}$ $\displaystyle=$
$\displaystyle\hat{f}_{g}(x;d):=\binom{n}{2}^{-1}h^{-2d}\frac{1}{2}\sum_{i\neq
j}w_{gi}(x_{1};d_{1})w_{gj}(x_{2};d_{2})$ $\displaystyle\hat{S}_{g}$
$\displaystyle=$
$\displaystyle\hat{S}_{g}(x;d):=\binom{n}{2}^{-1}h^{-2d}\frac{1}{2}\sum_{i\neq
j}\left(\frac{x-X_{gij}}{h}\right)w_{gi}(x_{1};d_{1})w_{gj}(x_{2};d_{2})$
$\displaystyle\hat{M}_{g}$ $\displaystyle=$
$\displaystyle\hat{M}_{g}(x;d):=\binom{n}{2}^{-1}h^{-2d}\frac{1}{2}\sum_{i\neq
j}\left(\frac{x-X_{gij}}{h}\right)\left(\frac{x-X_{gij}}{h}\right)^{\prime}w_{gi}(x_{1};d_{1})w_{gj}(x_{2};d_{2})$
$\displaystyle\hat{N}_{g}$ $\displaystyle=$
$\displaystyle\hat{N}_{g}(x;d):=\binom{n}{2}^{-1}h^{-2d}\frac{1}{2}\sum_{i\neq
j}\left(\frac{x-X_{gij}}{h}\right)W_{gij}w_{gi}(x_{1};d_{1})w_{gj}(x_{2};d_{2})$

Applying Lemma C.1 component by component, each of these terms converges
uniformly to its expectation. Specifically, denoting $b_{n}:=\left(\frac{\log
n}{nh^{d}}\right)^{1/2}$, and $\Sigma:=\int uu^{\prime}K(u)du$, standard
calculations for conditional moment estimation using local linear regression
(see also the proofs of Theorems 8 and 10 in Hansen (2008)) then yield

$\displaystyle\hat{m}_{g}(x;d)$ $\displaystyle=$ $\displaystyle
m_{g}(x;d)+O(h^{2})+O_{p}(b_{n})$ $\displaystyle\hat{f}_{g}(x;d)$
$\displaystyle=$ $\displaystyle f_{g}(x;d)+O(h^{2})+O_{p}(b_{n})$
$\displaystyle\hat{S}_{g}(x;d)$ $\displaystyle=$ $\displaystyle
h\Sigma\nabla_{x}f_{g}(x;d)+O(h^{2})+O_{p}(b_{n})$
$\displaystyle\hat{M}_{g}(x;d)$ $\displaystyle=$ $\displaystyle\Sigma
f_{g}(x;d)+O(h^{2})+O_{p}(b_{n})$ $\displaystyle\hat{N}_{g}(x;d)$
$\displaystyle=$ $\displaystyle
h\Sigma\nabla_{x}m_{g}(x;d)+O(h^{2})+O_{P}(b_{n})$

uniformly in $(x;d)$, where $m_{g}(x;d):=\Psi_{g}(x;d)f_{g}(x;d)$.

We can now confirm that $b_{n}^{1/3}h\equiv\left(\frac{n}{\log
n}\right)^{1/3}h^{d}\rightarrow\infty$ implies that

$\frac{h}{b_{n}}=\left(\frac{n}{\log
n}\right)h^{d+1}=\left[\left(\frac{n}{\log
n}\right)^{1/3}h^{d}\right]^{3}h^{1-2d}\rightarrow\infty$

for any $d\geq 1$. Moreover, by assumption,
$\delta_{Gn}/b_{n}\rightarrow\infty$. Since the log-derivative $\nabla_{x}\log
f_{g}(x)=\nabla_{x}f_{g}(x;d)/f_{g}(x;d)$ is uniformly bounded, so is
$\left(\nabla_{x}m_{g}(x;d)\right)/f_{g}(x;d)=\nabla_{x}\Psi_{g}(x;d)+\Psi_{g}(x;d)\nabla_{x}f_{g}(x;d)/f_{g}(x;d)$.
Hence, collecting terms,

$\displaystyle\left|\hat{S}_{g}(x;d)^{\prime}\hat{M}_{g}(x;d)^{-1}\hat{N}_{g}(x;d)-h^{2}\frac{\nabla_{x}f_{g}(x;d)^{\prime}\Sigma\nabla_{x}m_{g}(x;d)}{f_{g}(x;d)}\right|$
$\displaystyle=$ $\displaystyle
O_{p}\left(\left\\{\frac{h^{2}}{f_{g}(x;d)}+h\right\\}b_{n}\right)$
$\displaystyle\leq$ $\displaystyle
O_{p}\left(h^{2}\delta_{Gn}^{-1}b_{n}+hb_{n}\right)=hO_{p}\left(b_{n}\right)$

where convergence is uniform in $(x;d)$. Note also that
$h^{2}\frac{\nabla_{x}f_{g}(x;d)^{\prime}\Sigma\nabla_{x}m_{g}(x;d)}{f_{g}(x;d)}=O\left(\delta_{Gn}^{-1}h^{2}\right)$
which goes to zero more slowly than $O_{p}\left(\delta_{Gn}^{-1}b_{n}\right)$
since $h/b_{n}\rightarrow\infty$.

Hence, from standard rate calculations

$\displaystyle\hat{\Psi}_{g}(x;d)$ $\displaystyle=$
$\displaystyle\frac{\hat{m}_{g}+\hat{S}_{g}^{\prime}\hat{M}_{g}^{-1}\hat{N}_{g}}{\hat{f}_{g}-\hat{S}_{g}^{\prime}\hat{M}_{g}^{-1}\hat{S}_{g}}$
$\displaystyle=$
$\displaystyle\frac{\Psi_{g}(x;d)f_{g}(x;d)+O(h^{2})+O_{p}\left(b_{n}+hb_{n}\right)}{f_{g}(x;d)+O(h^{2})+O_{P}\left(b_{n}+hb_{n}\right)}$
$\displaystyle=$
$\displaystyle\Psi_{g}(x;d)+O_{P}\left(\delta_{Gn}^{-1}b_{n}\right)+O\left(\delta_{Gn}^{-1}h^{2}\right)$

Since units $i=1,\dots,n$ are sampled independently in each location, the
$O_{P}(\cdot)$ terms are independent across locations $g=1,\dots,G$ with
expectation of order $o(h^{2})$.

Hence, aggregating over $g=1,\dots,G$,

$\displaystyle\hat{\Psi}(x_{1},x_{2};d_{1},d_{2})$ $\displaystyle=$
$\displaystyle\frac{1}{G-1}\sum_{g=1}^{G}R_{g}\hat{\Psi}_{g}(x_{1},x_{2};d_{1},d_{2})$

Hence, by the triangle inequality

$\displaystyle\left|\hat{\Psi}(x_{1},x_{2};d_{1},d_{2})-\Psi(x_{1},x_{2};d_{1},d_{2})\right|$
$\displaystyle\leq$
$\displaystyle\left|\frac{1}{G-1}\sum_{g=1}^{G}R_{g}\hat{\Psi}_{g}(x_{1},x_{2};d_{1},d_{2})-\Psi_{g}(x_{1},x_{2};d_{1},d_{2})\right|$
$\displaystyle+\left|\frac{1}{G-1}\sum_{g=1}^{G}R_{g}|\Psi_{g}(x_{1},x_{2};d_{1},d_{2})-\Psi(x_{1},x_{2};d_{1},d_{2})\right|$
$\displaystyle=$ $\displaystyle
O_{p}\left(\frac{b_{n}+h^{2}}{G^{1/2}\delta{Gn}^{-1}}\right)+O_{P}\left(\frac{1}{G}\right)$

since $R_{g}$ is equal to zero for a single unit $g^{*}$ selected at random
and one otherwise, where we use unconfoundedness of location, Assumption 2.2
and bounded conditional moments in Assumption 2.4. By our previous arguments,
convergence is also uniform with respect to the arguments
$x_{1},x_{2};d_{1},d_{2}$.

We can immediately verify that $\hat{\mu}(x;d)$ and $\mu(x;d)$ correspond to
$\hat{\Psi}(x;d)$ and $\Psi(x;d)$, respectively, for the case $W_{1gi}=Y_{gi}$
and $W_{2gj}=1$, so that

$\sup_{\|x\|\leq
c_{n}}|\hat{\mu}(x;d)-\mu(x;d)|=O_{p}\left(\frac{b_{n}+h^{2}}{G^{1/2}\delta{Gn}^{-1}}\right)+O_{P}\left(\frac{1}{G}\right)$

For the covariance kernel $H(x_{1},x_{2};d_{1},d_{2})$, we can set
$W_{1gi}=Y_{gi}$ and $W_{2gj}=Y_{gj}$ so that the cluster-specific local
linear estimator

$\displaystyle\hat{\Psi}_{g}(x_{1},x_{2};d_{1},d_{2})$ $\displaystyle:=$
$\displaystyle\arg_{b_{0}^{(g)}}\min_{b_{0}^{(g)},b_{11}^{(g)},b_{12}^{(g)}}\sum_{i=1}^{n_{g}}\sum_{j\neq
i}\left\\{1\textnormal{l}\\{D_{gi}=d_{1},D_{gj}=d_{2}\\}\left(Y_{gi}Y_{gj}-b_{0}^{(g)}-b_{11}^{(g)}(X_{gi}-x_{1})-b_{12}^{(g)}(X_{gj}-x_{2})\right)\right\\}$
$\displaystyle\times w_{gi}(x_{1};d_{1})w_{gj}(x_{2};d_{2})$

is uniformly consistent for any $g=1,\dots,G$ with $R_{g}=1$ so that for

$\hat{H}(x_{1},x_{2};d_{1},d_{2})=\frac{1}{G-1}\sum_{g=1}^{G}R_{g}\hat{\Psi}_{d_{1}d_{2},g}(x_{1},x_{2};d_{1},d_{2})-\bar{\mu}(x_{1};d_{1})\bar{\mu}(x_{2};d_{2})$

we can conclude

$\sup_{\|x_{1}\|,\|x_{2}\|\leq
c_{n}}|\hat{H}(x_{1},x_{2};d_{1},d_{2})-H(x_{1},x_{2};d_{1},d_{2})|=O_{p}\left(\frac{b_{n}+h^{2}}{G^{1/2}\delta{Gn}^{-1}}\right)+O_{P}\left(\frac{1}{G}\right)$

establishing the convergence rates for a general choice of the bandwidth
sequence subject to the rate conditions in the theorem.

Since by standard arguments the bias is of the order $h^{2}$, the rate of the
root mean square error is minimized at bandwidth sequences such that
$h^{2}=\left(\frac{\log n}{Gnh^{d}}\right)^{1/2}$ so that such a sequence must
go to zero at a rate $h^{*}=O\left(\frac{\log n}{Gn}\right)^{\frac{1}{4+d}}$.∎

###  C.2. Proof of Theorem 4.2

For the following formal arguments we let $\|\cdot\|_{F}$ denote the trace
operator norm $\|T\|_{F}:=\textnormal{tr}(T^{*}T)$. We start by proving the
following Lemma:

######  Lemma C.2.

Suppose that $S$ is a compact, self-adjoint operator with eigenvalues
$\lambda_{1}\geq\lambda_{2},\dots,$ counted by their multiplicity, and
corresponding eigenfunctions $\phi_{1},\phi_{2},\dots$. Then for any sequence
$\hat{S}$ of compact, self-adjoint operators with eigenvalues
$\hat{\lambda}_{1}\geq\hat{\lambda}_{2},\dots$ such that
$\|\hat{S}-S\|_{F}=O_{P}(r_{Gn})$ we can choose the eigenfunctions
$\hat{\phi}_{1},\hat{\phi}_{2},\dots$ such that

$\displaystyle|\hat{\lambda}_{k}-\lambda_{k}|$ $\displaystyle=$ $\displaystyle
O_{p}\left(r_{Gn}\right)$ $\displaystyle\|\hat{\phi}_{k}-\phi_{k}\|$
$\displaystyle=$ $\displaystyle O_{p}\left(r_{n}\right)$ (C.3)

for each $k=1,2,\dots$ such that $\lambda_{k}>0$.

Proof: We follow closely the proofs for Proposition 4.2 in Gobet, Hoffmann,
and Reiß (2004) and Theorem 2 in Yao, Müller, and Wang (2005a). Since $S$ is
compact, existence of eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\dots$ and
eigenfunctions $\phi_{1},\phi_{2},\dots$ follows from Mercer’s Theorem (see
Happ and Greven (2018) for a multivariate generalization where $\mathcal{X}$
may be of dimension greater than one).

We define the resolvent maps of the operators $S$ and $\hat{S}$,