# causal_effect_of_functional_treatment__c267a3ed.pdf

Journal of Machine Learning Research 26 (2025) 1-39 Submitted 3/23; Revised 7/24; Published 5/25

Causal Effect of Functional Treatment

Ruoxu Tan ruoxut@tongji.edu.cn School of Mathematical Sciences and School of Economics and Management Tongji University Shanghai, China

Wei Huang wei.huang@unimelb.edu.au School of Mathematics and Statistics University of Melbourne Melbourne, VIC, Australia

Zheng Zhang zhengzhang@ruc.edu.cn Center for Applied Statistics, Institute of Statistics & Big Data Renmin University of China Beijing, China

Guosheng Yin gyin@hku.hk Department of Statistics and Actuarial Science University of Hong Kong Hong Kong SAR, China

Editor: Eric Laber

We study the causal effect with a functional treatment variable, where practical applications often arise in neuroscience, biomedical sciences, etc. Previous research concerning the effect of a functional variable on an outcome is typically restricted to exploring correlation rather than causality. The generalized propensity score, which is often used to calibrate the selection bias, is not directly applicable to a functional treatment variable due to a lack of definition of probability density function for functional data. We propose three estimators for the average dose-response functional based on the functional linear model, namely, the functional stabilized weight estimator, the outcome regression estimator and the doubly robust estimator, each of which has its own merits. We study their theoretical properties, which are corroborated through extensive numerical experiments. A real data application on electroencephalography data and disease severity demonstrates the practical value of our methods. Keywords: Causality, Double Robustness, Functional Data, Functional Linear Model, Treatment Effect

1. Introduction

In observational studies, a fundamental question is to investigate the causal effect of a treatment variable on an outcome. Much of the literature on estimating causal effects from observational data focus on a binary treatment variable (i.e., treatment versus control); see e.g., Rosenbaum and Rubin (1983); Hirano et al. (2003); Imai and Ratkovic (2014); Chan

. Wei Huang and Zheng Zhang are corresponding authors.

c 2025 Ruoxu Tan, Wei Huang, Zheng Zhang and Guosheng Yin.

License: CC-BY 4.0, see https://creativecommons.org/licenses/by/4.0/. Attribution requirements are provided at http://jmlr.org/papers/v26/23-0381.html.

Tan, Huang, Zhang and Yin

et al. (2016); Athey et al. (2018); Ding et al. (2019); Tan (2020); Guo et al. (2021); Lin et al. (2023). Recently, there is a growing interest in studying more complex treatment variables, e.g., categorical (Schulte et al., 2014; Lopez and Gutman, 2017a,b; Chen et al., 2018; Li and Li, 2019; Luckett et al., 2019) or continuous treatment variables (Hirano and Imbens, 2004; Moodie et al., 2014; Galvao and Wang, 2015; Laber and Zhao, 2015; Kennedy et al., 2017; Fong et al., 2018; Kennedy, 2019; Ai et al., 2021; Bonvini and Kennedy, 2022). Our interest is in estimating the average dose-response functional (ADRF) of a functional treatment variable. Functional data are collected repeatedly over a continuous domain, and are fundamentally different from scalar variables due to their infinite dimensionality (Wang et al., 2016). We focus on the causal relationship between a functional treatment variable and a scalar outcome variable in cross-sectional studies, where the treatment variable is real-valued over a continuous domain (e.g., space or frequency, rather than time) and the outcome variable does not vary with respect to that continuous domain. In particular, our work is motivated by an electroencephalography (EEG) dataset (Ciarleglio et al., 2022), which was collected via a standard headcap with multiple scalp electrodes. The main purpose there is to assess the EEG data as a potential moderator of treatment effect in a clinical trial. Here, we investigate possible causation between the EEG data and the severity of major depressive disorder. A random subsample of 20 frontal asymmetric curves, measuring the difference between the intensities of the human neuronal activities from two electrodes (right and left) placed on the scalp, is shown in Figure 1. For additional applications, see Morris (2015). As an illustration, we may be interested in the causal relationship between spectrometric data and fat content of certain types of food (Ferraty and Vieu, 2006) and the causal effect of body shape represented by the circumference over a given range on human visceral adipose tissue (Zhang et al., 2021). These problems are distinct from sequential decision problems in which multiple treatments are assigned longitudinally (Robins, 2000; Robins et al., 2000; Moodie and Stephens, 2012; Zhao et al., 2015; Kennedy, 2019). In addition, functional data have been used as one of the covariates in estimating optimal treatment regimes (Mc Keague and Qian, 2014; Ciarleglio et al., 2015, 2016, 2018; Laber and Staicu, 2018; Li and Kosorok, 2023; Park et al., 2023). To identify the causal effect in an observational study, certain assumptions are required. In addition to standard assumptions such as no interference, consistency and positivity, a routinely used condition for identification is the unconfoundedness assumption (Rosenbaum and Rubin, 1983), i.e., the set of potential outcomes is conditionally independent of the treatment given the observed confounders. For a continuous scalar treatment, under the unconfoundedness assumption, the generalized propensity score or the stabilized weight (Hirano and Imbens, 2004; Imai and van Dyk, 2004; Galvao and Wang, 2015; Kennedy et al., 2017; Ai et al., 2021; Bonvini and Kennedy, 2022) has been widely used to calibrate the selection bias via estimating equations. It relies upon the (conditional) density of the treatment variable. However, the probability density function of functional data generally does not exist (Delaigle and Hall, 2010). Therefore, in our context, neither the generalized propensity score nor the stabilized weight are well defined, and thus traditional methods are not directly applicable. To circumvent this issue, Zhang et al. (2021) proposed to approximate the functional treatment variable based on functional principal component (PC) analysis. They then defined a stabilized weight function by the (conditional) probability density of the PC scores to adjust for the selection bias. However, as the number of PCs

Functional Treatment Effect

5 10 15 20 25 30

Figure 1: A random subsample of 20 frontal asymmetric curves Z( ) on a given frequency domain from the EEG dataset.

approaches infinity, this stabilized weight becomes ill-defined, and the consistency of their estimator is not guaranteed. Under the unconfoundedness assumption in the context of a functional treatment variable, we propose three new estimation approaches to the ADRF: a functional stabilized weight (FSW) method, a partially linear outcome regression (OR) method, and a doubly robust (DR) approach combining the first two. The first method assumes a functional linear model for the ADRF and incorporates a new FSW. Our FSW is well-defined and guarantees the identification of the ADRF for any functional treatment variable. We propose a novel, consistent, nonparametric FSW estimator and estimate the ADRF using functional linear regression with a weighted outcome. The OR method assumes a more restrictive additive model that considers confounding variables. It quantifies the ADRF by the expectation of the parametric model with respect only to the confounding variables, which is computed through a backfitting algorithm. Our DR estimator is consistent if either of the first two estimators is consistent. If the partially linear model is correctly specified, the DR estimator converges at the same fast rate as the OR estimator. Our double robustness is different from those in the literature on a continuous scalar treatment variable (see e.g., Kennedy et al., 2017; Bonvini and Kennedy, 2022), where the estimator is consistent if one of two models for the conditional outcome and the generalized propensity score is correctly specified. In our case, the functional linear model for the ADRF is required for the consistency of all three estimators, while the FSW, serving the role of the generalized propensity score, is estimated nonparametrically. More recently, Wang et al. (2023) proposed a so-called weight-modified kernel ridge regression (WMKRR) estimator for the ADRF. They assumed that the ADRF lies in a functional reproducing kernel Hilbert space (RKHS) and utilized the kernel ridge regression to estimate the ADRF. Their method also requires an estimation of the FSW. To this end, they proposed to minimize a novel uniform balancing error derived from their estimation of the ADRF. Although the kernel ridge regression is more flexible than the functional linear model that is used in our methods, it still faces the risk of model misspecification since the

Tan, Huang, Zhang and Yin

model depends on the choice of the kernel. We will compare our methods with those from Zhang et al. (2021) and Wang et al. (2023) in simulation experiments. The rest of the paper is organized as follows. We introduce the basic setup for the identification of ADRF in Section 2, followed by the development of our three estimators in Section 3. We investigate the convergence rates of the estimators in Section 4. The selection of tuning parameters is discussed in Section 5. Simulation experiments and the real data analysis are presented in Sections 6 and 7, respectively. We conclude in Section 8. Proofs and additional technical details are provided in the Appendix.

2. Identification of the Average Dose-Response Functional

We consider a functional treatment variable such that the observed treatment variable, denoted by Z = Z( ) : T R, is a smooth real-valued random function defined on a compact interval T , with E{ R

T Z2(t) dt} < . Let Y (z) be the potential outcome given the treatment Z = z for z L2(T ), the Hilbert space of squared integrable functions on T . Let X = (X1, . . . , Xp) Rp be the p-dimensional observable confounding covariates that are related to both Y (z) and Z. For each individual, we only observe a confounding vector X, a random treatment Z and the corresponding outcome Y = Y (Z) R. We further assume that the functional variable Z is fully observed. Essentially, our methodology can be applied when Z is only sparsely observed as long as the trajectory of Z can be well reconstructed, while the development of asymptotic properties in such cases would be much more challenging (Zhang and Wang, 2016; Zhou et al., 2023). We are interested in estimating the average dose-response functional (ADRF), E{Y (z)}, for any z L2(T ). Because we never observe Y (z) simultaneously for all z L2(T ), to identify E{Y (z)} from the observed data, we impose the following identification assumptions.

Assumption 1

(i) Unconfoundedness: Given X, Z is independent of Y (z) for all z L2(T ), i.e., {Y (z)}z L2(T ) Z | X.

(ii) No interference: There is no interference among the units, i.e., each individual s outcome depends only on their own treatment.

(iii) Consistency: Y = Y (z) a.s. if Z = z.

(iv) Positivity: The conditional density f X|Z satisfies f X|Z(X|Z) > 0 a.s.

Assumption 1 is a natural extension of the classical identification assumption in the literature on the scalar treatment effect to the context of a functional treatment variable. In particular, we focus on cross-sectional studies where the interval T does not represent a time interval and there is no dynamic confounding effect. To circumvent the problem of the nonexistent probability density function of functional data, the positivity condition (iv) is imposed on f X|Z instead of f Z|X.

Functional Treatment Effect

Under Assumption 1, the ADRF can be identified from the observable variables (X, Y, Z) in the three ways as follows. First, it is easy to verify

E{Y (z)} = E[E{Y (z)|X}] = E[E{Y (z)|X, Z = z}] = E{E(Y |X, Z = z)}. (1)

Second, we can further write

E{Y (z)} = E{E(Y |X, Z = z)}

Rp E(Y |X = x, Z = z)f X(x) dx

Rp E(Y |X = x, Z = z) f X(x) f X|Z(x|z)f X|Z(x|z) dx

= E{π0(Z, X)Y |Z = z} , (2)

where π0(Z, X) = f X(X)/f X|Z(X|Z) is called the functional stabilized weight (FSW). If Z is a scalar random variable, we have f X(X)/f X|Z(X|Z) = f Z(Z)/f Z|X(Z|X), where the right-hand side is the classical stabilized weight. Third, let EX denote the expectation with respect to X. Noting that

E[EX{E(Y |X, Z)}|Z = z] = E{E(Y |X, Z = z)} ,

E{E(Y |X, Z)π0(X, Z)|Z = z} = E{Y π0(X, Z)|Z = z} ,

together with (1) and (2), we can also identify E{Y (z)} as

E{Y (z)} = E[{Y E(Y |X, Z)}π0(X, Z) + EX{E(Y |X, Z)}|Z = z] , (3)

because the second and third terms above cancel out. Based on the three identification strategies in (1), (2), and (3) respectively, we develop three estimators of E{Y (z)}. Each estimator is based on different model assumptions and has its own merits.

3. Models and Estimation Methods

For all three methods, we adopt a functional linear model on the ADRF,

E{Y (z)} = a + Z

T b(t)z(t) dt , (4)

where a is a scalar intercept, and b L2(T ) is a slope function. The functional linear model is widely used in scalar-on-function regression, whereas the model here is assumed for the potential outcome Y (z). Essentially, model (4) can be replaced by any other nonlinear regression models (e.g., M uller and Yao, 2008; Li et al., 2010; Wang et al., 2023) without affecting our main idea. We choose to use the functional linear model mainly due to its simplicity and interpretability. The slope function b = b( ) is a quantity of primary interest convenient for the causal inference because it measures the intensity of the functional treatment effect. For example, to compute the average treatment effect (ATE) between z1, z2 L2(T ), E{Y (z1) Y (z2)}, it suffices to compute R

T b(t){z1(t) z2(t)} dt. Additionally, the slope function provides useful causal information (e.g., a positive/negative causal relationship), even if the model is misspecified.

Tan, Huang, Zhang and Yin

3.1 Functional Stabilized Weight Estimator

We construct our functional stabilized weight (FSW) estimator of the ADRF based on equation (2) and model (4). In this construction, we propose a nonparametric estimator for the FSW. In particular, equation (2) suggests that we can estimate E{Y (z)} in (4) using the classical functional linear regression technique (e.g., Hall and Horowitz, 2007). Here, the response variable Y is replaced by bπ(Z, X)Y , provided an estimator of π0, bπ. To facilitate the development, suppose that an estimator bπ of π0 is available for the time being. Let G(t, s) = Cov{Z(t), Z(s)} be the covariance function of Z. Its eigenvalues λ1 λ2 0 and the corresponding orthonormal eigenfunctions {φj} j=1 are defined as follows, Z

T G(t, s)φj(s) ds = λjφj(t) , for all t T , (5)

where {φj} j=1, also called the PC basis, is a complete basis of L2(T ). According to the Karhunen Lo eve expansion, we have

Z(t) = µZ(t) +

j=1 ξjφj(t) , (6)

where µZ(t) = E{Z(t)} and ξj = R

T {Z(t) µZ(t)}φj(t) dt. With a slight abuse of notation, let G : L2(T ) L2(T ) denote a linear operator such that (G ψ)(t) = R

T G(t, s)ψ(s) ds, for a function ψ L2(T ). If we let µY,π = E{π0(Z, X)Y } = a + R

T b(t)µZ(t) dt, by combining (2), (4) and (6), we can write

E{π0(Z, X)Y |Z = z} µY,π = Z

T b(t){z(t) µZ(t)} dt .

It follows that (G b)(t) = e(t), where e(t) = E[{π0(Z, X)Y µY,π}{Z(t) µZ(t)}]. Noting that {φj} j=1 is an orthonormal complete basis of L2(T ), we can write b(t) = P j=1 bjφj(t) and e(t) = P j=1 ejφj(t), and combining them with (5), we obtain that, for all j,

T b(t)φj(t) dt , and ej = Z

T e(t)φj(t) dt = Z

T (G b)(t)φj(t) dt = bjλj .

The above arguments suggest we can estimate b(t) by bb FSW(t) = Pq j=1bb FSW,j bφj(t),

where q is a truncation parameter and bb FSW,j = bλ 1 j be FSW,j is an estimator of bj. The bλj s and bφj s are the eigenvalues and eigenfunctions of the empirical covariance function b G(s, t) = Pn i=1{Zi(s) bµZ(s)}{Zi(t) bµZ(t)}/n, respectively, where bµZ(t) = Pn i=1 Zi(t)/n is the empirical estimator of µZ(t). Moreover, by the definition of e(t), we can estimate it by

be FSW(t) =

{Yibπ(Zi, Xi) bµY,π}{Zi(t) bµZ(t)}

with bµY,π = Pn i=1 Yibπ(Zi, Xi)/n, and take be FSW,j = R

T be FSW(t)bφj(t) dt as an estimator of ej. Finally, the estimator of the intercept a is given by ba FSW = bµY,π R

T bb FSW(t)bµZ(t) dt.

Functional Treatment Effect

It remains to estimate the FSW, π0(Z, X) = f X(X)/f X|Z(X|Z). A naive approach would first estimate the densities f X and f X|Z and then take the ratio. However, this may lead to unstable estimates because an estimator of the denominator f X|Z is difficult to derive and may be too close to zero. Also, approximating Z by its PC scores would not alleviate the challenge but pose additional difficulty in theoretical justification. To circumvent these problems, we treat π0 as a whole and develop a robust nonparametric estimator. We can find the conditional moments that identify π0(z, ) for every fixed z,

E{π0(z, X)v(X)|Z = z} = Z

Rp f X(x) f X|Z(x|z)v(x)f X|Z(x|z)dx = E{v(X)}, (7)

which holds for any integrable function v(X).

Proposition 1 For any fixed z L2(T ),

E{π(z, X)v(X)|Z = z} = E{v(X)} ,

holds for all integrable functions v(X) if and only if π(z, X) = π0(z, X) a.s.

The proof is given in Appendix A. This result indicates that π0 is fully characterized by (7), and thus we can use it to define our estimator. Since (7) is defined for a fixed z, we need to estimate π0 at all sample points Zi, where (Xi, Yi, Zi), for i = 1, . . . , n, are i.i.d. copies of (X, Y, Z). We consider the leave-one-out index set S i = {j : j = i, j = 1, . . . , n} to estimate π0 at each Zi. Specifically, for a given sample point z = Zi, the right-hand side of (7) can be estimated by its empirical version P j S j v(Xj)/(n 1). We propose to estimate the left-hand side of (7) by the Nadaraya Watson estimator for a functional covariate (Ferraty and Vieu, 2006), P j S i π0(Zi, Xj)v(Xj)K{d(Zj, Zi)/h}

P j S i K{d(Zj, Zi)/h} = 1 n 1

j S j v(Xj),

where d( , ) denotes a semi-metric in L2(T ), h > 0 is a bandwidth and K is a kernel function quantifying the proximity of Zj and Zi. Commonly used choices for d include the L2 norm and the projection metric d(z1, z2) = [ R

T {Π(z1) Π(z2)}2 dt]1/2, where Π denotes a projection operation to a subspace of L2(T ). As it is not possible to solve (7) for an infinite number of v s from a given finite sample, we use a growing number of basis functions to approximate any suitable function v. Specifically, let νk(X) = v1(X), . . . , vk(X) be a set of known basis functions, which may serve as a finite-dimensional sieve approximation to the original function space of all integrable functions. We define our estimator of π0(Zi, Xj), for j S i, as the solution to the k equations, P j S i π(Zi, Xj)νk(Xj)K{d(Zj, Zi)/h}

P j S i K{d(Zj, Zi)/h} = 1 n 1

j S i νk(Xj) , (8)

which asymptotically identifies π0 as k and h 0. However, in practice, a finite number k of equations in (8) cannot identify a unique solution. Indeed, for any strictly

Tan, Huang, Zhang and Yin

increasing and concave function ρ, replacing π(Zi, Xj) by ρ {bη Ziνk(Xj)} would satisfy (8). Here, ρ ( ) is the first derivative of ρ( ), and bηZi is a k-dimensional vector that maximizes the strictly concave function Hhk,Zi, defined as

Hhk,Zi(η) =

P j S i ρ{η νk(Xj)}K{d(Zj, Zi)/h}

P j S i K{d(Zj, Zi)/h} η 1 n 1

j S i νk(Xj) . (9)

This can be verified by noting that the gradient Hhk,Zi(bηZi) = 0 corresponds to (8) with π(Zi, Xj) replaced by ρ {bη Ziνk(Xj)}. Therefore, for a given strictly increasing and concave function ρ together with bηZi defined above, we define our estimator of π0(Zi, Xj) as bπhk(Zi, Xj) = ρ {bη Ziνk(Xj)}. This also induces the estimator bπhk(Zi, Xi) = ρ {bη Ziνk(Xi)} with Xj replaced by Xi. Our estimator bπhk has a generalized empirical likelihood interpretation for a variety of choices for ρ. Specifically, we show in Appendix B that the estimator defined by (9) is the dual solution to minimizing some distance measure between the FSW and 1 locally for each Zi, subject to the constraint (8) (i.e., the constraint (12) in Appendix B). That is, our estimator bπhk is the desired weight closest to the baseline uniform weight 1 locally in a small neighbourhood of each Zi, subject to the finite sample version of the moment restriction in Proposition 1. The uniform weight can be considered as a baseline. This is because, when Z and Y (z) are unconditionally independent for all z L2(T ) (i.e., a randomized trial without confoundedness), the functional stabilized weights should be uniformly equal to 1. Different choices of ρ correspond to different distance measures (see Appendix B for other examples). We suggest to choose ρ(v) = exp( v 1) which guarantees a positive weight estimator bπhk(Zi, Xi).

3.2 Outcome Regression Estimator and Doubly Robust Estimator

To estimate the ADRF based on the outcome regression (OR) identification method in (1), an estimator of E(Y |X, Z) is required. Thus, based on the functional linear model in (4), we further assume a partially linear additive model,

T b(t)Z(t) dt + g(X; θ) + ϵ , (10)

where a and b( ) are the same as those in (4), g( ; θ) is a known function with an unknown parameter θ Rp and ϵ is an error variable with E(ϵ|Z, X) = 0 and E(ϵ2|Z, X) < uniformly for all Z and X. Note that X and Z are generally dependent. We mainly investigate the linear function g(X; θ) = θ X, where X may involve some transformations, e.g., the log-transformation, polynomials, and exponential of some elements of X. Without loss of generality, we assume that E(θ X) = 0; otherwise the non-zero mean can be absorbed into a. Model (10) may be extended to include interactions between Z and X. However, research on functional models involving interactions between functional and multivariate continuous variables is scarce, even in the literature on non-causal functional regression. Such extension is beyond the scope of this paper and deserves future research. Recalling E{Y (z)} = E{E(Y |X, Z = z)} from (1), it is clear that (10) implies (4), and thus model (10) imposes a stronger structure than model (4). Specifically, the linear

Functional Treatment Effect

part θ X can be replaced by any nonparametric function of X whose mean is zero, which still implies model (4). However, model (10) gives a simpler estimator of E{Y (z)} with a faster convergence rate and a better interpretability of the effects of the treatment and confounding variables on the outcome. The OR estimators ba OR, bb OR and bθOR can be obtained by adapting the backfitting algorithm (Buja et al., 1989). In the first step, we set θ to be zero and apply the same method as in Section 3.1 to regress Y on Z, which gives the initial estimators of a and b in (10). Specifically, we estimate b(t) initially by bbini(t) = Pq j=1bbini,j bφj(t) with bbini,j = bλ 1 j beini,j, where bλj s and bφj s are the same as those in Section 3.1, while beini(t) = Pn i=1(Yi

bµY ){Zi(t) bµZ(t)}/n with bµY = Pn i=1 Yi/n. We then define beini,j = R

T beini(t)bφj(t) dt. Finally, the initial estimator baini of a is given by bµY R

T bbini(t)bµZ(t) dt. In the second step, we use the traditional least square method to regress the residual Yi baini R

T bbini(t)Zi(t) dt on Xi, for i = 1, . . . , n. Specifically, we estimate θ by bθini = (Ξ Ξ) 1Ξ Yres, where Ξ is an n p matrix with the (i, j)-th entry being Xij, the j-th element of Xi, and

Yres = n Y1 baini Z

T bbini(t)Z1(t) dt, . . . , Yn baini Z

T bbini(t)Zn(t) dt o .

Subsequently, we can repeat the first step with Yi replaced by Yi bθ ini Xi to update the estimators of baini and bbini. This procedure is iterated until the outcome regression (OR) estimators of a, b and θ, denoted by ba OR, bb OR and bθOR, are stabilized. More specifically, let ba(t) OR,bb(t) OR and bθ(t) OR denote respectively the estimators of a, b and θ at the t-th iteration. We set our final estimators as ba OR = ba(t+1) OR ,bb OR = bb(t+1) OR and bθOR = bθ(t+1) OR if {|ba(t+1) OR ba(t) OR|2 + bb(t+1) OR bb(t) OR 2 L2 + bθ(t+1) OR bθ(t) OR 2 Rp}1/2/{|ba(t) OR|2 + bb(t) OR 2 L2 + bθ(t) OR 2 Rp}1/2 < 0.05 or t reaches 100 (a prespecfied maximum number of iterations). Here, L2 and Rp denote the L2 norm in L2(T ) and the Euclidean norm in Rp, respectively. The OR estimation is easy to implement but requires strong parametric assumptions, while the FSW estimator requires fewer modelling assumptions but is subject to a slow convergence rate. Thus, it is desirable to develop an estimator that possesses more attractive properties by combining these two. Note that (3) satisfies the so-called doubly robust (DR) property: for two generic functions e E(Y |X, Z) and eπ(X, Z), it holds that

E{Y (z)} = E[{Y e E(Y |X, Z)}eπ(X, Z) + EX{ e E(Y |X, Z)}|Z = z] , (11)

if either e E(Y |X, Z) = E(Y |X, Z) or eπ(X, Z) = π0(X, Z) but not necessarily both are satisfied; see Appendix C for the derivation. Under model (4), (11) equals a + R

T b(t)z(t) dt. To estimate a and b, it suffices to conduct the functional linear regression as in Section 3.1 to regress an estimator of {Y E(Y |X, Z)}π(X, Z) + EX{E(Y |X, Z)} on Z. Specifically, we estimate E(Y |X, Z) by the OR estimator, b E(Y |X, Z) = ba OR + R

T bb OR(t)Z(t) dt + bθ ORX, and π0(X, Z) by bπhk(X, Z). We define the DR estimator of b as bb DR(t) = Pq j=1bb DR,j bφj(t), where bb DR,j = bλ 1 j be DR,j with

be DR,j = R

T be DR(t)bφj(t) dt. The expression of be DR(t) is

{Yi b E(Y |Xi, Zi)}bπhk(Xi, Zi) + 1

j=1 { b E(Y |Xj, Zi)} bµY,DR

{Zi(t) bµZ(t)} ,

Tan, Huang, Zhang and Yin

{Yi b E(Y |Xi, Zi)}bπhk(Xi, Zi) + 1

j=1 { b E(Y |Xj, Zi)} .

The DR estimator of a is defined as ba DR = bµY,DR R

T bb DRbµZ(t) dt. Provided that model (4) is correctly specified, according to the DR property (11), the DR estimator b EDR{Y (z)} = ba DR + R

T bb DR(t)z(t) dt is consistent if either b E(Y |X, Z) or bπhk(X, Z) is consistent. The consistency of b E(Y |X, Z) mainly depends on the correct specification of the partially linear model (10), while the consistency of bπhk(X, Z), as a nonparametric estimator, relies on less restrictive (but more technical) assumptions.

4. Asymptotic Properties

We investigate the asymptotic convergence rates of the proposed three estimators, while assuming that Z is fully observed. Under the main model is E{Y (z)} = a + R

T b(t)z(t) dt in (4), although the estimators of a are based on those of b, they are essentially empirical means, which can achieve the n 1/2 convergence rate (Shin, 2009). Thus, our primary focus is to derive the convergence rates of our estimators of b, which, depending on the smoothness of Z and b, are slower than those in the finite-dimensional settings. Let X Rp denote the support of X and Z L2(T ) the support of Z. To derive the convergence rate of bb FSW, we first provide the uniform convergence rate of the estimator bπhk of the FSW π0 over X and Z, which requires the following conditions.

Condition A

(A1) The set X is compact. For all z Z and x X, π0(z, x) is strictly bounded away from zero and infinity.

(A2) For all z Z, ηz Rk and a constant α > 0 such that sup(z,x) Z X |ρ ( 1){π0(z, x)} η z νk(x)| = O(k α), where ρ ( 1) is the inverse function of ρ .

(A3) The eigenvalues of E{νk(X)νk(X) |Z = z} are bounded away from zero and infinity uniformly in k and z. There exists a sequence of constants ζ(k) satisfying supx X νk(x) ζ(k).

(A4) There exists a continuously differentiable function µ(h) > 0, for h > 0, such that for all z Z, P{d(Z, z) h}/µ(h) is strictly bounded away from zero and infinity. There exists h0 > 0 such that suph (0,h0) µ (h) is bounded.

(A5) The kernel K : R R+ is bounded, Lipschitz and supported on [0, 1] with K(1) > 0.

(A6) For gk1,k2,z(X) = ρ {η z νk(X)}vk1(X) and vk1(X)vk2(X) with k1, k2 = 1, . . . , k, there exists λ > 0 such that for all m 2, E[gk1,k2,z1(X)|Z = z1] E[gk1,k2,z2(X)|Z = z2] /dλ(z1, z2) and E[|gk1,k2,z1(X)|m|Z = z1] are uniformly bounded over k, z1 and z2.

(A7) Let ψZ(ϵ) be the Kolmogorov ϵ-entropy of Z, i.e., ψZ(ϵ) = log{Nϵ(Z)}, where Nϵ(Z) is the minimal number of open balls in L2(T ) of radius ϵ (with the semi-metric d)

Functional Treatment Effect

covering Z. It satisfies P n=1 exp[ (δ 1)ψZ{(log n)/n}] < for some δ > 1, and (log n)2/{nµ(h)} < ψZ{(log n)/n} < nµ(h)/ log n for n large enough.

Condition (A1) requires that X is compactly supported, which can be relaxed by restricting the tail behaviour of the distribution of X at the cost of more tedious technical arguments. The boundedness of π0 in Condition (A1) (or some equivalent condition) is also commonly required in the literature (Kennedy et al., 2017; Ai et al., 2021; D Amour et al., 2021). This condition can be relaxed if other smoothness assumptions are made on the potential outcome distribution (Ma and Wang, 2020). Condition (A2) essentially assumes that the convergence rate of the sieve approximation η z νk(x) for ρ ( 1){π0(z, x)} is polynomial. This can be satisfied, for example, with α = + if X is discrete, and with α = s/r if X has r continuous components and νk is a power series, where s is the degree of smoothness of ρ ( 1){π0(z, )} with respect to the continuous components in X for all z L2(T ) (Chen, 2007). Further, Condition (A3) ensures that the sieve approximation conditional on a functional variable is not degenerate. Similar conditions are routinely assumed in the literature of sieve approximation (Newey, 1997). Conditions (A4) to (A6) are standard in functional nonparametric regression (Ferraty and Vieu, 2006; Ferraty et al., 2010). In particular, the function µ(h) in Condition (A4), also referred to as the small ball probability, has been studied extensively in the literature (Li and Shao, 2001; Ferraty and Vieu, 2006); Condition (A5) requires that K is compactly supported, which is not satisfied for the Gaussian kernel but convenient for technical arguments; the variable gk1,k2,z(X) in Condition (A6) serves as the role of the response variable in the traditional regression setting. Condition (A7) is less standard: it regularizes the support Z by controlling its Kolmogorov entropy, which is used to establish the uniform convergence rate on Z. This condition is satisfied, for example, if Z is a compact set, d( , ) is a projection metric and (log n)2 = O{nµ(h)}; see Ferraty et al. (2010) for other examples.

Theorem 1 Under Conditions (A1) to (A7), assuming that k and h 0 as n , we have

sup (z,x) Z X |bπhk(z, x) π0(z, x)| = O ζ(k) k α + hλ

ψZ{(log n)/n}k

holds almost surely.

Theorem 1 provides the sup-norm convergence rate of bπhk over the support of X and Z. The proof is given in Appendix D. The term ζ(k)(k α + hλ

k) and the term ζ(k) p

ψZ{(log n)/n}k/{nµ(h)} can be viewed as the convergence rates of the bias and standard deviation, respectively. To provide the convergence rate of bb FSW, we recall that E{π0(Z, X)Y |Z = z} = a + R

T b(t)z(t) dt under model (4). Let ϵπ = π0(Z, X)Y a R

T b(t)Z(t) dt be the residual variable. Let γ > 1 and β > γ/2 + 1 be two constants and the cj s be some generic positive constants.

Condition B

(B1) The variable ϵπ has a finite fourth moment, i.e., E(ϵ4 π) < .

Tan, Huang, Zhang and Yin

(B2) The functional variable Z has a finite fourth moment, i.e., R

T E{Z4(t)} dt < . The PC score ξj defined in (6) satisfies E(ξ4 j ) c1λ2 j, and the eigenvalue λj defined in (5) satisfies λj λj+1 c 1 2 j γ 1 for all j 1.

(B3) The coefficient bj satisfies |bj| c3j β, for all j 1.

Condition (B1) makes a mild restriction on the moment of ϵπ. Condition (B2) imposes regularity conditions on the random process Z, which requires that the differences of the adjacent eigenvalues do not decay too fast. Condition (B3) assumes an upper bound for the decay rate of bj, the coefficients of b projected on {φj} j=1. The latter two conditions are adopted from Hall and Horowitz (2007) for technical purposes.

Theorem 2 Under Conditions (A1) to (A7), (B1) to (B3) and model (4), assuming that k and h 0 as n and choosing the truncation parameter q n1/(γ+2β), we have Z

T {bb FSW(t) b(t)}2 dt = Op

k 2α + h2λk + ψZ{(log n)/n}k

ζ2(k)n(γ+1)/(γ+2β) .

The proof of Theorem 2 is given in Appendix E. The convergence rate of R

T {bb FSW(t) b(t)}2 dt in Theorem 2 can be quite slow mainly because of the small ball probability function µ(h). In the case of Gaussian process endowed with a metric, µ(h) has an exponential decay rate of h (Li and Shao, 2001), which leads to a non-infinitesimal rate of R

T {bb FSW(t) b(t)}2 dt. However, µ(h) can be much larger by choosing a proper semi-metric d (Ferraty and Vieu, 2006; Ling and Vieu, 2018), so that the final convergence rate of R

T {bb FSW(t) b(t)}2 dt reaches a polynomial decay of n. Another way to improve the convergence rate is to impose additional structure assumptions. For example, if we assume that the slope function b can be fully characterized by a finite number of PC basis functions, or the random process Z essentially lies on a finite dimensional (not necessarily Euclidean) space (Lin and Yao, 2021), then the truncation parameter q does not need to tend to infinity and the convergence rate of R

T {bb FSW(t) b(t)}2 dt can be much faster. Without such assumptions, estimating b in an infinite dimensional space with a nonparametric estimator of π0 leads to a cumbersome convergence rate. The OR estimator bb OR relies on functional principal component analysis and the partially linear model (10). Additional conditions are needed to derive the convergence rate of bb OR.

Condition C

(C1) The covariate X has finite fourth moment, i.e., E( X 4) < . For m = 1, . . . , p and j 1, | Cov{Xm, R

T Z(t)φj(t) dt}| < c4j γ β.

(C2) Let gm(t) = P j=1 Cov{Xm, R

T Z(t)φj(t) dt}φj(t)/λj and uim = Xim E(Xm) R

T gm(t){Zi(t) µZ(t)} dt. For m = 1, . . . , p, E(u2 1m|Z1) < , and E(u 1 u1) is positive definite with u1 = (u11, . . . , u1p) .

Conditions (C1) and (C2) are adopted from Shin (2009), which make mild assumptions on covariate X to ensure n-consistency of the estimated coefficient bθOR.

Functional Treatment Effect

Theorem 3 Under Conditions (B2), (B3), (C1), (C2) and model (10), choosing the truncation parameter q n1/(γ+2β), we have Z

T {bb OR(t) b(t)}2 dt = Op(n( 2β+1)/(γ+2β)) .

Because bb OR is theoretically equivalent to the least square estimator proposed by Shin (2009) provided that the estimators of each component in the additive model are unique (Buja et al., 1989), the theorem above is a direct result from Theorem 3.2 in Shin (2009) and thus its proof is omitted. The convergence rate of bb OR is the same as that of the estimated slope function in the traditional functional linear regression (Hall and Horowitz, 2007), i.e., such rate remains unchanged despite of the additional estimation of θ. Compared with Theorem 2, the convergence rate here Op(n( 2β+1)/(γ+2β)) = Op(n 1) Op(n(γ+1)/(γ+2β)) is much faster. Specifically, the part of the convergence rate Op ζ2(k) k 2α + h2λk + ψZ{(log n)/n}k/{nµ(h)} , which is the convergence rate of sup(z,x) Z X |bπhk(z, x) π0(z, x)|2, is replaced by a faster one, Op(n 1), because the nonparametric estimator bπhk is not used. According to the discussion at the end of Section 3.2, under model (4) and the assumptions of Theorem 1, bb DR is consistent. If assumptions of Theorem 2 hold, it has the same convergence rate as bb FSW. If the stronger model assumption (10) holds, then bb DR enjoys the same convergence rate as bb OR, provided that Condition C is satisfied. This differs from classical doubly robust estimators, whose consistency relies on correctly specifying one of the two parametric models.

5. Selection of Tuning Parameters

For the FSW estimator, our methodology requires a kernel function K with a metric d, basis functions νk and tuning parameters h, k and q. Using a compactly supported kernel, K would make the denominator of the first term in (9) too small for some Zi and destabilize the maximization of Hhk,Zi in (9). Therefore, we suggest using the Gaussian kernel, K(u) = exp( u2), with the semi-metric d being the L2 norm. To develop a general estimation strategy for π0, we standardize the Xi s and use a unified set of choices for νk. Letting Xmin = (mini{Xi1}n i=1, . . . , mini{Xip}n i=1) , Xmax = (maxi{Xi1}n i=1, . . . , maxi{Xip}n i=1) , we transform the Xi s as follows,

Xi,st = 2(Xi Xmin)

Xmax Xmin 1 ,

so that they have support on [ 1, 1]. Let Xij,st be the j-th component of Xi,st and Pℓbe the ℓ-th standard Legendre polynomial on [ 1, 1]. If k is given, we define v1(Xi) = 1 and vp(ℓ 1)+j+1(Xi) = Pℓ(Xij,st), for i = 1, . . . , n, j = 1, . . . , p and ℓ= 1, . . . , (k 1)/p. We restrict the range of k such that (k 1)/p is a positive integer, and orthonormalize the resulting matrix {νk(X1), . . . , νk(Xn)} . One may select h, k and q jointly by an L-fold cross-validation (CV), but this may be computationally intensive. We suggest a two-stage CV procedure as follows. In the first stage, we select the truncation parameter q using the OR estimator. Specifically, we randomly split the dataset into L parts, S1, . . . , SL. Let S ℓdenote the remaining sample

Tan, Huang, Zhang and Yin

with Sℓexcluded. We define the CV loss as

CVOR L (q) =

j=1 bbj, ℓ,OR

T bφj(t)bµZ(t) dt + bξij

where ba ℓ,OR, bbj, ℓ,OR and bθ ℓ,OR are obtained by applying the method in Section 3.2 to S ℓ. We choose the truncation parameter as bq CV = argminq {1,...,q99} CVOR L (q), where q99 is the number of principal components that can explain at least 99% of the variance of Z. Noting that the rate requirements for q in Theorem 1 and Theorem 2 are the same, we thus use the same bq CV in the FSW estimator. It remains to select h and k in the second stage. We consider

CVFSW L (h, k)

Yibπhk, ℓ(Zi, Xi) ba ℓ,FSW

j=1 bbj, ℓ,FSW

T bφj(t)bµZ(t) dt + bξij

where bξij = R

T {Zi(t) bµZ(t)}bφj(t) dt is the estimated PC score, ba ℓ,FSW, bbj, ℓ,FSW and bπhk, ℓare obtained by applying the method in Section 3.1 to S ℓ. We choose the one from a candidate set of {h, k} such that CVL(h, k) is minimized. As for the DR estimator, we still need to select the truncation parameter q. This can be achieved by minimizing the CV loss,

CVDR L (q) =

{Yi b E ℓ(Y |Xi, Zi)}bπhk, ℓ(Xi, Zi) + 1 |Sℓ|

j Sℓ { b E ℓ(Y |Xj, Zi)}

j=1 bbj, ℓ,DR

T bφj(t)bµZ(t) dt + bξij

where |{ }| denotes the number of elements in the set { }, ba ℓ,DR, b E ℓ(Y |Xi, Zi), bπhk, ℓ and bbj, ℓ,DR are obtained by applying the methods in Section 3.2 to S ℓ.

6. Simulation Study

In addition to the methods of functional stabilized weight (FSW), outcome regression (OR) and double robustness (DR) developed in Section 3, we also consider the method using the nonparametric principal component weight (PCW) proposed by Zhang et al. (2021), which assumes the same functional linear model for the ADRF; the weight-modified kernel ridge regression (WMKRR), a nonlinear method proposed by Wang et al. (2023); as well as the direct functional linear regression (FLR) and kernel ridge regression (KRR) of Y on Z. The FLR and KRR are expected to be biased if there are confounding effects. We select the tuning parameters h, k and q following the procedure in Section 5, and choose the number of PCs used in estimating the PCW so as to explain 95% of the variance of Z following Zhang et al. (2021). For a fair comparison, we use the truncation parameter bq CV in Section 5 for estimating b in all methods that utilize the functional linear model.

Functional Treatment Effect

Table 1: Mean (standard deviation) of the empirical mean squared error (MSE) of the ADRF with the best values highlighted in boldface.

n = 200 (i) (ii) (iii) (iv) (v) (vi) FSW 1.14 (0.66) 2.15 (1.36) 8.40 (8.60) 9.69 (10.08) 7.94 (12.86) 180.24 (42.74) OR 1.11 (0.68) 6.96 (3.42) 6.85 (9.77) 11.55 (14.11) 29.52 (11.79) 192.79 (58.15) DR 1.08 (0.66) 7.56 (3.65) 6.65 (9.56) 11.34 (14.36) 29.27 (12.12) 186.10 (52.82) PCW 2.68 (1.90) 12.98 (2.93) 13.10 (7.74) 16.13 (8.31) 16.17 (9.72) 186.84 (46.58) FLR 4.93 (1.63) 9.78 (2.29) 21.53 (7.90) 22.98 (7.72) 37.98 (14.08) 205.45 (61.24) WMKRR 8.30 (3.23) 22.67 (4.72) 25.07 (10.54) 25.57 (10.67) 34.71 (14.18) 133.79 (40.30) KRR 6.41 (1.83) 11.84 (2.42) 25.94 (6.59) 27.63 (6.60) 79.56 (21.87) 343.84 (158.34) n = 500 (i) (ii) (iii) (iv) (v) (vi) FSW 0.50 (0.34) 1.36 (0.95) 5.82 (4.27) 4.57 (2.98) 3.60 (14.88) 173.47 (28.29) OR 0.42 (0.33) 6.42 (2.58) 1.61 (2.52) 5.66 (2.53) 25.33 (6.13) 169.20 (41.27) DR 0.38 (0.31) 6.86 (2.74) 1.52 (2.45) 5.64 (2.41) 25.39 (6.55) 166.22 (36.65) PCW 1.68 (1.08) 12.16 (1.85) 8.42 (3.54) 11.50 (3.05) 12.09 (7.21) 169.64 (34.34) FLR 4.31 (1.05) 9.00 (1.51) 17.74 (3.95) 18.83 (3.60) 33.14 (8.37) 174.16 (46.45) WMKRR 6.62 (2.37) 19.62 (3.52) 23.89 (9.77) 24.80 (10.01) 35.26 (11.86) 125.12 (32.51) KRR 5.14 (1.06) 11.10 (1.54) 25.06 (4.22) 26.38 (4.33) 79.08 (14.62) 302.76 (79.42)

The code to reproduce the simulation results for our methods and the FLR is available at https://github.com/ruoxut/Functional Treatment. We consider the sample sizes n = 200 and 500, and generate data for models (i) to (v) as follows. For i = 1, . . . , n, we generate the functional treatment Zi(t), for t [0, 1], by Zi(t) = P6 j=1 Aijφj(t), where Ai1 = 4Ui1, Ai2 = 2

3Ui2, Ai3 = 2

2Ui3, Ai4 = 2Ui4, Ai5 = Ui5, Ai6 = Ui6/

2 with the Uij s being the independent standard normal variables, φ2m 1(t) =

2 sin(2mπt) and φ2m(t) =

2 cos(2mπt) for m = 1, 2, 3. For the slope function b, we define b(t) = 2φ1(t) + φ2(t) + φ3(t)/2 + φ4(t)/2. The Zi s and b are the same as those in Zhang et al. (2021). We consider five models for the covariate Xi and the outcome variable Yi:

(i) Xi = Ui1 + ϵi1 and Yi = 1 + R 1 0 b(t)Zi(t) dt + 2Xi + ϵi2,

(ii) Xi = Ui1 + 0.25ϵi1 and Yi = 1 + R 1 0 b(t)Zi(t) dt + 5 sin(Xi) + ϵi2,

(iii) Xi = (Xi1, Xi2) = {(Ui1+1)2+ϵi1, Ui2} and Yi = 1+ R 1 0 b(t)Zi(t) dt+2Xi1+2Xi2+ ϵi2,

(iv) Xi = (Xi1, Xi2) = {(Ui1 + 1)2 + ϵi1, Ui2} and Yi = 1 + R 1 0 b(t)Zi(t) dt + 2Xi1 + 2 cos(Xi1) + 5.5 sin(Xi2) + ϵi2,

(v) Xi = (Xi1, Xi2) = {(Ui1+1)2+ϵi1, Ui2} and Yi = 1+ R 1 0 b(t)Zi(t) dt+{ R 1 0 b(t)Zi(t) dt}2/25+ 2Xi1 + 5.5 sin(Xi2) + ϵi2,

where ϵi1 N(0, 1) and ϵi2 N(0, 25) are generated independently. In addition, we also consider a nonlinear model from Wang et al. (2023),

(vi) Xi = (Ui11, Ui21, Ui31, Ui41) , Yi = 10{Ui21U2 i11 +U2 i41 sin(2Ui31)}+0.5A2 i1 +4 sin(Ai1), Zi(t) = P4 j=1 Aij

2 sin(2jπt), for t [0, 1], where Ai1 = 4Ui11+Ui12, Ai2 = 2

Tan, Huang, Zhang and Yin

Ui22, Ai3 = 2

2Ui31 + Ui32, Ai4 = 2Ui41 + Ui42, with the Uijk s being the independent standard normal variables.

The PC scores Aij of Zi(t) affect covariate Xi linearly in models (i) and (ii) and nonlinearly in models (iii) and (iv), and covariate Xi affects the outcome variable Yi linearly in models (i) and (iii) and nonlinearly in models (ii) and (iv). The functional linear model (4) for the ADRF is correctly specified for models (i) to (iv), while the partially linear model (10) is only correctly specified for models (i) and (iii). Under models (v) and (vi), the functional linear model (4) for the ADRF is misspecified. Due to the confounding effect of X, the FLR and KRR ignoring X are expected to be biased. For each combination of model and sample size, we replicate 200 simulations and evaluate the results by the empirical mean squared error (MSE) of the ADRF following Wang et al. (2023), MSE = Pn i=1[ b E{Y (Zi)} E{Y (Zi)}]2/n, where b E( ) denotes a generic estimator of the ADRF.

Table 1 summarizes the mean and the standard deviation of the MSEs for all configurations. Our proposed estimators FSW, OR and DR outperform the other methods under models (i) (iv). In particular, our three estimators are close to each other and are better than the others under models (i) and (iii), where covariate X affects the outcome variable Y linearly. Under models (ii) and (iv) where the partially linear model (10) is misspecified, the FSW performs the best, followed by the OR and DR, and it is worth noting that the latter two still perform better than the other methods. Under model (v) where the functional linear model is misspecified, the FSW remarkably performs the best, followed by the PCW, while the FLR and WMKRR perform similarly. Under model (v), the FSW is expected to be biased. The performance of PCW is better than FLR but inferior to our methods under models (i) (iv). The WMRKK performs the best under model (vi), but its performance is not satisfactory under models (i) (v).

To give visualization of representative estimated slope functions, we show the quartile curves of the samples corresponding to the first, second and third quartile values of 200 MSEs for model (i) in Figure 2. It can be seen that the estimated curves using the FSW, OR and DR are close to the truth when n = 500, while the PCW is more variable and the FLR is more biased. Figure 3 corresponds to the results under model (iv), i.e., a more challenging model structure. The differences among the estimators are more significant: FSW performs the best, followed by the DR, OR and PCW, while the FLR is more biased than the others.

We also conduct a small-scale running time comparison of our estimator FSW and the WMKRR proposed by Wang et al. (2023). Both estimators contain the same weight function π0 estimated by different nonparametric methods. We compare the running time of computing the two estimators of π0, including selection of the tuning parameters. The experiments were carried out on a PC with an i7-12700 CPU and 16 GB RAM. As an illustration, model (i) is used as the data generation model with n = 100, 500 and 1000. We repeat the experiments 10 times. The values of averaged running time in seconds are 6.19 (n = 100), 37.28 (n = 500) and 236.88 (n = 1000) for the FSW, while the corresponding values for the WMKRR are 36.27, 204.77 and 440.67. We see that, for the given model, computing the estimator of π0 using the FSW is much faster than that using the WMKRR.

Functional Treatment Effect

Figure 2: True curve ( ), first ( ), second ( ) and third ( ) quartile estimated slope functions under model (i) with n = 200 and 500 for all methods.

7. Real Data Analysis

We illustrate the estimation of the ADRF using the five methods, FSW, OR, DR, PCW and FLR, on the electroencephalography (EEG) dataset from Ciarleglio et al. (2022). The EEG is a relatively low-cost tool to measure human neuronal activities. The measurements are taken from multiple electrodes placed on the scalps of subjects, and they are then processed and transformed to produce current source density curves on the frequency domain, which can provide information on the intensity of neuronal activities. In particular, the

Tan, Huang, Zhang and Yin

Figure 3: True curve ( ), first ( ), second ( ) and third ( ) quartile estimated slope functions under model (iv) with n = 200 and 500 for all methods.

frontal asymmetry curves, considered potential biomarkers for major depressive disorder, are treated as the functional treatment variable Z(t), for t [4, 31.75] Hz. The outcome variable Y is defined as log(e Y + 1), where e Y is the quick inventory of depressive symptomatology (QIDS) score measuring the severity of depressive symptomatology. A larger value of Y indicates a more severe depressive disorder. We investigate the causal effect of neuronal activities represented by the frontal asymmetry curves on the severity of major depressive disorder. Potential confounding covariates X = (X1, X2, X3) include age X1, sex X2 (1 for female and 0 for male) and Edinburgh Handedness Inventory score X3 (rang-

Functional Treatment Effect

5 10 15 20 25 30 -0.08

5 10 15 20 25 30 -0.6

FSW OR DR PCW FLR 2.2

Figure 4: Left: the estimated slope functions using FSW ( ), OR ( ), DR ( ), PCW ( ) and FLR ( ); middle: the frontal asymmetry curves corresponding to the smallest ( ), median (+) and largest ( ) average values over 4 to 15 Hz; right: the estimated E{Y (Z)} of the three selected curves with marks corresponding to the middle panel.

ing from 100 to 100 corresponding to completely left to right handedness). Individuals with missing Z (29.6% of the total sample) are removed, which results in a sample of size 85 males and 151 females. The means (standard deviations) of Y , X1 and X3 are 2.73 (0.72), 35.97 (13.07) and 72.69 (48.76), respectively. For comparison of the estimated ADRF using different methods, the confounding variables X are centralized.

We use the number of PCs explaining 95% of the variance of Z (equal to 11 for this dataset) in estimating the PCW. We choose the truncation parameter q by bq CV to estimate the slope function b for all methods and choose h and k used in FSW by minimizing CVFSW L as in Section 5. As a result, bq CV = 2 is used, which can explain about 80% of the variance of Z.

The left panel of Figure 4 shows the estimated slope functions. All the estimated slope functions have the same increasing trend over frequencies, while they have a negative effect on the outcome variable for low frequencies. In other words, subjects with higher frontal asymmetry values on the low-frequency domain tend to be healthier. This is consistent with the observations in Ciarleglio et al. (2022), where it was found a negative association between the frontal asymmetry curves and the major depressive disorder status using their functional regression model. More specifically, in the left panel of Figure 4, the estimated slope function of FSW is the steepest and shows the largest effect in the low-frequency domain. In contrast, the estimated slope function of FLR has the overall smallest absolute values; the other three curves are quite close to each other.

To visualize the above negative effect in the estimated slope function, we compute the estimated ADRF E{Y (z)} of three curves, which correspond to the smallest, median and largest average frontal asymmetry values over low frequencies 4 to 15 Hz. The middle panel of Figure 4 exhibits the selected curves and the right panel shows the corresponding estimated E{Y (z)} using all methods, which are around 2.9, 2.8 and 2.3 for the three curves, respectively.

Tan, Huang, Zhang and Yin

8. Discussion

In the case of a functional treatment variable, we propose three estimators of the ADRF, namely, the FSW, OR and DR estimators, based on the functional linear model for the ADRF. The consistency of the FSW estimator relies on the functional linear model of E{Y (z)} by developing a nonparametric estimator of the weight; while the OR estimator requires a more restrictive linear model for Y regressed on (Z, X). The DR estimator is consistent if either of the first two estimators is consistent. It is of interest to construct the confidence band for the slope function to better quantify the ADRF or ATE, which, however, is difficult even in the simpler context of functional linear regression. As shown in Cardot et al. (2007), it is impossible for an estimator of the slope function to converge in distribution to a nondegenerate random element under the norm topology. For the OR estimator, it is possible to adapt the method proposed by Imaizumi and Kato (2019) to our context to construct an approximate confidence band. For the other two estimators, the construction of the confidence band is challenging due to the less restrictive modelling assumption and the nonparametric estimator of the FSW, which warrants further investigation.

Acknowledgements

We thank the Action Editor and two referees for their careful reviews that greatly improved our work. We also thank Dr. Jiayi Wang for sharing her code. The research of Zheng Zhang is supported by the funds from the National Key R&D Program of China [grant number 2022YFA1008300] and the National Natural Science Foundation of China [grant number 12371284]. Tan s research was supported by NSFC (No. 12401363 and 12471263), the Fundamental Research Funds for the Central Universities, and Key Laboratory of Intelligent Computing and Applications (Ministry of Education). Yin s research was supported by the Research Grants Council of Hong Kong (17308321), the Theme-based Research Scheme (TRS) Project T45-401/22-N, and the Patrick SC Poon Endowment Fund.

Functional Treatment Effect

Appendix A. Proof of Proposition 1

Proof The if part follows from (7), and thus we show the only if part below. Suppose

E{π(z, X)v(X)|Z = z} = E{v(X)} ,

for all integrable functions v. Since E{v(X)} = E{π0(z, X)v(X)|Z = z} for all integrable functions v, we have

E[{π(z, X) π0(z, X)}v(X)|Z = z}] = 0 .

Taking v(X) = exp(ia X) for all a Rp, we have

E[{π(z, X) π0(z, X)} exp(ia X)|Z = z}] = 0 ,

According to the uniqueness of the Fourier transformation, we conclude that π(z, ) = π0(z, ) a.s.

Appendix B. Generalized empirical likelihood interpretation of bπhk

To gain more insight on the function ρ, we investigate the generalized empirical likelihood interpretation of our estimator. We show that the estimator defined by (9) is the dual solution to the following local generalized empirical likelihood maximization problem. For each Zi, i = 1, . . . , n,

max {π(Zi,Xj)}j S i

P j S i D{π(Zi, Xj)}K{d(Zj, Zi)/h}

P j S i K{d(Zj, Zi)/h} ,

P j S i π(Zi, Xj)νk(Xj)K{d(Zj, Zi)/h}

P j S i K{d(Zj, Zi)/h} = 1 n 1

j S i νk(Xj) , (12)

where D(w) is a distance measure between w and 1 for w R. The function D(v) is continuously differentiable satisfying D(1) = 0 and ρ( u) = D{D ( 1)(u)} u D ( 1)(u), where D ( 1) is the inverse function of D , the first derivative of D. Different choices of ρ correspond to different distance measures. For example, if we take ρ(v) = exp( v 1), then D(v) = v log(v) is the information entropy and the weights correspond to the exponential tilting (Imbens et al., 1998). Other choices include ρ(v) = log(v) + 1 with D(v) = log(v), the empirical likelihood (Owen, 1988), and ρ(v) = (1 v)2/2 with D(v) = (1 v)2/2, yielding the implied weights of the continuous updating of the generalized method of moments (Hansen et al., 1996). Following Tseng and Bertsekas (1991), we show that the dual problem of (12) is to maximize Hhk,Zi in (9). For each i = 1, . . . , n, we let νk = P j S i νk(Xj)/(n 1), Kij = K{d(Zj, Zi)/h}/[P j S i K{d(Zj, Zi)/h}], πij = π(Zi, Xj) and wij = Kijπij, for j S i. Moreover, let Dij(ν) = D(ν/Kij), Vk = νk(Xj)

j S i Rk (n 1), wi = (wij) j S i, and F(wi) = P j S i Dij(wij)Kij. Then, (12) can be written as

min wi F(wi) , s.t. Vkwi = νk . (13)

Tan, Huang, Zhang and Yin

We define the conjugate convex function of F as

F (u) = sup wi

j S i {ujwij Kij Dij(wij)}

= sup {πij}j S i

j S i Kij{ujπij D(πij)}

j S i Kij{ujπ ij D(π ij)} ,

where the π ij s satisfy

uj = D (π ij) π ij = D ( 1)(uj), j = 1, . . . , n ,

by taking the first order condition. It follows that

j S i Kijρ( uj) , (14)

where ρ( u) = D{D ( 1)(u)} u D ( 1)(u). Following Tseng and Bertsekas (1991) and using (14), we conclude that the dual problem of (12) is

max η Rk{ F (η Vk) + η νk} = max η Rk X

Kijρ{ η νk(Xj)} + η νk(Xj)

= max η Rk X

Kijρ{η νk(Xj)} η νk(Xj)

= max η Rk Hhk,Zi(η) .

Appendix C. Proof of Double Robustness in (11)

Recall from Section 2 that E{Y (z)} = E{Y π0(X, Z)|Z = z} = E{E(Y |X, Z = z)}. First assuming that e E(Y |X, Z) = E(Y |X, Z), then we have

E[{Y e E(Y |X, Z)}eπ(X, Z) + EX{ e E(Y |X, Z)}|Z = z]

= E[{Y E(Y |X, Z)}eπ(X, Z)|Z = z] + E{E(Y |X, Z = z)}

= E{Y (z)} ,

where the second equality follows from the law of total expectation. Assuming that eπ(X, Z) = π0(X, Z) = f X(X)/f X|Z(X|Z), we have

E[{Y e E(Y |X, Z)}eπ(X, Z) + EX{ e E(Y |X, Z)}|Z = z]

= E{Y π0(X, Z)|Z = z} E{ e E(Y |X, Z)π0(X, Z)|Z = z} + EX{ e E(Y |X, Z)}|Z = z]

= E{Y (z)} Z e E(Y |X = x, Z = z)f X(x) dx + Z e E(Y |X = x, Z = z)f X(x) dx

= E{Y (z)} .

Therefore, the doubly robust property is proved.

Functional Treatment Effect

Appendix D. Proof of Theorem 1

Proof Recall from Section 3.1 that

bπhk(Zi, Xi) = ρ {bη Ziνk(Xi)} ,

where bηZi is a k-dimensional vector maximizing the strictly concave function Hhk,Zi defined by

Hhk,Zi(η) =

P j S i ρ{η νk(Xj)}K{d(Zj, Zi)/h}

P j S i K{d(Zj, Zi)/h} η 1 n 1

j S i νk(Xj) .

For a fixed z L2(T ), we have

bπhk(z, Xi) = ρ {bη z νk(Xi)} ,

where bηz is a k-dimensional vector maximizing Hhk,z defined by

P j S i ρ{η νk(Xj)}K{d(Zj, z)/h}

P j S i K{d(Zj, z)/h} η 1 n 1

j S i νk(Xj) . (15)

Here, S i is independent of Zi for i = 1, . . . , n. Let H hk,z(η) be the theoretical counterpart of Hhk,z(η),

H hk,z(η) = E[ρ{η νk(X)}K{d(Z, z)/h}]

E{K{d(Z, z)/h}} η E{νk(X)} ,

η z = argmaxη Rp H hk,z(η) and π hk(z, x) = ρ {η z νk(x)}. Note that

sup (z,x) Z X |bπhk(z, x) π0(z, x)|

sup (z,x) Z X |bπhk(z, x) π hk(z, x)| + sup (z,x) Z X |π hk(z, x) π0(z, x)|

To prove Theorem 1, it suffices to prove the convergence rates for the two terms on the right-hand side of the above inequality. Lemmas 2 and 3 provide such results. Lemma 1 is needed to prove Lemmas 2 and 3.

Lemma 1 Let

bek1,k2(z) =

P j S i gk1,k2,z(Xj)K{d(Zj, z)/h}

(n 1)E[K{d(Zj, z)/h}] ,

where gk1,k2,z is defined in (A6), and

P j S i K{d(Zj, z)/h}

(n 1)E[K{d(Zj, z)/h}] .

Tan, Huang, Zhang and Yin

Under Conditions (A4) to (A7), we have

sup z Z | bf(z) 1| = Oa.c.

ψZ{(log n)/n}

n=1 P inf z Z bf(z) < 1

sup z Z |E{bek1,k2(z)} E[gk1,k2,z(X)|Z = z]| = O(hλ) , k1, k2 {1, . . . , k},

sup z Z |bek1,k2(z) E{bek1,k2(z)}| = Oa.c.

ψZ{(log n)/n}

, k1, k2 {1, . . . , k},

hold uniformly over k.

Proof Note that for all k1, k2 = 1, . . . , k, gk1,k2,z is a scalar function, and Condition (A6) is assumed uniformly over k. Therefore, Lemmas 8 to 11 in Ferraty et al. (2010) can be used to derive the lemma above. Also, note the Oa.c. denotes the almost complete convergence, which implies the almost sure convergence and convergence in probability.

Lemma 2 Under Conditions (A1) to (A6), we have

sup (z,x) Z X |π hk(z, x) π0(z, x)| = O{ζ(k)(k α + hλ

Proof We first prove that

sup (z,x) Z X |π0(z, x) ρ {η z νk(x)}| = O(k α) , (16)

where ηz is given in Condition (A2). Recall that ρ is a strictly increasing and concave function, so ρ ( 1) is strictly decreasing. Thus, we have

ρ ( 1)(c2) inf (z,x) Z X ρ ( 1){π0(z, x)} sup (z,x) Z X ρ ( 1){π0(z, x)} ρ ( 1)(c1) ,

for some positive constants c1, c2, because π0(z, x) is bounded away from zero and infinity by Condition (A1). According to Condition (A2), there exists a constant C1 > 0 such that

sup (z,x) Z X |ρ ( 1){π0(z, x)} η z νk(x)| C1k α .

Thus, we have, z Z and x X,

η z νk(x) [ρ ( 1){π0(z, x)} C1k α, ρ ( 1){π0(z, x)} + C1k α]

[ρ ( 1)(c2) C1k α, ρ ( 1)(c1) + C1k α] (17)

Functional Treatment Effect

ρ {η z νk(x) + C1k α} ρ {η z νk(x)}

< π0(z, x) ρ {η z νk(x)}

< ρ {η z νk(x) C1k α} ρ {η z νk(x)} . (18)

By the mean value theorem, for k large enough, there exist

ξ1(z, x) (η z νk(x), η z νk(x) + C1k α) [ρ ( 1)(c2) C1k α, ρ ( 1)(c1) + 2C1k α] Γ1 ,

ξ2(z, x) (η z νk(x) C1k α, η z νk(x)) [ρ ( 1)(c2) 2C1k α, ρ ( 1)(c1) + C1k α] Γ1 ,

ρ {η z νk(x) + C1k α} ρ {η z νk(x)} = ρ {ξ1(z, x)}C1k α inf u Γ1 ρ (u)C1k α ,

ρ {η z νk(x) C1k α} ρ {η z νk(x)} = ρ {ξ2(z, x)}C1k α sup u Γ1 { ρ (u)}C1k α , (19)

Γ1 = [ρ ( 1)(c2) 1, ρ ( 1)(c1) + 1].

Note that max{ infu Γ1 ρ (u), supu Γ1{ ρ (u)}} is a positive and finite constant, because Γ1 is compact and ρ is continuous and strictly negative. Thus, the claim (16) is proved using (18) and (19). Next, we compute a bound for supz Z H hk,z(ηz) . We have

sup z Z H hk,z(ηz)

E[ρ {η z νk(X)}νk(X)K{d(Z, z)/h}]

E[K{d(Z, z)/h}] E{νk(X)}

E[ρ {η z νk(X)}νk(X)K{d(Z, z)/h}]

E[K{d(Z, z)/h}] E[ρ {η z νk(X)}νk(X)|Z = z]

E [ρ {η z νk(X)} π0(z, X)]νk(X)|Z = z

k =1 sup z Z

E[ρ {η z νk(X)}vk (X)K{d(Z, z)/h}]

E[K{d(Z, z)/h}] E[ρ {η z νk(X)}vk (X)|Z = z]

+ sup (z,x) Z X |ρ {η z νk(x)} π0(z, x)| = O(

khλ + k α) , (20)

where Lemma 1 and equation (16) are used for the last inequality. For all z Z and for some constant C2 > 0 (to be chosen later), define the set

Λz = {η Rk : η ηz C2(k α + hλ

By Condition (A3), we have η Λz,

sup x X |η νk(x) η z νk(x)| η ηz sup x X νk(x) < C2(k α + hλ

Tan, Huang, Zhang and Yin

so that x X and for sufficiently large k, by (17),

[ρ ( 1)(c2) C1k α C2(k α + hλ

k)ζ(k), ρ ( 1)(c1) + C1k α + C2(k α + hλ

k)ζ(k)] Γ1 . (21)

Based on (20), there exists a C4 > 0 that is independent of z such that supz Z H hk,z(ηz) < C4(k α + hλ

k). For any η Λz, i.e. η ηz = C2(k α + hλ

k), by the mean value theorem we have

H hk,z(η) H hk,z(ηz)

= (η ηz) H hk,z(ηz) + 1

2(η ηz) 2H hk,z( ηz)(η ηz)

η ηz H hk,z(ηz)

2(η ηz) E[ρ { η z νk(X)}K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}] (η ηz)

η ηz H hk,z(ηz) a1

2 (η ηz) [E{νk(X)νk(X) |Z = z} + O(hλ)Jk k](η ηz)

η ηz n H hk,z(ηz) a1{λmin + O(hλk)}

η ηz {C4(k α + hλ

k) a1{λmin + o(1)}C2(k α + hλ

where Jk k is the k k matrix of ones, ηz lies between η and ηz on Λz, a1 = infu Γ1{ ρ (u)} > 0 uniformly in z Z, λmin = infz Z λz,min with λz,min the smallest eigenvalue of E{νk(X)νk(X) |Z = z}, and we used O(hλk) = o(1), Lemma 1 and Condition (A3). Therefore, we choose

C2 > 2C4 a1λmin ,

so that H hk,z(η) < H hk,z(ηz), for η Λz. Since H hk,z is continuous, there is a local maximum of H hk,z in the interior of Λz. On the other hand, H hk,z is a strictly concave function with the unique maximum at η z. Therefore, we conclude that η z Λo z, i.e. η z ηz < C2(k α + hλ

k). Note that C2 is independent of z so that

sup z Z η z ηz < C2(k α + hλ

Recalling (21) and using (22) and Condition (A3), for large enough k, we have

sup (z,x) Z X |ρ {η z νk(x)} ρ {η z νk(x)}|

= sup (z,x) Z X |ρ {ξ (z, x)}||η z νk(x) η z νk(x)|

a2 sup z Z η z ηz sup x X νk(x) a2C2(k α + hλ

k)ζ(k) , (23)

Functional Treatment Effect

where ξ (z, x) Γ1 lies between η z νk(x) and η z νk(x) and a2 = supu Γ1 |ρ (u)|. Finally, using (16) and (23), we conclude that

sup (z,x) Z X |π hk(z, x) π0(z, x)|

sup (z,x) Z X |π hk(z, x) ρ {η z νk(x)}| + sup (z,x) Z X |ρ {η z νk(x)} π0(z, x)|

= O{ζ(k)(k α + hλ

Lemma 3 Under Conditions (A1) to (A7), we have

sup (z,x) Z X |bπhk(z, x) π hk(z, x)| = O ζ(k)

ψZ{(log n)/n}k

Proof The key step is to provide the convergence rate of supz Z bηz η z , for which we need to provide the convergence rate of supz Z Hhk,z(η z) . For all z Z, we write

P j S i K{d(Zj, z)/h}νk(Xj)νk(Xj)

P j S i K{d(Zj, z)/h}

= E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}] + (n 1)E[K{d(Z, z)/h}] P j S i K{d(Zj, z)/h}

P j S i K{d(Zj, z)/h}νk(Xj)νk(Xj)

(n 1)E[K{d(Z, z)/h}] E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}]

1 (n 1)E[K{d(Z, z)/h}] P j S i K{d(Zj, z)/h}

E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}]

E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}] + (n 1)E[K{d(Z, z)/h}] P j S i K{d(Zj, z)/h} A1,z,k

A2,z E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}] . (24)

It follows from Lemma 1 that

sup z Z A1,z,k = O s

ψZ{(log n)/n}

Jk k and sup z Z A2,z = O s

ψZ{(log n)/n}

Using Lemma 1 again, we have

(n 1)E[K{d(Z, z)/h}] P i S i K{d(Zi, z)/h} 1 ,

Tan, Huang, Zhang and Yin

and for any fixed k and all z Z

E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}] E[νk(X)νk(X) |Z = z] ,

as n , whose eigenvalues are bounded away from zero and infinity uniformly in k and z by Condition (A3).

Combing the results above, we deduce that for all z Z,

b Sn = E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}] + O s

ψZ{(log n)/n}

Jk k . (25)

By Conditions (A3) and (A5), for h sufficiently small, there exist two positive constants s1 and s2 that are independent of z such that

0 < s1 λmin

E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}]

E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}]

s2 < , (26)

where λmin(A) and λmax(A) denote the minimum and maximum eigenvalues of A, respectively. Consider the event set for all z Z,

En(z) = (η η z) b Sn(η η z)

> (η η z) E[K{d(Z, z)/h}νk(X)νk(X) ]

E[K{d(Z, z)/h}] s1

(η η z), η = η z

Using (25) and (26), we deduce that for h sufficiently small,

P En(z) = 1 . (27)

On the other hand, noting that by defition of η z,

E[ρ {η z νk(X)}K{d(Z, z)/h}νk(X)]

E[K{d(Z, z)/h}] = E{νk(X)} ,

Functional Treatment Effect

sup z Z Hhk,z(η z)

P j S i ρ {η z νk(Xj)}K{d(Zj, z)/h}νk(Xj) P j S i K{d(Zj, z)/h} 1 n 1

j S i νk(Xj)

ρ {η z νk(Xj)}K{d(Zj, z)/h}νk(Xj)

E[K{d(Z, z)/h}]

E[ρ {η z νk(X)}K{d(Z, z)/h}νk(X)]

E[K{d(Z, z)/h}]

+ E{νk(X)} 1 n 1

j S i νk(Xj)

P j S i ρ {η z νk(Xj)}K{d(Zj, z)/h}νk(Xj) P j S i K{d(Zj, z)/h}

P j S i K{d(Zj, z)/h}

(n 1)E[K{d(Z, z)/h}] 1

ψZ{(log n)/n}k

ψZ{(log n)/n}k

ψZ{(log n)/n}k

where Lemma 1 is used for the second equality. This means that there exists a constant C3 > 0 that is independent of z such that

sup z Z Hhk,z(η z) < C3

ψZ{(log n)/n}k

nµ(h) . (28)

For some C4 > 0 (to be chosen later), define the set

bΛz = η Rk : η η z C4C3

ψZ{(log n)/n}k

For all η bΛz, x X, and for sufficiently large n and k, by Condition (A3), (17) and (22), we have

η νk(x) η z νk(x) η η z νk(x) C4C3

ψZ{(log n)/n}k

η z νk(x) C4C3

ψZ{(log n)/n}k

nµ(h) ζ(k), η z νk(x) + C4C3

ψZ{(log n)/n}k

ρ ( 1)(c2) C1k α C2(k α + hλ

ψZ{(log n)/n}k

nµ(h) ζ(k),

ρ ( 1)(c1) + C1k α + C2(k α + hλ

ψZ{(log n)/n}k

Γ1 [ρ ( 1)(c2) 1, ρ ( 1)(c1) + 1] . (29)

Tan, Huang, Zhang and Yin

By the mean value theorem, for any η bΛz, we have, based on the fact that P(En(z)) = 1 from (27),

Hhk,z(η) Hhk,z(η z)

= (η η z) Hhk,z(η z) + 1

2(η η z) 2Hhk,z( η)(η η z)

η η z Hhk,z(η z) a

2(η η z) b Sn(η η z)

< η η z Hhk,z(η z) a

2(η η z) E[K{d(Z, z)/h}νk(X)νk(X) }

E[K{d(Z, z)/h}] s1

η η z Hhk,z(η z) a

2(η η z) s1Ik k s1

η η z Hhk,z(η z) a

4s1 η η z , (30)

where η bΛz lies between η and η z, a = infu Γ1{ ρ (u)} > 0 and the third inequality follows from (26). Therefore, we choose

C4 > 4 as1 ,

so that Hhk,z(η) < Hhk,z(ηz), for η bΛz. On the other hand, we know that bηz is the unique maximum point of Hhk,z(η). It follows that bηz bΛo z(ϵ), i.e. bηz η z C4C3 p

ψZ{(log n)/n}k/nµ(h). Noting that C4 and C3 are independent of z, we conclude that

sup z Z bηz η z = O s

ψZ{(log n)/n}k

Finally, we are able to prove the convergence rate for sup(z,x) Z X |bπhk(z, x) π hk(z, x)|. By the mean value theorem, we have

bπhk(z, x) π hk(z, x) = ρ {bη z νk(x)} ρ {η z νk(x)} = ρ { η z νk(x)}(bηz η z) νk(x) ,

where η z lies between bηz and η z. From (29) and (31), we have

sup (z,x) Z X |ρ { η z νk(x)}| = O(1) . (32)

It follows that by combining (31), (32) and Condition (A3),

sup (z,x) Z X |bπhk(z, x) π hk(z, x)|

sup (z,x) Z X |ρ { η z νk(x)}| sup z Z bηz η z sup x X νk(x) = O ζ(k)

ψZ{(log n)/n}k

Functional Treatment Effect

Appendix E. Proof of Theorem 2

To aid the proof of Theorem 2, we need the following lemma.

Lemma 4 Let Ed = Ed(n) = {λd/2 b G G O}, where M O = R R M2(s, t) ds dt for a symmetric bivariate function M. Under Conditions (B1) to (B3), we have

Proof This follows from Section 5.1 of Hall and Horowitz (2007).

We give the formal proof of Theorem 2. Proof From (5.2) of Hall and Horowitz (2007), we know that supj |bλj λj| b G G O, where O is defined in Lemma 4. Let Eq = {λq/2 b G G O}. We conclude that

1 2λj bλj 3

for j = 1, . . . , q, provided that the event Eq holds. We see from Lemma 4 that P(Eq) 1 as n . Therefore, since our result is probabilistic, we can argue under the assumption that the event Eq holds. Recall that bb FSW = Pq j=1bb FSW,j bφj, where bb FSW,j = bλ 1 j be FSW,j. Write bb FSW,j = ˇbj + bλ 1 j (Sj1 + Sj2 + Sj3), where ˇbj = bλ 1 j R

T e(t)φj(t) dt, Sj1 = R

T {be FSW(t) e(t)}φj(t) dt,

T e(t){bφj(t) φj(t)} dt and Sj3 = R

T {be FSW(t) e(t)}{bφj(t) φj(t)} dt. Using (33) and the fact that S2 j3 be FSW e 2 bφj φj 2, we have

j=1 (bb FSW,j ˇbj)2 3

j=1 bλ 2 j (S2 j1 + S2 j2 + S2 j3)

j=1 λ 2 j (S2 j1 + S2 j2 + S2 j3)

j=1 λ 2 j (S2 j1 + S2 j2) + 12

j=1 λ 2 j be FSW e 2 bφj φj 2 , (34)

Also, we have

j=1 bj bφj(t) b(t) 2 dt 2 Z

j=1 bj{bφj(t) φj(t)} 2 dt + 2

j=1 b2 j bφj φj 2 + 2

j=q+1 b2 j . (35)

Tan, Huang, Zhang and Yin

Using (34) and (35), we have

T {bb FSW(t) b(t)}2 dt

j=1 bb FSW,j bφj(t)

j=1 ˇbj bφj(t) +

j=1 ˇbj bφj(t)

j=1 bj bφj(t) +

j=1 bj bφj(t) b(t) 2 dt

j=1 (bb FSW,j ˇbj)2 + 3

j=1 (ˇbj bj)2 + 3 Z

j=1 bj bφj(t) b(t) 2 dt

j=1 λ 2 j (S2 j1 + S2 j2) + 36

j=1 (λ 2 j be FSW e 2 + qb2 j) bφj φj 2 + 6

+ op(n (2β 1)/(γ+2β)) , (36)

where Pq j=1(ˇbj bj)2 = op(n (2β 1)/(γ+2β)) is used for the last inequality. This can be proved following the same argument of proving (5.11) in Hall and Horowitz (2007) with the response variable Y replaced by Y π0(Z, X). Also, using q n1/(γ+2β) and Condition (B3), we have P j=q+1 b2 j R q t 2β dt = O(n (2β 1)/(γ+2β)). From Condition (B2), we know that there exists a constant C1 > 1 such that λj C1j γ. From (5.21) and (5.22) of Hall and Horowitz (2007), we know that bφj φj 2 = Op(j2n 1). Finally, using (5.30) and (5.31) of Hall and Horowitz (2007), we have Pq j=1 λ 2 j S2 j2 = Op(n (2β 1)/(γ+2β)). Combining these results with (36), we deduce that there exists a constant C2 > 0 such that

T {bb FSW(t) b(t)}2 dt C2

j=1 λ 2 j (S2 j1 + be FSW e 2j2n 1) + Op(n (2β 1)/(γ+2β))

It remains to provide the convergence rates of S2 j1 and be FSW e 2. Recall that

be FSW(t) = 1

bπhk(Zi, Xi)Yi bµY,π {Zi(t) bµ(t)} ,

where bµY,π = Pn i=1 bπhk(Zi, Xi)Yi/n. Let

π0(Zi, Xi)Yi eµY,π {Zi(t) bµ(t)} ,

where eµY,π = Pn i=1 π0(Zi, Xi)Yi/n. Also, recall from the Karhunen-Lo eve expansion that Zi(t) = µ(t) + P j=1 ξijφj(t), where the PC scores ξij satisfy E(ξij) = 0, E(ξ2 ij) = λj and E(ξijξij ) = 0 for j = j .

Functional Treatment Effect

Defining ϵ = Pn i=1 ϵi/n and ξj = Pn i=1 ξij/n, we can write

T {be FSW(t) ee(t) + ee(t) e(t)}φj(t) dt

i=1 {Yibπhk(Zi, Xi) Yiπ0(Zi, Xi) bµY,π + eµY,π} Z

T {Zi(t) bµ(t)}φj(t) dt

i=1 {Yiπ0(Zi, Xi) eµY,π} Z

T {Zi(t) bµ(t)}φj(t) dt E[{Y π0(Z, X) µY,π}ξj]

i=1 Yi{bπhk(Zi, Xi) π0(Zi, Xi)}(ξij ξj)

T b(t){Zi(t) bµ(t)} dt + ϵi ϵ (ξij ξj) E[{Y π0(Z, X) µY,π}ξj]

= Aj1 + Aj2 + Aj3 + Aj4 , (38)

i=1 Yi{bπhk(Zi, Xi) π0(Zi, Xi)}(ξij ξj) ,

T b(t)Zi(t) dt E ξij

T b(t)Zi(t) dt ,

i=1 {ξijϵi E(ξijϵi)} ,

T b(t)bµ(t) dt ξj ϵ .

We next provide the convergence rates for the Ajk s sequentially. First note that E(ξij ξj)2 = (n 1)E(ξ2 j )/n = O(λj) by Condition (B2). It follows that

i=1 (ξij ξj)2 = Op(λj) . (39)

Using the Cauchy-Schwarz inequality, we have

|Aj1| sup (z,x) Z X |bπhk(z, x) π0(z, x)| 1

i=1 (ξij ξj)2 1/2

= O ζ(k) k α + hλ

ψZ{(log n)/n}k

Op(1) Op(λ1/2 j )

λ1/2 j ζ(k) k α + hλ

ψZ{(log n)/n}k

Tan, Huang, Zhang and Yin

based on Theorem 1 and (39). For Aj2, we have

T b2(t)Z2 i (t) dt 1

n{E(ξ4 j )}1/2 E Z

T b4(t)Z4(t) dt 1/2 = O λj

where the last equality follows from Condition (B2). Similarly, we have

n E(ξ2 j ϵ2) 1

E(ξ4 j )E(ϵ4) = O λj

where the last equality follows from Conditions (B2) and (C1). It follows that

Note that bµ µ and ϵ E(ϵ) a.s. as n by the law of large numbers, and µ and E(ϵ) are bounded. It follows that R

T b(t)bµ(t) dt and ϵ are bounded for large n. Thus, there exists a constant C3 > 0 such that for n sufficiently large,

|Aj4| C3 ξj = Op nq

E( ξ2 j ) o = Op

Combing (38), (39), (41) and (42), we conclude that

λjζ2(k) k 2α + h2λk + ψZ{(log n)/n}k

As for be FSW e 2, we have

be FSW e 2 2 be FSW ee 2 + 2 ee e 2 ,

|be FSW(t) ee(t)|

i=1 Yi{bπhk(Zi, Xi) π0(Zi, Xi)}{Zi(t) bµ(t)}

sup (z,x) Z X |bπhk(z, x) π0(z, x)| 1

i=1 {Zi(t) bµ(t)}2 1/2 .

It follows by using Theorem 1 and the fact that Pn i=1 R

T {Zi(t) bµ(t)}2 dt/n R

T var{Z(t)} dt = O(1) as n that

be FSW ee 2 = Op

ζ2(k) k 2α + h2λk + ψZ{(log n)/n}k

Also, the standard result on the convergence rate of the empirical mean estimator gives ee e 2 = Op(n 1). Therefore, we conclude that

be FSW e 2 = Op

ζ2(k) k 2α + h2λk + ψZ{(log n)/n}k

The proof is completed by combining (37), (43) and (44).

Functional Treatment Effect

Chunrong Ai, Oliver Linton, Kaiji Motegi, and Zheng Zhang. A unified framework for efficient estimation of general treatment models. Quantitative Economics, 12(3):779 816, 2021.

Susan Athey, Guido W Imbens, and Stefan Wager. Approximate residual balancing: debiased inference of average treatment effects in high dimensions. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 80(4):597 623, 2018.

Matteo Bonvini and Edward H Kennedy. Fast convergence rates for dose-response estimation. ar Xiv preprint ar Xiv:2207.11825, 2022.

Andreas Buja, Trevor Hastie, and Robert Tibshirani. Linear smoothers and additive models. The Annals of Statistics, pages 453 510, 1989.

Herv e Cardot, Andr e Mas, and Pascal Sarda. Clt in functional linear regression models. Probability Theory and Related Fields, 138(3):325 361, 2007.

Kwun Chuen Gary Chan, Sheung Chi Phillip Yam, and Zheng Zhang. Globally efficient non-parametric inference of average treatment effects by empirical balancing calibration weighting. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 78(3):673 700, 2016.

Jingxiang Chen, Haoda Fu, Xuanyao He, Michael R Kosorok, and Yufeng Liu. Estimating individualized treatment rules for ordinal treatments. Biometrics, 74(3):924 933, 2018.

Xiaohong Chen. Large sample sieve estimation of semi-nonparametric models. Handbook of Econometrics, 6:5549 5632, 2007.

Adam Ciarleglio, Eva Petkova, R Todd Ogden, and Thaddeus Tarpey. Treatment decisions based on scalar and functional baseline covariates. Biometrics, 71(4):884 894, 2015.

Adam Ciarleglio, Eva Petkova, Thaddeus Tarpey, and R Todd Ogden. Flexible functional regression methods for estimating individualized treatment rules. Stat, 5(1):185 199, 2016.

Adam Ciarleglio, Eva Petkova, Todd Ogden, and Thaddeus Tarpey. Constructing treatment decision rules based on scalar and functional predictors when moderators of treatment effect are unknown. Journal of the Royal Statistical Society Series C: Applied Statistics, 67(5):1331 1356, 2018.

Adam Ciarleglio, Eva Petkova, and Ofer Harel. Elucidating age and sex-dependent association between frontal eeg asymmetry and depression: An application of multiple imputation in functional regression. Journal of the American Statistical Association, 117 (537):12 26, 2022.

Aurore Delaigle and Peter Hall. Defining probability density for a distribution of random functions. The Annals of Statistics, 38(2):1171 1193, 2010.

Tan, Huang, Zhang and Yin

Peng Ding, Avi Feller, and Luke Miratrix. Decomposing treatment effect variation. Journal of the American Statistical Association, 114(525):304 317, 2019.

Alexander D Amour, Peng Ding, Avi Feller, Lihua Lei, and Jasjeet Sekhon. Overlap in observational studies with high-dimensional covariates. Journal of Econometrics, 221(2): 644 654, 2021.

Fr ed eric Ferraty and Philippe Vieu. Nonparametric functional data analysis: theory and practice, volume 76. Springer, 2006.

Fr ed eric Ferraty, Ali Laksaci, Amel Tadj, and Philippe Vieu. Rate of uniform consistency for nonparametric estimates with functional variables. Journal of Statistical planning and inference, 140(2):335 352, 2010.

C. Fong, C. Hazlett, and K. Imai. Covariate balancing propensity score for a continuous treatment: Application to the efficacy of political advertisements. The Annals of Applied Statistics, 12:156 177, 2018.

Antonio F Galvao and Liang Wang. Uniformly semiparametric efficient estimation of treatment effects with a continuous treatment. Journal of the American Statistical Association, 110(512):1528 1542, 2015.

Wenchuan Guo, Xiao-Hua Zhou, and Shujie Ma. Estimation of optimal individualized treatment rules using a covariate-specific treatment effect curve with high-dimensional covariates. Journal of the American Statistical Association, 116(533):309 321, 2021.

Peter Hall and Joel L Horowitz. Methodology and convergence rates for functional linear regression. The Annals of Statistics, 35(1):70 91, 2007.

Lars Peter Hansen, John Heaton, and Amir Yaron. Finite-sample properties of some alternative gmm estimators. Journal of Business & Economic Statistics, 14(3):262 280, 1996.

Keisuke Hirano and Guido W Imbens. The propensity score with continuous treatments. Applied Bayesian Modeling and Causal Inference from Incomplete-data Perspectives, 226164: 73 84, 2004.

Keisuke Hirano, Guido W Imbens, and Geert Ridder. Efficient estimation of average treatment effects using the estimated propensity score. Econometrica, 71(4):1161 1189, 2003.

K. Imai and D. A. van Dyk. Causal inference with general treatment regimes: Generalizing the propensity score. Journal of the American Statistical Association, 99:854 866, 2004.

Kosuke Imai and Marc Ratkovic. Covariate balancing propensity score. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):243 263, 2014.

Masaaki Imaizumi and Kengo Kato. A simple method to construct confidence bands in functional linear regression. Statistica Sinica, 29(4):2055 2081, 2019.

Guido W Imbens, Richard H Spady, and Phillip Johnson. Information theoretic approaches to inference in moment condition models. Econometrica, 66(2):333, 1998.

Functional Treatment Effect

Edward H Kennedy. Nonparametric causal effects based on incremental propensity score interventions. Journal of the American Statistical Association, 114(526):645 656, 2019.

Edward H Kennedy, Zongming Ma, Matthew D Mc Hugh, and Dylan S Small. Nonparametric methods for doubly robust estimation of continuous treatment effects. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(4):1229 1245, 2017.

Eric B Laber and Ana-Maria Staicu. Functional feature construction for individualized treatment regimes. Journal of the American Statistical Association, 113(523):1219 1227, 2018.

Eric B Laber and Ying-Qi Zhao. Tree-based methods for individualized treatment regimes. Biometrika, 102(3):501 514, 2015.

Fan Li and Fan Li. Propensity score weighting for causal inference with multiple treatments. The Annals of Applied Statistics, 13(4):2389 2415, 2019.

Wenbo V Li and Q-M Shao. Gaussian processes: inequalities, small ball probabilities and applications. Handbook of Statistics, 19:533 597, 2001.

Xinyi Li and Michael R Kosorok. Functional individualized treatment regimes with imaging features. ar Xiv preprint ar Xiv:2304.13003, 2023.

Yehua Li, Naisyin Wang, and Raymond J Carroll. Generalized functional linear models with semiparametric single-index interactions. Journal of the American Statistical Association, 105(490):621 633, 2010.

Zhenhua Lin and Fang Yao. Functional regression on the manifold with contamination. Biometrika, 108(1):167 181, 2021.

Zhexiao Lin, Peng Ding, and Fang Han. Estimation based on nearest neighbor matching: from density ratio to average treatment effect. Econometrica, 91(6):2187 2217, 2023.

Nengxiang Ling and Philippe Vieu. Nonparametric modelling for functional data: selected survey and tracks for future. Statistics, 52(4):934 949, 2018.

Michael J Lopez and Roee Gutman. Estimation of causal effects with multiple treatments: a review and new ideas. Statistical Science, pages 432 454, 2017a.

Michael J Lopez and Roee Gutman. Estimating the average treatment effects of nutritional label use using subclassification with regression adjustment. Statistical methods in medical research, 26(2):839 864, 2017b.

Daniel J Luckett, Eric B Laber, Anna R Kahkoska, David M Maahs, Elizabeth Mayer Davis, and Michael R Kosorok. Estimating dynamic treatment regimes in mobile health using v-learning. Journal of the American Statistical Association, 2019.

Xinwei Ma and Jingshen Wang. Robust inference using inverse probability weighting. Journal of the American Statistical Association, 115(532):1851 1860, 2020.

Tan, Huang, Zhang and Yin

Ian W Mc Keague and Min Qian. Estimation of treatment policies based on functional predictors. Statistica Sinica, 24(3):1461, 2014.

Erica EM Moodie and David A Stephens. Estimation of dose response functions for longitudinal data using the generalised propensity score. Statistical methods in medical research, 21(2):149 166, 2012.

Erica EM Moodie, Nema Dean, and Yue Ru Sun. Q-learning: Flexible learning about useful utilities. Statistics in Biosciences, 6:223 243, 2014.

Jeffrey S Morris. Functional regression. Annual Review of Statistics and Its Application, 2: 321 359, 2015.

Hans-Georg M uller and Fang Yao. Functional additive models. Journal of the American Statistical Association, 103(484):1534 1544, 2008.

Whitney K Newey. Convergence rates and asymptotic normality for series estimators. Journal of Econometrics, 79(1):147 168, 1997.

Art B Owen. Empirical likelihood ratio confidence intervals for a single functional. Biometrika, 75(2):237 249, 1988.

Hyung Park, Eva Petkova, Thaddeus Tarpey, and R Todd Ogden. Functional additive models for optimizing individualized treatment rules. Biometrics, 79(1):113 126, 2023.

James M Robins. Marginal structural models versus structural nested models as tools for causal inference. In Statistical Models in Epidemiology, the Environment, and Clinical Trials, pages 95 133. Springer, 2000.

James M Robins, Miguel Angel Hernan, and Babette Brumback. Marginal structural models and causal inference in epidemiology. Epidemiology, 11(5):550 560, 2000.

P. R. Rosenbaum and D. B. Rubin. The central role of the propensity score in observational studies for causal effects. Biometrika, 70:45 55, 1983.

Phillip J Schulte, Anastasios A Tsiatis, Eric B Laber, and Marie Davidian. Q-and a-learning methods for estimating optimal dynamic treatment regimes. Statistical Science: a Review Journal of the Institute of Mathematical Statistics, 29(4):640, 2014.

Hyejin Shin. Partial functional linear regression. Journal of Statistical Planning and Inference, 139(10):3405 3418, 2009.

Zhiqiang Tan. Model-assisted inference for treatment effects using regularized calibrated estimation with high-dimensional data. The Annals of Statistics, 48(2):811 837, 2020.

Paul Tseng and Dimitri P Bertsekas. Relaxation methods for problems with strictly convex costs and linear constraints. Mathematics of Operations Research, 16(3):462 481, 1991.

Jane-Ling Wang, Jeng-Min Chiou, and Hans-Georg M uller. Functional data analysis. Annual Review of Statistics and its Application, 3:257 295, 2016.

Functional Treatment Effect

Jiayi Wang, Raymond KW Wong, Xiaoke Zhang, and Kwun Chuen Gary Chan. Flexible functional treatment effect estimation. ar Xiv preprint ar Xiv:2309.08039, 2023.

Xiaoke Zhang and Jane-Ling Wang. From sparse to dense functional data and beyond. The Annals of Statistics, 44(5):2281 2321, 2016.

Xiaoke Zhang, Wu Xue, and Qiyue Wang. Covariate balancing functional propensity score for functional treatments in cross-sectional observational studies. Computational Statistics & Data Analysis, 163:107303, 2021.

Ying-Qi Zhao, Donglin Zeng, Eric B Laber, and Michael R Kosorok. New statistical learning methods for estimating optimal dynamic treatment regimes. Journal of the American Statistical Association, 110(510):583 598, 2015.

Hang Zhou, Fang Yao, and Huiming Zhang. Functional linear regression for discretely observed data: from ideal to reality. Biometrika, 110(2):381 393, 2023.