Causal Effect of Functional Treatment

Our functional stabilized weight (FSW) estimator does not rely on any model assumption other than (3.1). In particular, equation (2.2) suggests that we can estimate $E\{Y^{*}(z)\}$ in (3.1) using the classical functional linear regression technique (e.g., Hall and Horowitz, 2007) with the response variable $Y$ replaced by $\widehat{\pi}(Z,\mathbf{X})Y$, provided an estimator $\widehat{\pi}$ of $\pi_{0}$. To facilitate the development, suppose that an estimator $\widehat{\pi}$ of $\pi_{0}$ is available for now. Let $G(t,s)=\operatorname{Cov}\{Z(t),Z(s)\}$ be the covariance function of $Z$. Its eigenvalues $\lambda_{1}\geq\lambda_{2}\geq\cdots\geq 0$ and the corresponding orthonormal eigenfunctions $\{\phi_{j}\}_{j=1}^{\infty}$ are defined as follows,

where $\{\phi_{j}\}_{j=1}^{\infty}$, also called the PC basis, is a complete basis of $L^{2}(\mathcal{T})$. According to the Karhunen–Loéve expansion, we have

where $\mu_{Z}(t)=E\{Z(t)\}$ and $\xi_{j}=\int_{\mathcal{T}}\{Z(t)-\mu_{Z}(t)\}\phi_{j}(t)\,dt$.

With no ambiguity, let $G$ denote the linear operator such that $(G\circ\psi)(t)=\int_{\mathcal{T}}G(t,s)\psi(s)\,ds$, for a function $\psi\in L^{2}(\mathcal{T})$. Combining (2.2), (3.1) and (3.3), we can write

where $\mu_{Y,\pi}=E\{\pi_{0}(Z,\mathbf{X})Y\}=a+\int_{\mathcal{T}}b(t)\mu_{Z}(t)\,dt$. It follows that $(G\circ b)(t)=e(t)$, where $e(t)=E[\{\pi_{0}(Z,\mathbf{X})Y-\mu_{Y,\pi}\}\{Z(t)-\mu_{Z}(t)\}]$. Noting that $\{\phi_{j}\}_{j=1}^{\infty}$ is an orthonormal complete basis of $L^{2}(\mathcal{T})$, we can write $b(t)=\sum_{j=1}^{\infty}b_{j}\phi_{j}(t)$ and $e(t)=\sum_{j=1}^{\infty}e_{j}\phi_{j}(t)$, and combining with (3.2), we obtain for all $j$ that

The above argument suggests that we can estimate $b(t)$ by $\widehat{b}_{\mathrm{FSW}}(t)=\sum_{j=1}^{q}\widehat{b}_{\mathrm{FSW},j}\widehat{\phi}_{j}(t)$ with $q$ being a truncation parameter and $\widehat{b}_{\mathrm{FSW},j}=\widehat{\lambda}_{j}^{-1}\widehat{e}_{\mathrm{FSW},j}$, where the $\widehat{\lambda}_{j}$’s and $\widehat{\phi}_{j}$’s are the eigenvalues and the eigenfunctions of the empirical covariance function $\widehat{G}(s,t)=\sum_{i=1}^{n}\{Z_{i}(s)-\widehat{\mu}_{Z}(s)\}\{Z_{i}(t)-\widehat{\mu}_{Z}(t)\}/n$, respectively, and $\widehat{\mu}_{Z}(t)=\sum_{i=1}^{n}Z_{i}(t)/n$ is the empirical estimator of $\mu_{Z}(t)$. Also, the estimator of $e(t)$ is given by

with $\widehat{\mu}_{Y,\pi}=\sum_{i=1}^{n}Y_{i}\widehat{\pi}(Z_{i},\mathbf{X}_{i})/n$. We then define $\widehat{e}_{\mathrm{FSW},j}=\int_{\mathcal{T}}\widehat{e}_{\mathrm{FSW}}(t)\widehat{\phi}_{j}(t)\,dt$. Finally, the estimator of $a$ is given by $\widehat{a}_{\mathrm{FSW}}=\widehat{\mu}_{Y,\pi}-\int_{\mathcal{T}}\widehat{b}_{\mathrm{FSW}}(t)\widehat{\mu}_{Z}(t)\,dt$.

It remains to estimate the FSW, $\pi_{0}(Z,\mathbf{X})=f_{\mathbf{X}}(\mathbf{X})/f_{\mathbf{X}|Z}(\mathbf{X}|Z)$. In particular, we need to estimate the weights at the sample points, $\pi_{0}(Z_{i},\mathbf{X}_{i})$, for $i=1,\ldots,n$. A naive approach would first estimate the densities $f_{\mathbf{X}}$ and $f_{\mathbf{X}|Z}$ and then take the ratio, which, however, may lead to unstable estimates because an estimator of the denominator $f_{\mathbf{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 $\pi_{0}$ as a whole and develop a robust nonparametric estimator. The idea of estimating the density ratio directly has also been exploited in Ai et al. (2021) for a scalar treatment variable. However, their method is not applicable to the functional treatment.

Specifically, when the treatment $Z$ was a scalar, Ai et al. (2021) found that the moment equation

holds for any integrable functions $v(\mathbf{X})$ and $u(Z)$ and it identifies $\pi_{0}(\cdot,\cdot)$. Denoting the support of $\mathbf{X}$ by $\mathcal{X}$, they further estimated the function $\pi_{0}(\cdot,\cdot):\mathbb{R}\times\mathcal{X}\to\mathbb{R}$ by maximizing a generalized empirical likelihood, subject to the restriction of the empirical counterpart of (3.4), where the spaces of integrable functions $v(\mathbf{X})$ and $u(Z)$ are approximated by a growing number of basis functions; see Section 4 of Ai et al. (2021) for details. However, those restrictions are not applicable in our context where $Z$ is a random function in $L^{2}(\mathcal{T})$. In particular, it is unclear whether (3.4) still identifies $\pi_{0}(\cdot,\cdot):L^{2}(\mathcal{T})\times\mathcal{X}\to\mathbb{R}$. Moreover, the space of functionals $\{u(\cdot):L^{2}(\mathcal{T})\to\mathbb{R}\}$ is generally nonseparable, which impedes a consistent approximation of $\pi_{0}(\cdot,\cdot)$ using the basis expansion.

To overcome the challenge of estimating $\pi_{0}$, we derive a novel expanding set of equations that can identify $\pi_{0}$ from functional data and avoid approximating a general functional $u(\cdot):L^{2}(\mathcal{T})\to\mathbb{R}$. Specifically, instead of estimating the function $\pi_{0}(\cdot,\cdot):L^{2}(\mathcal{T})\times\mathcal{X}\to\mathbb{R}$, we estimate its projection $\pi_{0}(z,\cdot):\mathcal{X}\to\mathbb{R}$ for every fixed $z\in L^{2}(\mathcal{T})$, and find that

holds for any integrable function $v(\mathbf{X})$. In particular, the following proposition shows that (3.5) indeed identifies the function $\pi_{0}(z,\cdot)$ for every fixed $z\in L^{2}(\mathcal{T})$.

The proof is given in Supplement A. This result indicates that $\pi_{0}$ is fully characterized by (3.5), and thus we can use it to define our estimator. Since (3.5) is defined for a fixed $z$ and we need to estimate $\pi_{0}$ at all sample points $Z_{i}$, we consider the leave-one-out index set $\mathcal{S}_{-i}=\{j:j\neq i,j=1,\ldots,n\}$ to estimate $\pi_{0}$ at each $Z_{i}$. Specifically, for a given sample point $z=Z_{i}$, the right-hand side of (3.5) can be estimated by its empirical version $\sum_{j\in\mathcal{S}_{-j}}v(\mathbf{X}_{j})/(n-1)$, and we propose to estimate the left-hand side of (3.5) by the Nadaraya–Watson estimator for a functional covariate (see e.g., Ferraty and Vieu, 2006),

where $d(\cdot,\cdot)$ denotes a semi-metric in $L^{2}(\mathcal{T})$, $h>0$ is a bandwidth and $K$ is a kernel function quantifying the proximity of $Z_{j}$ and $Z_{i}$. Commonly used choices for $d$ include the $L^{2}$ norm and the projection metric $d(z_{1},z_{2})=[\int_{\mathcal{T}}\{\Pi(z_{1})-\Pi(z_{2})\}^{2}\,dt]^{1/2}$, where $\Pi$ denotes a projection operation to a subspace of $L^{2}(\mathcal{T})$. The projection metric may yield better asymptotic properties; see Section 4 for more details, and see Section 5 on how to choose $h$ and $K$.

As it is not possible to solve (3.5) 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 $\nu_{k}(\mathbf{X})=\big{(}v_{1}(\mathbf{X}),\ldots,v_{k}(\mathbf{X})\big{)}^{\top}$ 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 $\pi_{0}(Z_{i},\mathbf{X}_{j})$, for $j\in\mathcal{S}_{-i}$, as the solution to the $k$ equations,

which asymptotically identifies $\pi_{0}$ as $k\rightarrow\infty$ and $h\rightarrow 0$. However in practice, a finite number $k$ of equations as in (3.6) are not able to identify a unique solution. Indeed, for any strictly increasing and concave function $\rho$, replacing $\pi(Z_{i},\mathbf{X}_{j})$ by $\rho^{\prime}\{\widehat{\eta}_{Z_{i}}^{\top}\nu_{k}(\mathbf{X}_{j})\}$ would satisfy (3.6), where $\rho^{\prime}(\cdot)$ is the first derivative of $\rho(\cdot)$ and $\widehat{\eta}_{Z_{i}}$ is a $k$-dimensional vector maximizing the strictly concave function $H_{hk,Z_{i}}$ defined as

This can be verified by noting that the gradient $\nabla H_{hk,Z_{i}}(\widehat{\eta}_{Z_{i}})=0$ corresponds to (3.6) with $\pi(Z_{i},\mathbf{X}_{j})$ replaced by $\rho^{\prime}\{\widehat{\eta}_{Z_{i}}^{\top}\nu_{k}(\mathbf{X}_{j})\}$. Therefore, for a given strictly increasing and concave function $\rho$ together with $\widehat{\eta}_{Z_{i}}$ defined above, we define our estimator of $\pi_{0}(Z_{i},\mathbf{X}_{j})$ as $\widehat{\pi}_{hk}(Z_{i},\mathbf{X}_{j})=\rho^{\prime}\{\widehat{\eta}_{Z_{i}}^{\top}\nu_{k}(\mathbf{X}_{j})\}$, which also induces the estimator $\widehat{\pi}_{hk}(Z_{i},\mathbf{X}_{i})=\rho^{\prime}\{\widehat{\eta}_{Z_{i}}^{\top}\nu_{k}(\mathbf{X}_{i})\}$ by replacing $\mathbf{X}_{j}$ by $\mathbf{X}_{i}$.

To gain more insight on the function $\rho$, we investigate the generalized empirical likelihood interpretation of our estimator. As shown in Supplement B, the estimator defined by (3.7) is the dual solution to the following local generalized empirical likelihood maximization problem. For each $Z_{i}$, $i=1,\ldots,n$,

where $D(w)$ is a distance measure between $w$ and 1 for $w\in\mathbb{R}$. The function $D(v)$ is continuously differentiable satisfying $D(1)=0$ and $\rho(-u)=D\{D^{{}^{\prime}(-1)}(u)\}-u\cdot D^{{}^{\prime}(-1)}(u)$, where $D^{{}^{\prime}(-1)}$ is the inverse function of $D^{{}^{\prime}}$, the first derivative of $D$. The sample moment restriction stabilizes our weight estimator.

Our estimator $\widehat{\pi}_{hk}$ is the desired weight closest to the baseline uniform weight 1 locally around a small neighbourhood of each $Z_{i}$, subject to the finite sample version of the moment restriction in Proposition 3.1. The uniform weight can be considered as a baseline because when $Z$ and $Y^{*}(z)$ are unconditionally independent for all $z\in L^{2}(\mathcal{T})$ (i.e., a randomized trial without confoundedness), the functional stabilized weights should be uniformly equal to 1. Different choices of $\rho$ correspond to different distance measures. For example, if we take $\rho(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 $\rho(v)=\log(v)+1$ with $D(v)=-\log(v)$, the empirical likelihood (Owen, 1988), and $\rho(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). We suggest to choose $\rho(v)=-\exp(-v-1)$ which guarantees a positive weight estimator $\widehat{\pi}_{hk}(Z_{i},\mathbf{X}_{i})$.

To estimate the ADRF based on the outcome regression (OR) identification method in (2.1), an estimator of $E(Y|\mathbf{X},Z)$ is required. Thus, based on the functional linear model in (3.1), we further assume the following partially linear model,

where $a$ and $b$ are the same as those in (3.1), $\theta\in\mathbb{R}^{p}$ is the slope coefficient for $\mathbf{X}$ and $\epsilon$ is the error variable with $E(\epsilon|Z,\mathbf{X})=0$ and $E(\epsilon^{2}|Z,\mathbf{X})<\infty$ uniformly for all $Z$ and $\mathbf{X}$. Without loss of generality, we assume that $E(\theta^{\top}\mathbf{X})=0$; otherwise the non-zero mean can be absorbed into $a$.

Recalling $E\{Y^{*}(z)\}=E\{E(Y|\mathbf{X},Z=z)\}$ from (2.1), we can see that (3.9) implies (3.1). Indeed, model (3.9) imposes a stronger structure than model (3.1). For example, the linear part $\theta^{\top}\mathbf{X}$ can be replaced by any nonparametric function of $\mathbf{X}$, which still implies model (3.1). However, model (3.9) gives a simpler estimator of $E\{Y^{*}(z)\}$ with a faster convergence rate and a better interpretability of the effects of the treatment and the confounding variables on the outcome.

It remains to estimate the unknown parameters $a$ and $b$ in (3.9), for which the backfitting algorithm (Buja et al., 1989) can be adapted to our context. In the first step, we set $\theta$ 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 (3.9). Specifically, we estimate $b(t)$ initially by $\widehat{b}_{\mathrm{ini}}(t)=\sum_{j=1}^{q}\widehat{b}_{\mathrm{ini},j}\widehat{\phi}_{j}(t)$ with $\widehat{b}_{\mathrm{ini},j}=\widehat{\lambda}_{j}^{-1}\widehat{e}_{\mathrm{ini},j}$, where the $\widehat{\lambda}_{j}$’s and $\widehat{\phi}_{j}$’s are the same as those in Section 3.1, while $\widehat{e}_{\mathrm{ini}}(t)=\sum_{i=1}^{n}(Y_{i}-\widehat{\mu}_{Y})\{Z_{i}(t)-\widehat{\mu}_{Z}(t)\}/n$ with $\widehat{\mu}_{Y}=\sum_{i=1}^{n}Y_{i}/n$. We then define $\widehat{e}_{\mathrm{ini},j}=\int_{\mathcal{T}}\widehat{e}_{\mathrm{ini}}(t)\widehat{\phi}_{j}(t)\,dt$. Finally, the initial estimator $\widehat{a}_{\mathrm{ini}}$ of $a$ is defined as $\widehat{\mu}_{Y}-\int_{\mathcal{T}}\widehat{b}_{\mathrm{ini}}(t)\widehat{\mu}_{Z}(t)\,dt$.

In the second step, we use the traditional least square method to regress the residual $Y_{i}-\widehat{a}_{\mathrm{ini}}-\int_{\mathcal{T}}\widehat{b}_{\mathrm{ini}}(t)Z_{i}(t)\,dt$ on $\mathbf{X}_{i}$, for $i=1,\ldots,n$. Specifically, we estimate $\theta$ by $\widehat{\theta}=(\Xi^{\top}\Xi)^{-1}\Xi^{\top}\mathbf{Y}_{\mathrm{res}}$, where $\Xi$ is an $n\times p$ matrix with the $(i,j)$-th entry being $X_{ij}$, the $j$-th element of $\mathbf{X}_{i}$, and

Subsequently, we can repeat the first step with $Y_{i}$ replaced by $Y_{i}-\widehat{\theta}^{\top}\mathbf{X}_{i}$ to update the estimators of $\widehat{a}_{\mathrm{ini}}$ and $\widehat{b}_{\mathrm{ini}}$. This procedure is iterated until the outcome regression (OR) estimators of $a$, $b$ and $\theta$, denoted by $\widehat{a}_{\mathrm{OR}}$, $\widehat{b}_{\mathrm{OR}}$ and $\widehat{\theta}_{\mathrm{OR}}$, are stabilized.

The OR estimator is easy to implement but requires strong parametric assumptions, while the FSW estimator requires less modelling assumptions but subject to a slow convergence rate. It is thus of interest to develop an estimator by combining these two that possesses more attractive properties.

Note that (2.3) satisfies the so-called doubly robust (DR) property: for two generic functions $\widetilde{E}(Y|\mathbf{X},Z)$ and $\widetilde{\pi}(\mathbf{X},Z)$, it holds that

if either $\widetilde{E}(Y|\mathbf{X},Z)=E(Y|\mathbf{X},Z)$ or $\widetilde{\pi}(\mathbf{X},Z)=\pi_{0}(\mathbf{X},Z)$ but not necessarily both are satisfied; see Supplement C for the derivation.

Under model (3.1), (3.10) equals $a+\int_{\mathcal{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|\mathbf{X},Z)\}\pi(\mathbf{X},Z)+E_{\mathbf{X}}\{E(Y|\mathbf{X},Z)\}$ on $Z$. Specifically, we estimate $E(Y|\mathbf{X},Z)$ by the OR estimator, $\widehat{E}(Y|\mathbf{X},Z)=\widehat{a}_{\mathrm{OR}}+\int_{\mathcal{T}}\widehat{b}_{\mathrm{OR}}(t)Z(t)\,dt+\widehat{\theta}_{\mathrm{OR}}^{\top}\mathbf{X}$, and $\pi_{0}(\mathbf{X},Z)$ by $\widehat{\pi}_{hk}(\mathbf{X},Z)$. We define the DR estimator of $b$ as $\widehat{b}_{\mathrm{DR}}(t)=\sum_{j=1}^{q}\widehat{b}_{\mathrm{DR},j}\widehat{\phi}_{j}(t)$, where $\widehat{b}_{\mathrm{DR},j}=\widehat{\lambda}_{j}^{-1}\widehat{e}_{\mathrm{DR},j}$ with $\widehat{e}_{\mathrm{DR},j}=\int_{\mathcal{T}}\widehat{e}_{\mathrm{DR}}(t)\widehat{\phi}_{j}(t)\,dt$. Here, $\widehat{e}_{\mathrm{DR}}(t)$ is defined as

where

Also, the DR estimator of $a$ is defined as $\widehat{a}_{\mathrm{DR}}=\widehat{\mu}_{Y,\mathrm{DR}}-\int_{\mathcal{T}}\widehat{b}_{\mathrm{DR}}\widehat{\mu}_{Z}(t)\,dt$.

Provided that model (3.1) is correctly specified, according to the DR property (3.10), the DR estimator $\widehat{E}_{\mathrm{DR}}\{Y^{*}(z)\}=\widehat{a}_{\mathrm{DR}}+\int_{\mathcal{T}}\widehat{b}_{\mathrm{DR}}(t)z(t)\,dt$ is consistent if either $\widehat{E}(Y|\mathbf{X},Z)$ or $\widehat{\pi}_{hk}(\mathbf{X},Z)$ is consistent. The consistency of $\widehat{E}(Y|\mathbf{X},Z)$ mainly depends on correct specification of the partially linear model (3.9), while the consistency of $\widehat{\pi}_{hk}(\mathbf{X},Z)$, as a nonparametric estimator, relies on less restrictive (but technical) assumptions; see Section 4 for details.