# Causal Effect Estimation for Multivariate Continuous Treatments

Juan Chen  Key Laboratory of Advanced Theory and Application in Statistics and
Data Science-MOE, School of Statistics, East China Normal University.
Yingchun Zhou Corresponding author<EMAIL_ADDRESS>Key Laboratory of
Advanced Theory and Application in Statistics and Data Science-MOE, School of
Statistics, East China Normal University.

###### Abstract

Causal inference is widely used in various fields, such as biology, psychology
and economics, etc. In observational studies, we need to balance the
covariates before estimating causal effect. This study extends the one-
dimensional entropy balancing method to multiple dimensions to balance the
covariates. Both parametric and nonparametric methods are proposed to estimate
the causal effect of multivariate continuous treatments and theoretical
properties of the two estimations are provided. Furthermore, the simulation
results show that the proposed method is better than other methods in various
cases. Finally, we apply the method to analyze the impact of the duration and
frequency of smoking on medical expenditure. The results show that the
frequency of smoking increases medical expenditure significantly while the
duration of smoking does not.

Keywords: causal inference; causal effect; multivariate continuous treatments;
entropy balancing

##  1 INTRODUCTION

For decades, causal inference has been widely used in many fields, such as
biology, psychology and economics, etc. Most of the current research is based
on univariate treatment (binary treatment, multivalued treatment, continuous
treatment) (Imai and Ratkovic, (2014); Zhu et al., (2015); Fong et al.,
(2018); Zubizarreta, (2015); Chan et al., (2016); Xiong et al., (2017); Yiu
and Su, (2018); Dong et al., (2021); Hsu et al., (2020)). Some research are
focused on multivariate categorical treatments, such as factorial designs and
conjoint analysis to estimate the main or interaction effect of any
combination level of treatments (Hainmueller et al., (2014); Dasgupta et al.,
(2015)). However, sometimes decision-makers are interested in the causal
effects of multivariate continuous treatments in real life. For example, when
considering the impact of the export and import volume on a country’s GDP, one
is interested in a bivariate continuous treatment. The methods for
multivariate categorical treatments are not suitable to multivariate
continuous treatments and there have been few research on this. The goal of
this paper is to develop a new method to estimate the causal effect function
for multivariate continuous treatments.

A major challenge for inferring the causal effect in observational studies is
to balance the confounding covariates, which affect both the treatment and
outcome variables. The covariate balancing propensity score method is widely
used in controlling for confounding (Rosenbaum and Rubin, (1983); Rosenbaum
and Rubin, (1984); Rosenbaum and Rubin, (1985); Robins et al., (2000); Hirano
and Imbens, (2004)). When using the parametric method to model the propensity
score, the estimation bias will be large if the model is mis-specified.
Therefore, some nonparametric methods for estimating the propensity score have
been proposed, such as the kernel density estimation (Robbins et al., (2020)).
In addition, in recent years, some studies have used optimized weighting
methods to directly optimize the balance of covariates (Hainmueller, (2012);
Imai and Van Dyk, (2004); Vegetabile et al., (2020)). These methods avoid the
direct construction of the propensity scores, therefore the obtained estimate
achieves higher robustness. One of the methods, the entropy balancing method,
has been established as being doubly robust, in that a consistent estimate can
still be obtained when one of the two models, either the treatment assignment
model or the outcome model, is correctly specified (Zhao and Percival,
(2017)). Furthermore, this method can be easily implemented by solving a
convex optimization problem. Therefore, we extend this method to multivariate
continuous treatments to balance the covariates in this study.

This study has the following three contributions: First, it extends the
univariate entropy balancing method to multivariate continuous treatments to
balance the covariates. Second, both parametric and nonparametric causal
effect estimation methods for multivariate continuous treatments are proposed.
Under the parametric framework, a weighted optimization estimation is defined
and its theoretical properties are provided. Under the nonparametric
framework, B-splines are used to approximate the causal effect function and
the convergence rate of the estimation is provided. Third, we apply the
proposed method to explore the impact of the duration and frequency of smoking
on medical expenditure. The results reveal their relationship, in particular,
the frequency of smoking increases medical costs significantly.

The remainder of this paper is organized as follows: In Section 2, we
introduce the motivating example in this study. In Section 3, the entropy
balancing for multivariate treatment (EBMT) method is proposed. In Section 4,
theoretical properties of the parametric and nonparametric estimation are
shown. In section 5, methods for variance estimation and confidence interval
construction are provided. In section 6, a numerical simulation is performed
to evaluate the properties of the EBMT method. In Section 7, the EBMT method
is applied to the real data analysis. The conclusions and discussions are
summarized in Section 8.

##  2 MOTIVATING EXAMPLE

In this section, we introduce an observational study that contains two-
dimensional treatments that motivates our methodology. The causal relationship
between smoking and medical expenditure has long been a hot topic of research.
Most studies are confined to a binary treatment (smoker/non-smoker), and
estimate its effect on medical costs (Rubin, (2000); Larsen, (1999); Zeger et
al., (2000) ). However, the binary treatment is too rough in describing the
status of a smoker. A better way is to regard the duration and frequency of
smoking as bivariate continuous treatments and study their causal effect on
medical costs, which motivates this study.

The data we used is extracted from the 1987 National Medical Expenditure
Survey (NMES), which is originally studied by Elizabeth et al., (2003). The
origional study does not directly estimate the causal effect of smoking on
medical expenditure, instead, it first estimates the effect of smoking on
certain diseases and then examines the increase of medical costs due to the
diseases. Elizabeth et al., (2003) proposed a variable, called $packyear$, to
measure the cumulative smoking that combines the duration and frequency of
smoking, which is defined as

$packyear=\frac{\text{number of cigrattes per day}}{20}\times\text{number of
years smoked.}$

Imai and Dyk, (2004) directly estimates the causal effect of smoking on
medical expenditures with a univariate continuous treatment $packyear$ and
with a bivariate continuous treatment (the duration and frequency of smoking),
respectively. However, when considering the bivariate treatment, the authors
estimate their causal effects on medical costs separately, which may lose
important information.

Motivated by this example, we develop a method to estimate the causal effect
jointly when the treatment is multivariate.

##  3 ENTROPY BALANCING FOR MULTIVARIATE TREATMENTS

In this section, we introduce the EBMT method to obtain weights that can
balance covariates and both parametric and nonparametric approaches are
developed to estimate the causal effect function.

###  3.1 Notation and assumptions

Suppose that the treatment for subject $i$ is $\textbf{T}_{i}$, whose support
is $\mathcal{T}\subset\mathcal{R}^{p}$. $\textbf{X}_{i}\in\mathcal{R}^{q}$
denotes the observed covariates, where $p,q$ denote the dimensions of the
treatments and covariates, respectively. Suppose that for each subject, there
exists a potential outcome $Y_{i}(\textbf{t})$ for all
$\textbf{t}\in\mathcal{T}$. The observed outcome is defined as
$Y_{i}=Y_{i}(\textbf{t})$ if $\textbf{T}_{i}=\textbf{t}$. Assume that a sample
of observations $\\{Y_{i},\textbf{T}_{i},\textbf{X}_{i}\\}$ for
$i\in\\{1,....,n\\}$ is independently drawn from a joint distribution
$f(Y,\textbf{T},\textbf{X})$. For notational convenience, the treatments and
covariates are assumed to be standardized.

To perform causal inference with the observational data, the following three
standard assumptions are made (Hirano and Imbens, (2004); Imai and Van Dyk,
(2004)).

Assumption 1 (Strong Ignorability):

$\textbf{T}_{i}\perp Y_{i}(\textbf{t})\mid\textbf{X}_{i}$, which means that
the treatment assignment is independent of the counterfactual outcomes,
conditional on covariates. This implies that there is no unmeasured
confounding.

Assumption 2 (Positivity):

$f_{\textbf{T}\mid\textbf{X}}(\textbf{T}_{i}=\textbf{t}\mid\textbf{X}_{i})>0$
for all $\textbf{t}\in\mathcal{T}$, and the conditional density
$f(\textbf{T}_{i}\mid\textbf{X}_{i})$ is called the generalized propensity
score (Imbens, (2000)).

Assumption 3 (SUVTA):

Assume that there is no interference among the units, which means that each
individual’s outcome depends only on their own level of treatment intensity.

The goal of this paper is to estimate the causal effect function
$\mathbb{E}(Y(\mathbf{t}))$ for a multivariate continuous treatment T based on
the above assumptions. In order to estimate the causal effect function with
observational data, we first define the stabilized weight as

$w_{i}=\frac{f(\textbf{T}_{i})}{f(\textbf{T}_{i}\mid\textbf{X}_{i})},$

then under the strong ignorability assumption, one can estimate the causal
effect function based on the stabilized weight by using parametric or
nonparametric method.

###  3.2 Entropy balancing for multivariate treatments

The entropy balancing method (Hainmueller, (2012)) is used to determine the
optimal weight for inferring causal effects. It has been used for univariate
treatment and here we extend this method to multivariate treatments.

Note that the stabilized weight

$w_{i}=\frac{f(\textbf{T}_{i})}{f(\textbf{T}_{i}\mid\textbf{X}_{i})},$ (1)

satisfies the following conditions:

$\mathbb{E}(w_{i})=\iint\frac{f(\textbf{T}_{i})}{f(\textbf{T}_{i}\mid\textbf{X}_{i})}f(\textbf{T}_{i},\textbf{X}_{i})d\textbf{T}_{i}d\textbf{X}_{i}=1,$
(2)
$\begin{split}\mathbb{E}(w_{i}\textbf{T}_{i}\textbf{X}_{i}^{{}^{\prime}})&=\iint\frac{f(\textbf{T}_{i})}{f(\textbf{T}_{i}\mid\textbf{X}_{i})}\textbf{T}_{i}\textbf{X}_{i}^{{}^{\prime}}f(\textbf{T}_{i},\textbf{X}_{i})d\textbf{T}_{i}d\textbf{X}_{i}\\\
&=\int\\{\frac{f(\textbf{T}_{i})}{f(\textbf{T}_{i}\mid\textbf{X}_{i})}\textbf{T}_{i}f(\textbf{T}_{i}\mid\textbf{X}_{i})d\textbf{T}_{i}\\}\textbf{X}_{i}^{{}^{\prime}}f(\textbf{X}_{i})d\textbf{X}_{i}\\\
&=\mathbb{E}(\textbf{T}_{i})\mathbb{E}(\textbf{X}_{i}^{{}^{\prime}})=\mathbf{0}.\end{split}$
(3)

Similarly, it can be shown that weighting with $w_{i}$ also preserves the
marginal means of $\textbf{T}_{i}$ and $\textbf{X}_{i}$, which means

$\mathbb{E}(w_{i}\textbf{T}_{i})=\mathbb{E}(\textbf{T}_{i})=\mathbf{0},$ (4)
$\mathbb{E}(w_{i}\textbf{X}_{i})=\mathbb{E}(\textbf{X}_{i})=\mathbf{0}.$ (5)

Hence we can set the sample conditions as

$\sum_{i=1}^{n}w_{i}g(\textbf{T}_{i},\textbf{X}_{i})=\mathbf{0},\
\sum_{i=1}^{n}w_{i}=1,\ w_{i}>0\ \ \ \ \forall i=1,...,n,$ (6)

where
$g(\textbf{T}_{i},\textbf{X}_{i})=[\text{vec}((\textbf{T}_{i}\textbf{X}_{i}^{{}^{\prime}}))^{{}^{\prime}},\textbf{T}_{i},\textbf{X}_{i}]^{{}^{\prime}}$,
whose dimensions are $pq+p+q$. Actually, this framework can be generalized to
multivariate treatments with categorical variables or combinations of
categorical and continuous variables, which is beyond the scope of this
article.

We now choose weights that satisfy the conditions in Equation (6), while
minimizing the Kullback-Leibler divergence:

$\displaystyle\text{min}_{w}\sum_{i=1}^{n}w_{i}log(\frac{w_{i}}{v_{i}})$

s.t.

$\displaystyle\sum_{i=1}^{n}w_{i}g(\textbf{T}_{i},\textbf{X}_{i})=\mathbf{0},\
\sum_{i=1}^{n}w_{i}=1,\ w_{i}>0\ \forall i=1,...,n.$ (7)

In Equation (7), $v_{i}$ denotes the base weights, which is equal to
$\frac{1}{n}$ in this study.

We adopt a method which is similar to the standard Lagrange multiplier
technique for solving this optimization problem. Using this method, we obtain
the weighting function in terms of the Lagrange multipliers $\gamma$ as

$w_{i}=\frac{v_{i}exp(-\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i}))}{\sum_{i=1}^{n}v_{i}exp(-\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i}))}.$
(8)

The proof can be found in Appendix A.2. Based on these weights, a new
objective function is obtaind,

$\text{min}_{\gamma}\
log(\sum_{i=1}^{n}v_{i}exp(-\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i}))).$
(9)

This new dual objective function can be optimized by using an efficient convex
optimization algorithm. Note that when the dimension of covariates or
treatments is high, the numerical algorithm might fail to find a solution.
Hence, this method is most suitable to low-dimensional treatments and
covariates. High-dimensional covariates or treatments case will be considered
in future work.

###  3.3 CAUSAL EFFECT ESTIMATION

In this subsection, both parametric and nonparametric approaches are developed
to estimate the causal effect function. A weighted optimization estimation is
defined under the parametric framework and B-splines are used to approximate
the causal effect function under the nonparametric framework.

####  3.3.1 Parametric method

The causal effect function is parametrized as $s(\textbf{t};\theta)$, we
assume that it has a unique solution $\theta^{*}\in\mathcal{R}^{J}$ (with
$J\in\mathbb{N}$) defined as

$\theta^{*}=\text{agrmin}_{\theta}\int_{\mathcal{T}}\mathbb{E}[Y(\textbf{t})-s(\textbf{t};\theta)]^{2}f_{\textbf{T}}(\textbf{t})d\textbf{t}.$
(10)

The difficulty for solving Equation (10) is that the potential outcome
$Y(\textbf{t})$ is not observed for all t. Hence, Proposition 1 is proposed to
connect the potential outcome with the observed outcome.

Proposition 1

Under Assumption 1, it can be shown that

$\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}]=\int_{\mathcal{T}}\mathbb{E}[Y(\textbf{t})-s(\textbf{t};\theta)]^{2}f_{\textbf{T}}(\textbf{t})d\textbf{t}.$
(11)

The proof of Proposition 1 can be found in Appendix A.1. Note that
$Y(\textbf{t})$ on the right hand side of Equation (11) represents the
potential outcome and $Y$ on the left hand side represents the observed
outcome. Proposition 1 indicates that by inserting $w$ on the left hand side
of Equation (11), one can represent the objective function with the potential
outcome by that with the observed outcome. Therefore, the true value
$\theta^{*}$ is also a solution for the weighted optimization problem:

$\theta^{*}=\text{argmin}_{\theta}\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}].$
(12)

This result implies that the true value $\theta^{*}$ can be identified from
the observational data. One can obtain the estimator based on the sample,
which is

$\hat{\theta}=\text{argmin}_{\theta}\sum_{i=1}^{n}\hat{w_{i}}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}.$
(13)

####  3.3.2 Nonparametric method

Suppose $\mathbb{E}(Y(\textbf{t}))=s(\textbf{t})$. In a similar manner to the
proof of Proposition 1, it can be shown that

$\mathbb{E}(wY\mid\textbf{T}=\textbf{t})=\mathbb{E}(Y(\textbf{t})).$

In this paper, B-splines are used to approximate $s(\textbf{t})$. For
notational convenience, assume that without loss of generality, all B-splines
of order $r$ are defined on an extended partition associated with a uniform
partition of $m$ knots. Following Schumaker, (2007), we denote the B-spline
basis functions on the $i$th component of T as
$B_{j}(\textbf{t}_{k})(j=1,...,M;k=1,...,p;M=m+r)$. Furthermore, define

$B_{k_{1},...,k_{p}}(\textbf{t})=\prod_{j=1}^{p}B_{k_{j}}(\textbf{t}_{j}),\forall\
1\leq k_{1}....,k_{p}\leq M.$ (14)

Let $B(\textbf{t})$ be the $Q\equiv M^{p}$-dimensional vector consisting of
all product functions of the form (14) and
$Z_{n}=(B(\textbf{T}_{1}),...,B(\textbf{T}_{n}))^{{}^{\prime}}$, then the
B-spline estimation of $s(\textbf{t})$ is determined by

$\hat{s}(\textbf{t})=B(\textbf{t})^{{}^{\prime}}\hat{\beta},$ (15)

where

$\hat{\beta}=(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}}\hat{W}Y,\
\hat{W}=diag(\hat{w}_{1},.....,\hat{w}_{n}).$

##  4 LARGE SAMPLE PROPERTIES

To establish the large sample properties of the proposed estimators in section
3.3, we show the consistency and limiting distribution of the optimization
estimator $\hat{\theta}$ under the parametric framework and the convergence
rate of the nonparametric estimator $\hat{s}(\textbf{t})$. The following
assumptions are made:

Assumption 4

  1. (i)

The support $\mathcal{T}$ of T is a compact subset of $\mathcal{R}^{p}$ and
the support $\mathcal{X}$ of X is a compact subset of $\mathcal{R}^{q}$.

  2. (ii)

There exist two positive constants $a_{1}$ and $a_{2}$, such that
$0<a_{1}\leq w\leq
a_{2}<\infty,\forall(\textbf{t},\textbf{x})\in\mathcal{T}\times\mathcal{X}$.

Assumption 5

  1. (i)

The parameter space $\Theta$ is a compact set and the true parameter
$\theta_{0}$ is in the interior of $\Theta$.

  2. (ii)

$(Y-s(\textbf{t};\theta))^{2}$ is continuous in $\theta$ ,
$\mathbb{E}[\text{sup}_{\theta}(Y-s(\textbf{t};\theta))^{2}]<\infty$ and
$\text{sup}_{\theta}\mathbb{E}[(Y-s(\textbf{t};\theta))^{4}]<\infty$.

Assumption 6

  1. (i)

$s(\mathbf{t};\theta)$ is twice continuously differentiable in
$\theta\in\Theta$ and let
$h(\mathbf{t};\theta)\equiv\bigtriangledown_{\theta}s(\mathbf{t};\theta)$.

  2. (ii)

$\mathbb{E}\\{w(Y-s(\mathbf{T};\theta))h(\mathbf{T};\theta)\\}$ is
differentiable with respect to $\theta$ and
$U\equiv-\bigtriangledown_{\theta}\mathbb{E}\\{w(Y-s(\mathbf{T};\theta))h(\mathbf{T};\theta)\\}\mid_{\theta=\theta^{*}}$
is nonsingular.

  3. (iii)

$\mathbb{E}[\text{sup}_{\theta}\mid
Y-s(\textbf{t};\theta)\mid^{2+\delta}]<\infty$ for some $\delta>0$ and there
extists some finite positive constants $a$ and $b$ such that
$\mathbb{E}[\text{sup}_{\theta_{1}:\mid\mid\theta_{1}-\theta\mid\mid<\delta_{1}}\mid
s(\textbf{t};\theta_{1})-s(\textbf{t};\theta)\mid^{2}]^{1/2}<a\cdot\delta_{1}^{b}$
for any $\theta\in\Theta$ and any small $\delta_{1}>0$.

Assumption 4(i) restricts the treatments and covariates within a bounded
range. This is convenient for calculating the convergence rate under the
$L_{\infty}$ norm. Assumption 4(ii) restricts the weight function to be
bounded away from zero. Given Assumption 4(i), Assumption 4(ii) is equivalent
to the positivity assumption. Assumption 5(i) is a commonly used assumption in
nonparametric regression. Assumption 5(ii) is an envelope condition applicable
to the uniform law of large numbers. Assumption 6 (i) and (ii) impose
sufficient regularity conditions on the causal effect function and its
derivative function. Assumption 6 (iii) is a stochastic equicontinuity
condition, which is needed for establishing weak convergence (Andrews,
(1994)). Based on these assumptions, the following theorems are established.

Theorem 1

  1. (i)

Under Assumptions 1-5, $\mid\mid\hat{\theta}-\theta^{*}\mid\mid\ \to_{p}0$.

  2. (ii)

Under Assumptions 1-6, $\sqrt{n}(\hat{\theta}-\theta^{*})\to_{d}N(0,V)$, where

$V=4U^{-1}\cdot\mathbb{E}\\{w^{2}(Y-s(\mathbf{T};\theta^{*}))^{2}h(\mathbf{T};\theta^{*})h(\mathbf{T};\theta^{*})^{{}^{\prime}}\\}\cdot
U^{-1}$

Under the nonparametric framework, Theorem 2 is established to obtain the
convergence rate of the estimate
$\hat{s}(\textbf{t})=B(\textbf{t})^{{}^{\prime}}\hat{\beta}$.

Theorem 2 Suppose $\text{sup}_{\textbf{t}\in\mathcal{T}}\mid
s(\textbf{t})-B(\textbf{t})^{{}^{\prime}}\beta^{*})\mid=O(Q^{-\tilde{\alpha}})$
holds for some $\tilde{\alpha}$ and $\beta^{*}\in\mathbb{R}^{Q}$, and
Assumptions 1-4 hold, then

$\begin{split}&\int\mid\hat{s}(\textbf{t})-s(\textbf{t})\mid^{2}f_{\textbf{T}}(\textbf{t})d\textbf{t}=O_{p}(Q^{-2\tilde{\alpha}}+\frac{Q}{n})\\\
&\text{sup}_{\textbf{t}}\mid\hat{s}(\textbf{t})-s(\textbf{t})\mid=O_{p}(Q^{-\tilde{\alpha}}+\sqrt{\frac{Q}{n}})\end{split}$

The proofs of Theorem 1 and Theorem 2 can be found in Appendix A.3.

##  5 VARIANCE ESTIMATION AND CONFIDENCE INTERVAL

###  5.1 Variance estimation

To conduct statiatical inference, a consistent estimator of the covariance
matrix $V$ is needed, which can be obtained by replacing $w$ and $\theta^{*}$
by their consistent estimators. Specifically, one can write $U$ as

$\begin{split}U&=-\bigtriangledown_{\theta}\mathbb{E}\\{w(Y-s(\mathbf{T};\theta))h(\mathbf{T};\theta)\\}\mid_{\theta=\theta^{*}}\\\
&=-\mathbb{E}\\{w\bigtriangledown_{\theta}[(Y-s(\mathbf{T};\theta))h(\mathbf{T};\theta)]\\}\mid_{\theta=\theta^{*}}\\\
&=-\mathbb{E}\\{w[-h(\mathbf{T};\theta)h(\mathbf{T};\theta)^{{}^{\prime}}+(Y-s(\mathbf{T};\theta))\bigtriangledown_{\theta}h(\mathbf{T};\theta)]\\}\mid_{\theta=\theta^{*}}\\\
&=\mathbb{E}\\{w[h(\mathbf{T};\theta^{*})h(\mathbf{T};\theta^{*})^{{}^{\prime}}-(Y-s(\mathbf{T};\theta^{*}))\bigtriangledown_{\theta}h(\mathbf{T};\theta^{*})]\\}.\end{split}$

Hence,

$\hat{U}=\frac{1}{n}\sum_{i=1}^{n}\hat{w}_{i}[h(\mathbf{T}_{i};\hat{\theta})h(\mathbf{T}_{i};\hat{\theta})^{{}^{\prime}}-(Y_{i}-s(\mathbf{T}_{i};\hat{\theta}))\bigtriangledown_{\theta}h(\mathbf{T}_{i};\hat{\theta})].$

Then we have

$\hat{V}=4\hat{U}^{-1}\cdot\\{\frac{1}{n}\sum_{i=1}^{n}\hat{w}_{i}^{2}(Y_{i}-s(\mathbf{T}_{i};\hat{\theta}))^{2}h(\mathbf{T}_{i};\hat{\theta})h(\mathbf{T}_{i};\hat{\theta})^{{}^{\prime}}\\}\cdot\hat{U}^{-1}.$

According to Theorem 1(i) and Theorem 3 in Ai et al., (2021), we have
$\mid\mid\hat{\theta}-\theta^{*}\mid\mid\to_{p}0$ and
$\text{sup}\mid\hat{w}-w\mid=o_{p}(1)$, which implies that $\hat{V}$ is
consistent.

###  5.2 Confidence interval

According to the asymptotic normality of
$\hat{\theta}=(\hat{\theta}_{1},\dots,\hat{\theta}_{J})^{{}^{\prime}}$, one
can construct the 95$\%$ confidence interval of $\hat{\theta}_{j}$ as

$[\hat{\theta}_{j}-1.96\cdot\hat{SE}_{j},\hat{\theta}_{j}+1.96\cdot\hat{SE}_{j}],\
\forall j=1,\dots,J,$ (16)

where $\hat{SE}_{j}=\hat{V}_{jj}/\sqrt{n}$ is the standard error of
$\hat{\theta}_{j}$, $\hat{V}_{jj}$ is the $(j,j)$-element of $\hat{V}$.

Alternatively, one can also construct the confidence interval using bootstrap
method. Suppose the $b^{\text{th}}$ bootstrap sample
$\\{Y_{i}^{(b)},\mathbf{T}_{i}^{b},\mathbf{X}_{i}^{(b)}\\},i=1,\dots,n;b=1,\dots,B$
is sampled with replacement from the origional sample
$\\{Y_{i},\mathbf{T}_{i},\mathbf{X}_{i}\\},i=1,\dots,n$ with uniform
distribution. Denote $\hat{\theta}^{(b)}$ as the estimator of $\theta^{*}$
based on the $b^{\text{th}}$ bootstrap sample. Arrange the $B$ bootstrap
estimators in order from smallest to largest, then the 95$\%$ confidence
interval can be defined as

$[\hat{\theta}_{j}^{(0.025B)},\hat{\theta}_{j}^{(0.975B)}],\ \forall
j=1,\dots,J.$ (17)

The difference of the two methods is that the first method (16) relies on the
asymptotic normality result while the second method (17) does not.

##  6 SIMULATION

To examine the properties of the proposed estimators under finite samples,
simulation studies are performed under different data settings. The main
motivation in designing the simulation settings is to compare the proposed
method with three other methods when the treatment assignment model and the
outcome model are linear and nonlinear in various ways.

###  6.1 Assessment criteria of covariate balance and effect estimation

Assume that the treatments follow the multiple multivariate linear regression
model:

$\mathbf{T}_{i}=B^{{}^{\prime}}\mathbf{X}_{i}+\epsilon_{i},i=1,....,n,$ (18)

where $\mathbf{T}_{i}$ and $\mathbf{X}_{i}$ denote the $p$-dimensional
treatments and $q$-dimensional covariates, respectively. $B_{q\times
p}=(\bm{\beta}_{1}:\bm{\beta}_{2}:...:\bm{\beta}_{q})$ represents the
coefficient matrix, and $\epsilon_{i}\sim N_{p}(\mathbf{0},\Sigma)$. The MLEs
of $B$ are $\hat{B}_{q\times
p}=(\mathbf{X}^{{}^{\prime}}\mathbf{X})^{-1}\mathbf{X}^{{}^{\prime}}\mathbf{T}$,
where
$\mathbf{T}=(\mathbf{T}_{1}^{{}^{\prime}},.....,\mathbf{T}_{n}^{{}^{\prime}})^{{}^{\prime}},\
\mathbf{X}=(\mathbf{X}_{1}^{{}^{\prime}},.....,\mathbf{X}_{n}^{{}^{\prime}})^{{}^{\prime}}$.

To assess the performance of the covariate balance, the following hypothesis
test is considered:

$H_{0}:\bm{\beta}_{1}=\bm{\beta}_{2}=.....=\bm{\beta}_{q}=\mathbf{0}.$

The likelihood ratio statistic is:

$\Lambda=\frac{\mid SSE\mid}{\mid SSE+SSH\mid},$

where $SSE=\mathbf{T}^{T}\mathbf{T}-\hat{B}^{T}\mathbf{X}^{T}\mathbf{T}$ and
$SSH=\hat{B}^{T}\mathbf{X}^{T}\mathbf{T}-n\bar{\mathbf{T}}\bar{\mathbf{T}}^{T}$
denote the the residual sum of squares and the predicted sum of squares,
respectively.

When $n$ is large,

$-2log(\Lambda)\to\chi_{qp}^{2}.$

Hence, it is harder to reject the null hypothesis when $-2log(\Lambda)$ is
closer to zero, which implies that the covariate balance performs well.

Next, we focus on the outcome model, which is linear and nonlinear in the
treatments, respectively.

$\displaystyle Y_{i}=r(\mathbf{X}_{i})+\mathbf{T}_{i}\beta+\varepsilon_{i},$
(19) $\displaystyle
Y_{i}=r(\mathbf{X}_{i})+s(\mathbf{T}_{i})+\varepsilon_{i},$ (20)

where $\varepsilon_{i}\sim N(0,\sigma^{2}),i=1,...,n$.

For the linear outcome model, the bias and RMSE of the coefficients are used
to assess the performance of the different methods in various cases . For the
nonparametric outcome model, we use

$\displaystyle\text{RMSE}=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(\hat{s}(\mathbf{t})_{i}-s(\mathbf{t})_{i})^{2}}$

as the assessment criteria.

###  6.2 Data generating process

The data generating process is described in this subsection. For the treatment
assignment model, two major cases are considered, which are T-2d-L, T-2d-NL,
respectively, where ”T” stands for treatment, ”2d” stands for two-dimension,
”L” and ”NL” stands for linear and nonlinear in the covariates. For the
outcome model, linear and nonlinear cases are studied, which are considered
for each treatment assignment model. For notational convenience,
”Specification $YK$” is used to denote the $K\text{th}$ outcome model.

Since the covariates are shared across all scenarios, their data generating
process is first described. Specifically, we independently draw 5 covariates
from a multivariate normal distribution with mean 0, variance 1 and covariance
0.2, that is,

$\mathbf{X}=(X_{1},....,X_{5})^{{}^{\prime}}\sim N_{5}(\mu,\Sigma)\
\text{with}\ \mu=\begin{pmatrix}0\\\ \vdots\\\ 0\end{pmatrix}\text{and}\
\Sigma=\begin{pmatrix}1&0.2&\dots&0.2\\\ 0.2&1&0.2\dots&0.2\\\ {4}\\\
0.2&0.2&\dots&1\end{pmatrix}_{5\times 5.}$

In the first simulation setting, assume that both the treatment assignment
model ”T-2d-L” and the outcome model are linear in the covariates, with true
data generating process given by

$\textbf{T-2d-L}:\mathbf{T}_{i}=B_{1}^{{}^{\prime}}\mathbf{X}_{i}+\mathbf{\epsilon}_{i},$
(21)

Specification $Y1$:

$\displaystyle
Y_{i}=T_{i1}+T_{i2}+X_{i1}+0.1X_{i2}+0.1X_{i5}+\varepsilon_{i}.$

Specification $Y2$:

$\displaystyle
Y_{i}=T_{i1}+T_{i2}+0.2T_{i1}T_{i2}+X_{i1}+0.1X_{i2}+0.1X_{i5}+\varepsilon_{i}.$

Specification $Y3$:

$\displaystyle
Y_{i}=T_{i1}+T_{i2}+(T_{i1}-T_{i2})^{2}+X_{i1}+0.1X_{i2}+0.1X_{i5}+\varepsilon_{i},$

where

$\displaystyle\mathbf{T}_{i}=(T_{i1},T_{i2})^{{}^{\prime}},\
B_{1}=\begin{bmatrix}1&1\\\ 0&0.2\\\ 0.2&0\\\ 0&0\\\ 0&0\end{bmatrix}_{5\times
2},$ $\displaystyle varepsilon_{i}\sim N(0,2^{2}),$
$\displaystyle\mathbf{\epsilon}_{i}\sim N_{2}(\mu,\Sigma)\ \text{with}\
\mu=\begin{pmatrix}0\\\ 0\end{pmatrix}\text{and}\
\Sigma=\begin{pmatrix}3&0.8\\\ 0.8&3\end{pmatrix}.$

In the second simulation setting, consider the case when the treatment
assignment model is misspecified and the outcome model remains the same as in
the first simulation setting, that is, the treatment model is nonlinear in the
covariates, denoted as ”T-2d-NL”:

$\textbf{T-2d-NL}:\mathbf{T}_{i}=B_{1}^{{}^{\prime}}\mathbf{X}_{i}+B_{2}^{{}^{\prime}}(\mathbf{X}_{i}*\mathbf{X}_{i})+\epsilon_{i},$
(22)

where $B_{1}$ and $\mathbf{\epsilon}_{i}$ are the same as in Equation (22),
$\textbf{X}*\textbf{X}$ represents the Hadamard product of matrices X and X,
and its corresponding element is
$(\textbf{X}*\textbf{X})_{ij}=(x_{ij}x_{ij})$, $B_{2}=\begin{bmatrix}1&1\\\
0&0\\\ \vdots&\vdots\\\ 0&0\end{bmatrix}_{5\times 2}$.

In the third simulation setting, we consider misspecification of the outcome
model by including nonlinear covariate term in the outcome model. In this
case, the generating process for the treatment model remains the same as in
the first simulation setting. Specifically, the misspecified outcome model is
defined as

Specification$Y4$:

$\displaystyle
Y_{i}=T_{i1}+T_{i2}+(X_{i1}+1)^{3}+0.1X_{i2}+0.1X_{i5}+\varepsilon_{i}.$

Specification $Y5$:

$\displaystyle
Y_{i}=T_{i1}+T_{i2}+0.2T_{i1}T_{i2}+(X_{i1}+1)^{3}+0.1X_{i2}+0.1X_{i5}+\varepsilon_{i}.$

Specification $Y6$:

$\displaystyle
Y_{i}=T_{i1}+T_{i2}+(T_{i1}-T_{i2})^{2}+(X_{i1}+1)^{3}+0.1X_{i2}+0.1X_{i5}+\varepsilon_{i},$

where $\varepsilon_{i}\sim N(0,2^{2})$.

The last simulation setting considers the misspecification of both the
treatment assignment model and the outcome model, which are defined in the
second simulation setting and the third simulation setting, respectively.

###  6.3 Simulation results

In this subsection, simulation results are compared between the proposed
method (EBMT), the regression covariate adjustment method (RCAM), entropy
balancing for univariate treatment method (EBUT) and the multivariate
generalized propensity score method (mvGPS) (Williams and Crespi, (2020)),
where RCAM refers to estimating the causal effect function by adjusting for
all covariates in the outcome model, EBUT refers to handling each single
treatment seperately, mvGPS refers to the method that estimates the causal
effect function assuming that the treatments follow a multivariate normal
distribution.

Specifically, the performance of covariate balancing are compared among EBMT,
mvGPS and Unweighted methods, since RCAM does not balance covariates and EBUT
does not fit multiple multivariate regression so that the balancing effect
cannot be measured. When there is no interaction effect bewteen treatments,
the estimation accuracy of all four methods is compared. But when there is
interaction between treatments, only three methods (EBMT, RCAT, mvGPS) are
compared since the method EBUT can not deal with this case. In addition, the
95$\%$ coverage probability and the average width of the confidence interval
are considered. These are obtained using the bootstrap method with $B=500$
iterations based on Equation (17) for each effect estimation. For each
numerical setting, 1000 independent simulation experiments are run.

Figure 1 shows the results for the covariate balancing of the linear two-
dimensional treatment assignment model (top) and nonlinear two-dimensional
treatment assignment model (bottom) under different settings. The statistics
$-2log(\Lambda)$ produced by the EBMT method are almost zero, while those of
the mvGPS and Unweighted methods are far away from zero. These indicate that
the EBMT method balances the covariates well and is robust when the treatment
assignment model is linear and nonlinear in the covariates.

Table 1 shows the results of the causal effect estimation for the linear
outcome model without interaction between the two-dimensional treatments under
different simulation settings. It can be seen that when both the treatment
assignment model and outcome model are correctly specified, the mean bias of
the EBMT method is the smallest among the four methods. The RMSE of the EBMT
method is the second smallest while that of RCAT is the smallest. The RMSE of
all methods decreases as the sample size increases. Among the four methods
mvGPS performs the worst. When one of the two models is misspecified, the mean
bias of EBMT is the smallest in both cases, and the RMSE of EBMT is the
smallest when the treatment model is misspecified. When both the two models
are misspecified, all methods fail.

Table 2 shows the results of the causal effect estimation for the linear
outcome model with interaction between the treatments. As can be seen, the
results for the main effect are similar to those in Table 1, that is, EBMT
performs the best in most cases. For the interaction effect, the mean biases
are slightly larger in most cases. Similarly, all methods fail when both the
two models are misspecified.

Figure 1: Performance comparison of covariate balancing for two-dimensional treatments. Table 1: Mean bias and RMSE of the coefficients in the linear outcome model without interaction between the two-dimensional treatments. Method | Both $E(T\mid X)$ and $E(Y\mid X)$ correctly specified
---|---
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=1$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | -0.002/0.086 | 0.001/0.052 | 0.005/0.188 | 0.004/0.119
RCAM | -0.005/0.058 | -0.004/0.043 | 0.011/0.104 | 0.010/0.081
EBUT | -0.003/0.084 | -0.002/0.058 | 0.008/0.206 | 0.007/0.133
mvGPS | -0.012/0.132 | 0.006/0.082 | 0.151/0.319 | 0.035/0.171
| $E(T\mid X)$ correctly specified, $E(Y\mid X)$ misspecified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=1$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | -0.008/0.125 | -0.001/0.087 | 0.019/0.276 | 0.028/0.178
RCAM | -0.011/0.151 | 0.007/0.105 | 0.031/0.333 | -0.047/0.194
EBUT | 0.014/0.189 | 0.022/0.144 | 0.023/0.284 | 0.034/0.180
mvGPS | -0.075/0.767 | -0.016/0.591 | 1.035/1.408 | 0.762/1.051
|  $E(T\mid X)$ misspecified, $E(Y\mid X)$ correctly specified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | -0.001/0.063 | 0.003/0.048 | 0.005/0.118 | -0.004/0.074
RCAM | -0.003/0.054 | 0.007/0.043 | 0.009/0.079 | -0.005/0.058
EBUT | -0.013/0.105 | 0.003/0.067 | 0.008/0.137 | -0.007/0.108
mvGPS | -0.063/0.616 | -0.015/0.321 | -0.369/0.795 | -0.178/0.482
| Both $E(T\mid X)$ and $E(Y\mid X)$ misspecified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | 0.006/0.119 | 0.009/0.090 | 1.032/1.058 | 0.965/0.994
RCAM | 0.001/0.116 | 0.006/0.085 | 2.144/2.174 | 2.139/2.154
EBUT | 0.008/0.147 | 0.020/0.100 | 1.082/1.105 | 1.080/1.091
mvGPS | -0.076/0.814 | 0.228/1.762 | -1.359/1.742 | -1.937/2.650
Table 2: Mean bias and RMSE of the coefficients in the linear outcome model with interaction between the two-dimensional treatments.  Method | Both $E(T\mid X)$ and $E(Y\mid X)$ correctly specified
---|---
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$) | $\hat{\beta}_{12}$ ($\beta_{12}^{*}=0.2$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | 0.002/0.093 | -0.004/0.064 | 0.003/0.184 | 0.002/0.137 | -0.013/0.086 | -0.011/0.062
RCAM | 0.009/0.067 | -0.007/0.43 | 0.006/0.113 | 0.009/0.081 | -0.006/0.036 | -0.005/0.022
mvGPS | 0.011/0.139 | -0.017/0.123 | 0.153/0.317 | 0.086/0.271 | -0.008/0.076 | -0.007/0.068
|  $E(T\mid X)$ correctly specified, $E(Y\mid X)$ misspecified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$) | $\hat{\beta}_{12}$ ($\beta_{12}^{*}=0.2$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | 0.008/0.111 | -0.007/0.085 | -0.005/0.283 | -0.017/0.194 | -0.056/0.170 | -0.031/0.143
RCAM | 0.024/0.151 | 0.012/0.106 | -0.044/0.306 | 0.027/0.228 | -0.589/0.604 | -0.606/0.620
mvGPS | 0.176/0.796 | -0.067/0.430 | 0.575/0.431 | 0.621/0.387 | -0.066/0.258 | -0.004/0.193
|  $E(T\mid X)$ misspecified, $E(Y\mid X)$ correctly specified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$) | $\hat{\beta}_{12}$ ($\beta_{12}^{*}=0.2$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | 0.003/0.022 | 0.002/0.016 | -0.005/0.037 | -0.003/0.025 | -0.004/0.009 | -0.002/0.008
RCAM | 0.010/0.030 | 0.009/0.028 | -0.012/0.025 | -0.008/0.018 | -0.002/0.003 | -0.001/0.004
mvGPS | 0.002/0.043 | 0.006/0.032 | 0.109/0.128 | -0.093/0.100 | -0.009/0.014 | -0.007/0.012
| Both $E(T\mid X)$ and $E(Y\mid X)$ misspecified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) |  $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$) | $\hat{\beta}_{12}$ ($\beta_{12}^{*}=0.2$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE | Bias/RMSE
EBMT | 0.093/0.177 | 0.153/0.207 | 1.435/1.507 | 1.617/1.673 | -0.104/0.186 | -0.156/0.198
RCAM | -0.263/0.278 | -0.254/0.266 | 0.951/0.967 | 0.995/1.004 | 0.266/0.269 | 0.254/0.257
mvGPS | 0.417/0.947 | 0.571/0.994 | -0.888/1.807 | -0.989/1.474 | -0.146/0.304 | -0.168/0.245
Table 3: Confidence interval estimation of the coefficients in the linear outcome model without interaction between the two-dimensional treatments .  Method | Both $E(T\mid X)$ and $E(Y\mid X)$ correctly specified
---|---
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=1$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW
EBMT | 0.951/0.325 | 0.956/0.222 | 0.946/0.681 | 0.953/0.478
RCAM | 0.926/0.222 | 0.942/0.162 | 0.934/0.438 | 0.952/0.318
EBUT | 0.944/0.0.306 | 0.953/0.226 | 0.877/0.734 | 0.914/0.527
mvGPS | 0.948/0.434 | 0.954/0.340 | 0.758/0.723 | 0.812/0.598
|  $E(T\mid X)$ correctly specified, $E(Y\mid X)$ misspecified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=1$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW
EBMT | 0.943/0.429 | 0.952/0.290 | 0.933/0.870 | 0.942/0.644
RCAM | 0.914/0.630 | 0.936/0.414 | 0.925/1.307 | 0.941/0.819
EBUT | 0.933/0.606 | 0.941/0.418 | 0.904/0.942 | 0.917/0.690
mvGPS | 0.887/1.478 | 0.912/1.087 | 0.527/2.243 | 0.634/1.822
|  $E(T\mid X)$ misspecified, $E(Y\mid X)$ correctly specified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW
EBMT | 0.943/0.263 | 0.953/0.189 | 0.948/0.381 | 0.952/0.284
RCAM | 0.923/0.221 | 0.945/0.163 | 0.927/0.299 | 0.941/0.218
EBUT | 0.945/0.339 | 0.956/0.240 | 0.914/0.523 | 0.923/0.381
mvGPS | 0.932/1.240 | 0.935/1.232 | 0.847/1.568 | 0.866/1.389
| Both $E(T\mid X)$ and $E(Y\mid X)$ misspecified
|  $\hat{\beta}_{1}$ ($\beta_{1}^{*}=1$) | $\hat{\beta}_{2}$ ($\beta_{2}^{*}=0.8$)
| $n=500$ | $n=1000$ | $n=500$ | $n=1000$
| 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW | 95 $\%$ CP/AW
EBMT | 0.914/0.379 | 0.943/0.278 | 0.029/0.717 | 0.034/0.526
RCAM | 0.961/0.436 | 0.965/0.299 | 0.001/1.169 | 0.001/0.772
EBUT | 0.886/0.454 | 0.911/0.334 | 0.012/0.726 | 0.020/0.532
mvGPS | 0.771/2.254 | 0.892/2.074 | 0.042/3.293 | 0.038/2.683
Table 4: Mean RMSE comparison of the nonparametric outcome model of the two-dimensional treatments.  Method | Both $E(T\mid X)$ and $E(Y\mid X)$ correctly specified
---|---
| $n=500$ | $n=1000$
EBMT | 0.933 | 0.828
mvGPS | 1.630 | 1.298
|  $E(T\mid X)$ misspecified, $E(Y\mid X)$ correctly specified
| $n=500$ | $n=1000$
EBMT | 3.219 | 2.937
mvGPS | 4.802 | 4.480
|  $E(T\mid X)$ misspecified, $E(Y\mid X)$ correctly specified
| $n=500$ | $n=1000$
EBMT | 1.258 | 0.968
mvGPS | 1.788 | 2.784
| Both $E(T\mid X)$ and $E(Y\mid X)$ misspecified
| $n=500$ | $n=1000$
EBMT | 6.043 | 5.703
mvGPS | 7.286 | 7.537

Table 3 shows the results of the 95$\%$ coverage probability and average width
of the bootstrap confidence interval of the causal effect estimation. It can
be seen that the 95$\%$ coverage probability of the EBMT method is nearly
identical to 0.95 and it’s larger than the other methods under different
scenarios except for the case that both the treatment assignment model and
outcome model are misspecified. For all methods, the average width of the
confidence interval shrinks as the sample size grows.

Table 4 shows the mean RMSE of the estimated causal effect function for the
nonparametric outcome model under different scenarios. Since only the EBMT and
mvGPS methods can be applied to such problem, the results of those two methods
are compared. It can be seen that EBMT method enjoys a better performance
compared with mvGPS method when both the treatment assignment model and the
outcome model are correctly specified or either one of the two models is
correctly specified. Both methods fail when both the two models are
misspecificed.

The simulation results can be summarized as follows. (i) In the case that both
the treatment assignment model and outcome model are correctly specified, EBMT
has the smallest mean bias of the main effect and slightly larger RMSE than
RCAT. (ii) In the case that the treatment assignment model is correctly
specified and the outcome model is misspecified, EBMT has the smallest mean
bias and RMSE of all effects. (iii) In the case that both the treatment
assignment model and outcome model are misspecified, all methods fail. (iv) In
the case of nonparametric outcome model, EBMT has the smallest mean RMSE.

##  7 APPLICATION

We now apply the EBMT method to the motivating example described in Section 2.
The goal is to analyze the impact of the duration and frequency of smoking on
medical expenditures.

The pre-treatment covariates $\mathbf{X}$ includes age at the time of the
survey, age when the individual started smoking, gender, race, marital status,
education level, census region, poverty status and seat belt usage, which are
the covariates used in Imai et al. (2004). The treatments we are interested in
are the duration and frequency of smoking, and the outcome is log(Total
Medical Expenditure). Similar to the setting of Imai et al. (2004), we assume
a linear outcome model and conduct a complete-case analysis by discarding all
missing data, yielding a sample of 7847 smokers.

Before estimating the causal effect, we first examine the covariate balancing
of EBMT and mvGPS methods by conducting a multiple multivariate regression.
The statistic $-2log(\Lambda)$ defined in section 6.1 is $3.95e^{-8}$ for EBMT
method, 2.27 for mvGPS method and 2.31 for Unweighted method, which shows that
EBMT balances covariates very well while mvGPS does not balance well.

The causal effect function is estimated by the parametric method assuming a
linear model of the outcome. Bootstrap method is used to obtain the standard
error and confidence interval for the parameter estimates. Each of the 1000
bootstrap replicates estimates weight using EBMT, mvGPS and EBUT methods,
since RACM does not estimate weight. Then the outcome model is estimated by
the four methods EBMT, RACM,EBUT and mvGPS.

Table 5 shows the estimated causal effect of the duration and frequency of
smoking on medical expenditure as well as its standard error and confidence
interval. All methods indicate that the duration of smoking has no significant
impact on medical expenditure since their confidence intervals all cover zero
while most methods agree that the frequency of smoking increases medical
expenditure significantly except mvGPS. Besides, the width of confidence
interval of the estimated causal effect based on EBMT method is smallest.
Compared with univariate analysis using a single variable $packyear$, the
analysis of the bivariate treatments provides more information that the
significant effect of $packyear$ attributes mostly to the frequency of smoking
rather than to its duration.

Table 5: Causal effect estimation of the duration and frequency of smoking on meadical expenditure.  Method | Duration | Frequency
---|---|---
| Estimate | SE | 95$\%$ CI | Estimate | SE | 95$\%$ CI
EBMT | 0.003 | 0.002 | (-0.003,0.004) | 0.007 | 0.002 | (0.003,0.009)
RCAM | -0.001 | 0.004 | (-0.007,0.007) | 0.006 | 0.003 | (0.001,0.011)
EBUT | 0.002 | 0.004 | (-0.009,0.006) | 0.007 | 0.002 | (0.004,0.011)
mvGPS | -0.073 | 0.169 | (-0.079,0.104) | 0.018 | 0.047 | (-0.079,0.104)

##  8 CONCLUSION AND DISCUSSION

In this study, we extend the one-dimensional entropy balancing method to
multidimensional treatments. In addition, parametric and nonparametric methods
are developed to estimate the causal effect and their theoretical properties
are provided. The simulation results show that the proposed method balances
covariates well and produces a smaller mean bias compared with other methods.
In the real data analysis, the EBMT method is applied to investigate causal
relationship between the duration and frequency of smoking on medical
expenditure, which indicates that the frequency of smoking inceases medical
costs significantly.

In this paper, we mainly consider the causal effect function
$\mathbb{E}(Y(t))$ as the estimand for continuous treatment as it is general.
The average treatment effect or average partial effect can also be considered,
but they can be easily estimated based on the estimates of causal effect
function $\mathbb{E}(Y(t+\triangle t)-Y(t))$, $\frac{\mathbb{E}Y(t+\triangle
t)-\mathbb{E}Y(t)}{\triangle t}$ (Viet, (2008)). Indeed, the causal effect
function provides a complete description of the causal effect, rather than a
summary measure.

Future works include: (1) Estimating causal effect in the case of high-
dimensional covariates or treatments; (2) Considering treatments involves
other types of complex data, such as longitudinal data and functional data;
(3) Extending this method to the multidimensional outcome variables.

## ACKNOWLEDGEMENT

The work was supported by National Natural Science Foundation of China
(project number: 11771146, 11831008), the 111 Project (B14019) and Program of
Shanghai Subject Chief Scientist (14XD1401600). The work was supported by
National Natural Science Foundation of China (project number: 11771146,
11831008), the 111 Project (B14019) and Program of Shanghai Subject Chief
Scientist (14XD1401600).

## References

  * Ai et al., (2021) Ai, C., Linton, O., Motegi, K., and Zhang, Z. (2021).  A unified framework for efficient estimation of general treatment models.  Quantitative Economics, 12(3):779–816.
  * Andrews, (1994) Andrews, D. (1994).  Chapter 37 empirical process methods in econometrics.  Handbook of Econometrics, 4(05):2247–2294.
  * Bauer and Bauer, (1994) Bauer, P. and Bauer, M. M. (1994).  Testing equivalence simultaneously for location and dispersion of two normally distributed populations.  Biometrical journal, 36(6):643–660.
  * Chan et al., (2016) Chan, K. C. G., Yam, S. C. P., and Zhang, Z. (2016).  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.
  * Dasgupta et al., (2015) Dasgupta, T., Pillai, N. S., and Rubin, D. B. (2015).  Causal inference from 2k factorial designs by using potential outcomes.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 77(4):727–753.
  * Dewey et al., (2002) Dewey, D., Kaplan, B. J., Crawford, S. G., and Wilson, B. N. (2002).  Developmental coordination disorder: associated problems in attention, learning, and psychosocial adjustment.  Human movement science, 21(5-6):905–918.
  * Dong et al., (2021) Dong, Y., Lee, Y.-Y., and Gou, M. (2021).  Regression discontinuity designs with a continuous treatment.  Journal of the American Statistical Association, (just-accepted):1–31.
  * Elizabeth et al., (2003) Elizabeth, Johnson, , , Francesca, Dominici, , , Michael, Griswold, , , and Scott (2003).  Disease cases and their medical costs attributable to smoking: an analysis of the national medical expenditure survey.  Journal of Econometrics.
  * Farrington and Andrews, (2003) Farrington, P. and Andrews, N. (2003).  Application to infectious.  Monitoring the Health of Populations: Statistical Principles and Methods for Public Health Surveillance, page 203.
  * Fong et al., (2018) Fong, C., Hazlett, C., and Imai, K. (2018).  Covariate balancing propensity score for a continuous treatment: Application to the efficacy of political advertisements.  The Annals of Applied Statistics, 12(1):156–177.
  * Hainmueller, (2012) Hainmueller, J. (2012).  Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies.  Political analysis, 20(1):25–46.
  * Hainmueller et al., (2014) Hainmueller, J., Hopkins, D. J., and Yamamoto, T. (2014).  Causal inference in conjoint analysis: Understanding multidimensional choices via stated preference experiments.  Social Science Electronic Publishing.
  * Hirano and Imbens, (2004) Hirano, K. and Imbens, G. W. (2004).  The propensity score with continuous treatments.  Applied Bayesian modeling and causal inference from incomplete-data perspectives, 226164:73–84.
  * Hsu et al., (2020) Hsu, Y.-C., Lai, T.-C., and Lieli, R. P. (2020).  Counterfactual treatment effects: Estimation and inference.  Journal of Business & Economic Statistics, pages 1–16.
  * Huber et al., (2020) Huber, M., Hsu, Y.-C., Lee, Y.-Y., and Lettry, L. (2020).  Direct and indirect effects of continuous treatments based on generalized propensity score weighting.  Journal of Applied Econometrics, 35(7):814–840.
  * Imai and Dyk, (2004) Imai, K. and Dyk, D. V. (2004).  Causal inference with general treatment regimes.  Publications of the American Statistical Association, 99(467):854–866.
  * Imai and Ratkovic, (2014) Imai, K. and Ratkovic, M. (2014).  Covariate balancing propensity score.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 76(1):243–263.
  * Imai and Van Dyk, (2004) Imai, K. and Van Dyk, D. A. (2004).  Causal inference with general treatment regimes: Generalizing the propensity score.  Journal of the American Statistical Association, 99(467):854–866.
  * Imbens, (2000) Imbens, G. W. (2000).  The role of the propensity score in estimating dose-response functions.  Biometrika, 87(3):706–710.
  * Kanioglou et al., (2005) Kanioglou, A., Tsorbatzoudis, H., and Barkoukis, V. (2005).  Socialization and behavioral problems of elementary school pupils with developmental coordination disorder.  Perceptual and Motor Skills, 101(1):163–173.
  * Larsen, (1999) Larsen, M. D. (1999).  An analysis of survey data on smoking using propensity scores.  Sankhyā: The Indian Journal of Statistics, Series B, pages 91–105.
  * Miller et al., (2021) Miller, H. L., Sherrod, G. M., Mauk, J. E., Fears, N. E., Hynan, L. S., and Tamplain, P. M. (2021).  Shared features or co-occurrence? evaluating symptoms of developmental coordination disorder in children and adolescents with autism spectrum disorder.  Journal of Autism and Developmental Disorders, pages 1–13.
  * Newey and McFadden, (1994) Newey, W. K. and McFadden, D. (1994).  Large sample estimation and hypothesis testing.  Handbook of econometrics, 4:2111–2245.
  * Rasmussen and Gillberg, (2000) Rasmussen, P. and Gillberg, C. (2000).  Natural outcome of adhd with developmental coordination disorder at age 22 years: a controlled, longitudinal, community-based study.  Journal of the American Academy of Child & Adolescent Psychiatry, 39(11):1424–1431.
  * Robbins et al., (2020) Robbins, M. W., Griffin, B. A., Shih, R. A., and Slaughter, M. E. (2020).  Robust estimation of the causal effect of time-varying neighborhood factors on health outcomes.  Statistics in medicine, 39(5):544–561.
  * Robins et al., (2000) Robins, J. M., Hernan, M. A., and Brumback, B. (2000).  Marginal structural models and causal inference in epidemiology.
  * Rosenbaum and Rubin, (1983) Rosenbaum, P. R. and Rubin, D. B. (1983).  The central role of the propensity score in observational studies for causal effects.  Biometrika, 70(1):41–55.
  * Rosenbaum and Rubin, (1984) Rosenbaum, P. R. and Rubin, D. B. (1984).  Reducing bias in observational studies using subclassification on the propensity score.  Journal of the American statistical Association, 79(387):516–524.
  * Rosenbaum and Rubin, (1985) Rosenbaum, P. R. and Rubin, D. B. (1985).  Constructing a control group using multivariate matched sampling methods that incorporate the propensity score.  The American Statistician, 39(1):33–38.
  * Rubin, (2000) Rubin, D. B. (2000).  Statistical issues in the estimation of the causal effects of smoking due to the conduct of the tobacco industry.  In Statistical science in the courtroom, pages 321–351. Springer.
  * Schumaker, (2007) Schumaker, L. (2007).  Spline functions: basic theory.  Cambridge University Press.
  * Spitzer et al., (1980) Spitzer, R. L., Md, K. K., and Williams, J. B. (1980).  Diagnostic and statistical manual of mental disorders.  In American psychiatric association. Citeseer.
  * Van der Vaart, (2000) Van der Vaart, A. W. (2000).  Asymptotic statistics, volume 3.  Cambridge university press.
  * Vegetabile et al., (2020) Vegetabile, B. G., Gillen, D. L., and Stern, H. S. (2020).  Optimally balanced gaussian process propensity scores for estimating treatment effects.  Journal of the Royal Statistical Society: Series A (Statistics in Society), 183(1):355–377.
  * Viet, (2008) Viet, C. N. (2008).  A note on estimation of the average treatment effect and average partial effect in nonlinear models.  Economics Bulletin, 15.
  * Williams and Crespi, (2020) Williams, J. R. and Crespi, C. M. (2020).  Causal inference for multiple continuous exposures via the multivariate generalized propensity score.
  * Xiong et al., (2017) Xiong, C., Yu, M., and Shao, J. (2017).  Treatment recommendation and parameter estimation under single-index contrast function.  Statistical Theory and Related Fields, 1(2):171–181.
  * Yiu and Su, (2018) Yiu, S. and Su, L. (2018).  Covariate association eliminating weights: a unified weighting framework for causal effect estimation.  Biometrika, 105(3):709–722.
  * Zeger et al., (2000) Zeger, S. L., Wyant, T., Miller, L. S., and Samet, J. (2000).  Statistical Testimony on Damages in Minnesota v. Tobacco Industry.  Statistical Science in the Courtroom.
  * Zhao and Percival, (2017) Zhao, Q. and Percival, D. (2017).  Entropy balancing is doubly robust.  Journal of Causal Inference, 5(1).
  * Zhu et al., (2015) Zhu, Y., Coffman, D. L., and Ghosh, D. (2015).  A boosting algorithm for estimating generalized propensity scores with continuous treatments.  Journal of causal inference, 3(1):25–40.
  * Zubizarreta, (2015) Zubizarreta, J. R. (2015).  Stable weights that balance covariates for estimation with incomplete outcome data.  Journal of the American Statistical Association, 110(511):910–922.

## APPENDIX

### A.1. Proof of Proposition 1

Using the law of total expectation and Assumption 1, we can deduce that

$\displaystyle\begin{split}\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}]&=E[\frac{f(\textbf{T})}{f(\textbf{T}\mid\textbf{X})}(Y-s(\textbf{T};\theta))^{2}]\\\
&=\mathbb{E}(\mathbb{E}[\frac{f(\textbf{T})}{f(\textbf{T}\mid\textbf{X})}(Y-s(\textbf{T};\theta))^{2}]\mid\textbf{T}=\textbf{t},\textbf{X}=\textbf{x})\\\
&=\mathbb{E}(\frac{f(\textbf{t})}{f(\textbf{t}\mid\textbf{x})}\mathbb{E}([(Y-s(\textbf{T};\theta))^{2}]\mid\textbf{T}=\textbf{t},\textbf{X}=\textbf{x}))\\\
&=\int_{\mathcal{T}\times\mathcal{X}}\frac{f(\textbf{t})}{f(\textbf{t}\mid\textbf{x})}\mathbb{E}[(Y(\textbf{T})-s(\textbf{T};\theta))^{2}\mid\textbf{T}=\textbf{t},\textbf{X}=\textbf{x}]f(\textbf{t}\mid\textbf{x})d\textbf{t}d\textbf{x}\\\
&=\int_{\mathcal{T}\times\mathcal{X}}\mathbb{E}[(Y(\textbf{T})-s(\textbf{T};\theta))^{2}\mid\textbf{T}=\textbf{t},\textbf{X}=\textbf{x}]f(\textbf{t})f(\textbf{x})d\textbf{t}d\textbf{x}\\\
&=\int_{\mathcal{T}\times\mathcal{X}}\mathbb{E}[(Y(\textbf{t})-s(\textbf{t};\theta))^{2}\mid\textbf{X}=\textbf{x}]f(\textbf{t})f(\textbf{x})d\textbf{t}d\textbf{x}\\\
&=\int_{\mathcal{T}\times\mathcal{X}}\mathbb{E}[(Y(\textbf{T})-s(\textbf{T};\theta))^{2}f(\textbf{t})d\textbf{t}\quad(\text{using
Assumption 1}).\end{split}$

Hence, we complete the proof of Proposition 1.

### A.2. Derivation of optimization procedure

Considering the optimization problem defined in Section 3.2:

$\displaystyle\text{min}_{w}\sum_{i=1}^{n}w_{i}log(\frac{w_{i}}{v_{i}})\qquad$

s.t.

$\displaystyle\sum_{i=1}^{n}w_{i}g(\textbf{T}_{i},\textbf{X}_{i})=\mathbf{0},\
\sum_{i=1}^{n}w_{i}=1,\ w_{i}>0\ \forall i.$ (23)

Using the standard Lagrange multiplier method to solve this optimization
problem, we construct the Lagrangian as

$\mathcal{L}(w_{i},\alpha,\beta)=\sum_{i=1}^{n}w_{i}log(\frac{w_{i}}{v_{i}})+(\alpha-1)(\sum_{i=1}^{n}w_{i}-1)+\gamma^{{}^{\prime}}\sum_{i=1}^{n}w_{i}g(\textbf{T}_{i},\textbf{X}_{i}),$
(24)

where $\alpha$ and $\beta$ are Langrage multipliers. Taking derivatives with
respect to $w_{i}$, the first order condition is

$\displaystyle\begin{split}\frac{\partial}{\partial
w_{i}}\mathcal{L}(w_{i},\alpha,\beta)&=(log(\frac{w_{i}}{v_{i}})+1)+\alpha-1+\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})\\\
&=log(\frac{w_{i}}{v_{i}})+\alpha+\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})\\\
&=0.\end{split}$

Then we have

$\displaystyle w_{i}\propto
v_{i}exp(-[\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})+\alpha]).$

Enforcing the constraint that the sum of weights must be one, it follows that

$\displaystyle\begin{split}w_{i}&=\frac{v_{i}exp(-[\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})+\alpha])}{\sum_{i=1}^{n}v_{i}exp(-[\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})+\alpha])}\\\
&=\frac{v_{i}exp(-[\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})])}{\sum_{i=1}^{n}v_{i}exp(-[\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})])}.\end{split}$
(25)

Next, for notational simplicity, we redefine the weights as a parameter
$\lambda_{i}$ divided by a normalizing constant
$\lambda=\sum_{i=1}^{n}\lambda_{i}$, that is,

$\displaystyle
w_{i}=\frac{v_{i}exp(-[\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})])}{\sum_{i=1}^{n}v_{i}exp(-[\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})])}\equiv\frac{\lambda_{i}}{\lambda}.$

Taking the log transformation, we obtain that

$\displaystyle
log(\lambda_{i})=log(v_{i})-\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i}),$

implying

$\displaystyle\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i})=log(v_{i})-log(\lambda_{i}).$

Using these relations, we now insert the weights into the right-hand side of
equation (24) and obtain

$\displaystyle\begin{split}&\sum_{i=1}^{n}w_{i}log(\frac{w_{i}}{v_{i}})+\gamma^{{}^{\prime}}\sum_{i=1}^{n}w_{i}g(\textbf{T}_{i},\textbf{X}_{i})+(\alpha-1)(\sum_{i=1}^{n}w_{i}-1)\\\
&=\frac{1}{\lambda}(\sum_{i=1}^{n}\lambda_{i}[log(\lambda_{i})-log(v_{i})+log(\lambda)]+\sum_{i=1}^{n}\lambda_{i}[log(v_{i})-log(\lambda_{i})])\\\
&=\frac{log(\lambda)}{\lambda}(\sum_{i=1}^{n}\lambda_{i})\\\
&=log(\lambda)=log(\sum_{i=1}^{n}v_{i}exp(-\gamma^{{}^{\prime}}g(\textbf{T}_{i},\textbf{X}_{i}))).\end{split}$
(26)

This final line of Equation (26) is no longer a function of the weights. It
can therefore be solved by convex algorithms. Finally, because we set $v_{i}$
to $\frac{1}{n}$ throughout, we can obtain the objective function presented in
Section 3.2.

### A.3.2 Large Sample Properties

#### A.3.1 Proof of Theorem 1

We first show that the conclusion of Theorem(i).

Since $\hat{\theta}$ (as a estimator of $\theta^{*}$) is a unique minimizer of
$\frac{1}{n}\sum_{i=1}^{n}\hat{w_{i}}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}$(regarding
$\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}]$, according to the theory of
M-estimation (Van der Vaart, (2000), Theorem 5.7), if

$\displaystyle\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}\hat{w_{i}}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}-\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}])\mid\to_{p}0,$

then $\hat{\theta}\to_{p}\theta^{*}$. Note that

$\displaystyle\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}\hat{w_{i}}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}-\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}])\mid$
$\displaystyle\leq\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}(\hat{w_{i}}-w_{i})(Y_{i}-s(\textbf{T}_{i};\theta))^{2}\mid$
$\displaystyle+\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}w_{i}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}-\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}])\mid.$
(27)

We first show that
$\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}(\hat{w_{i}}-w_{i})(Y_{i}-s(\textbf{T}_{i};\theta))^{2}\mid$
is $o_{p}(1)$. Using the Causchy-Schwarz inequality and the fact that
$\hat{w}\to^{L^{2}}w$ (Ai, et al.(2020), Theorem 3), we have

$\displaystyle\begin{split}\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}(\hat{w_{i}}-w_{i})(Y_{i}-s(\textbf{T}_{i};\theta))^{2}\mid&\leq\\{\frac{1}{n}\sum_{i=1}^{n}(\hat{w_{i}}-w_{i})^{2}\\}^{1/2}\text{sup}_{\theta\in\Theta}\\{\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}\\}^{1/2}\\\
&\leq
o_{p}(1)\\{\text{sup}_{\theta\in\Theta}\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}]+o_{p}(1)\\}^{1/2}\\\
&=o_{p}(1).\end{split}$

Thereafter, under Assumption 6, we can conclude that
$\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}w_{i}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}-\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}])\mid$
is also $o_{p}(1)$ (Newey and McFadden, (1994), Lemma 2.4). Hence, we complete
the proof for Theorem 1(i). Next, we give the proof of Theorem 1(ii).
Define

$\hat{\theta}^{*}=\text{argmin}_{\theta}\sum_{i=1}^{n}w_{i}(Y_{i}-s(\mathbf{T}_{i};\theta))^{2}.$

Assume that
$\frac{1}{n}\sum_{i=1}^{n}w_{i}(Y_{i}-s(\mathbf{T}_{i};\hat{\theta}^{*}))h(\mathbf{T}_{i};\hat{\theta}^{*}))=o_{p}(n^{-1/2})$
holds with probablility to one as $n\to\infty$

By Assumption 5 and the uniform law of large number, one can get that

$\frac{1}{n}\sum_{i=1}^{n}w_{i}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}\to\mathbb{E}\\{w(Y-s(\textbf{T};\theta))^{2}\\}\
\text{in probability uniformly over }\ \theta,$

which implies $\mid\mid\hat{\theta}^{*}-\theta^{*}\mid\mid\to_{p}0$. Let

$r(\theta)=2\mathbb{E}\\{w(Y-s(\textbf{T};\theta))h(\textbf{T};\theta)\\},$

which ia a differentiable function in $\theta$ and $r(\theta^{*})=0$. By mean
value theorem, we have

$\sqrt{n}r(\hat{\theta}^{*})-\bigtriangledown_{\theta}r(\zeta)\cdot\sqrt{(}n)(\hat{\theta}^{*}-\theta^{*})=\sqrt{n}r(\theta^{*})=0$

where $\zeta$ lies on the line joining $\hat{\theta}^{*}$ and $\theta^{*}$.
Since $\bigtriangledown_{\theta}r(\theta)$ is continuous at $\theta^{*}$ and
$\mid\mid\hat{\theta}^{*}-\theta^{*}\mid\mid\to_{p}0$, then

$\sqrt{n}(\hat{\theta}^{*}-\theta^{*})=\bigtriangledown_{\theta}r(\theta^{*})^{-1}\cdot\sqrt{n}r(\hat{\theta}^{*})+o_{p}(1)$

Define the empirical process

$G_{n}(\theta)=\frac{2}{\sqrt{n}}\sum_{i=1}^{n}\\{w_{i}(Y_{i}-s(\textbf{T}_{i};\theta))h(\textbf{T}_{i};\theta)-\mathbb{E}\\{w(Y-s(\textbf{T};\theta))h(\textbf{T};\theta)\\}\\}.$

Then we have

$\begin{split}&\sqrt{n}(\hat{\theta}^{*}-\theta^{*})\\\
&=\bigtriangledown_{\theta}r(\theta^{*})^{-1}\cdot\\{\sqrt{n}r(\hat{\theta}^{*})-\frac{2}{\sqrt{n}}\sum_{i=1}^{n}\\{w_{i}(Y_{i}-s(\textbf{T}_{i};\hat{\theta}^{*}))h(\textbf{T}_{i};\hat{\theta}^{*})+\frac{2}{\sqrt{n}}\sum_{i=1}^{n}\\{w_{i}(Y_{i}-s(\textbf{T}_{i};\hat{\theta}^{*}))h(\textbf{T}_{i};\hat{\theta}^{*})\\}\\\
&=-\bigtriangledown_{\theta}r(\theta^{*})^{-1}\cdot
G_{n}(\hat{\theta}^{*})+o_{p}(1)\\\
&=U^{-1}\cdot\\{G_{n}(\hat{\theta}^{*})-G_{n}(\theta^{*})+G_{n}(\theta^{*})\\}+o_{p}(1).\end{split}$

By Assumption 5, 6, Theorem 4 and 5 of Andrews(1994), we have
$G_{n}(\hat{\theta}^{*})-G_{n}(\theta^{*})\to_{p}o$. Thus,

$\sqrt{n}(\hat{\theta}^{*}-\theta^{*})=U^{-1}\frac{2}{\sqrt{n}}\sum_{i=1}^{n}\\{w_{i}(Y_{i}-s(\textbf{T}_{i};\theta^{*}))h(\textbf{T}_{i};\theta^{*})\\}+o_{p}(1),$

then we can get that the asymptotic variance of
$\sqrt{n}(\hat{\theta}^{*}-\theta^{*})$ is $V$. Therefore,
$\sqrt{n}(\hat{\theta}^{*}-\theta^{*})\to_{d}N(0,V)$. Next, we will prove
$\hat{\theta}\to_{p}\hat{\theta}^{*}$. Since

$\displaystyle\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}\hat{w_{i}}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}-\frac{1}{n}\sum_{i=1}^{n}w_{i}(Y_{i}-s(\textbf{T}_{i};\theta))^{2})\mid$
$\displaystyle\leq\text{sup}_{\theta\in\Theta}\mid\frac{1}{n}\sum_{i=1}^{n}(\hat{w_{i}}-w_{i})(Y_{i}-s(\textbf{T}_{i};\theta))^{2}\mid$
$\displaystyle\leq\\{\frac{1}{n}\sum_{i=1}^{n}(\hat{w_{i}}-w_{i})^{2}\\}^{1/2}\text{sup}_{\theta\in\Theta}\\{\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-s(\textbf{T}_{i};\theta))^{2}\\}^{1/2}$
$\displaystyle\leq
o_{p}(1)\\{\text{sup}_{\theta\in\Theta}\mathbb{E}[w(Y-s(\textbf{T};\theta))^{2}]+o_{p}(1)\\}^{1/2}$
$\displaystyle=o_{p}(1),$

which implies $\hat{\theta}^{*}\to_{p}\hat{\theta}$. Then by Slutskey’s
Theorem, we can draw the conclusion that
$\sqrt{n}(\hat{\theta}-\theta^{*})\to_{d}N(0,V)$. Therefore, we have completed
the proof of Theorem 1.

#### A.2 Proof of Theorem 2

To obtain the convergence rate of the estimate
$\hat{s}(\textbf{t})=B(\textbf{t})^{{}^{\prime}}\hat{\beta}$, we first obtain
the convergence rate of $\hat{\beta}$.

Note that

$\displaystyle\begin{split}\hat{\beta}-\beta^{*}&=(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}}\hat{W}Y-\beta^{*}\\\
&=(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}}(\hat{W}-W)Y+(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}}(WY-E(WY\mid\textbf{T}))\\\
&+(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}}(E(WY\mid\textbf{T})-Z_{n}\beta^{*})\\\
&\equiv B_{1}+B_{2}+B_{3}.\end{split}$ (28)

First, we compute the convergence order of $B_{1}$. Let $H=(\hat{W}-W)Y$ and
$\hat{\Lambda}=\frac{1}{n}Z_{n}^{{}^{\prime}}Z_{n}$; then,

$\displaystyle\begin{split}\mid\mid
B_{1}\mid\mid^{2}&=\mid\mid(n\hat{\Lambda})^{-1}Z_{n}^{{}^{\prime}}H\mid\mid^{2}\\\
&=n^{-2}tr(H^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-2}Z_{n}^{{}^{\prime}}H)\\\
&=n^{-2}tr((\hat{\Lambda})^{-1}Z_{n}^{{}^{\prime}}HH^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1})\\\
&=n^{-2}tr((\hat{\Lambda})^{-1/2}Z_{n}^{{}^{\prime}}HH^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1/2}(\hat{\Lambda})^{-1})\\\
&\leq\lambda_{max}((\hat{\Lambda})^{-1})n^{-2}tr((\hat{\Lambda})^{-1/2}Z_{n}^{{}^{\prime}}HH^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1/2})\\\
&=\lambda_{max}((\hat{\Lambda})^{-1})n^{-2}tr(HH^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1}Z_{n}^{{}^{\prime}})\\\
&=\lambda_{max}((\hat{\Lambda})^{-1})n^{-1}tr(HH^{{}^{\prime}}Z_{n}(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}})\\\
&\leq\lambda_{max}((\hat{\Lambda})^{-1})n^{-1}\mid\mid H\mid\mid^{2}\\\
&=\lambda_{max}((\hat{\Lambda})^{-1})n^{-1}Y^{{}^{\prime}}(\hat{W}-W)^{{}^{\prime}}(\hat{W}-W)Y\\\
&=\lambda_{max}((\hat{\Lambda})^{-1})sup_{(\textbf{t},\textbf{x})}\mid\hat{W}-W\mid^{2}\frac{1}{n}Y^{{}^{\prime}}Y\\\
&\leq O_{p}(1)O_{p}(\frac{1}{n})O_{p}(1)\\\ &=O_{p}(\frac{1}{n}),\end{split}$
(29)

where the first inequality follows from $tr(AB)\leq\lambda_{max}(B)tr(A)$ for
any symmetric matrix $B$ and positive semidefinite matrix $A$, the second
inequality follows from the fact that
$Z_{n}(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}}$ is a projection
matrix with a maximum eigenvalue of 1, and the last inequality follows from
the fact that $\mid\lambda_{max}((\hat{\Lambda})^{-1})\mid=O_{p}(1)$, Ai et
al., 2020, Lemma 3.1, Corollary 3.3, and
$\frac{1}{n}Y^{{}^{\prime}}Y=O_{p}(1)$.

Thereafter, we compute the convergence order of $B_{2}$. Let
$E=WY-E(WY\mid\textbf{T})$, then

$\displaystyle\begin{split}\mid\mid
B_{2}\mid\mid^{2}&=\mid\mid(n\hat{\Lambda})^{-1}Z_{n}^{{}^{\prime}}E\mid\mid^{2}\\\
&=n^{-2}tr(E^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-2}Z_{n}^{{}^{\prime}}E)\\\
&=n^{-2}tr((\hat{\Lambda})^{-1}Z_{n}^{{}^{\prime}}EE^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1})\\\
&=n^{-2}tr(Z_{n}^{{}^{\prime}}EE^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-2})\\\
&\leq\lambda_{max}((\hat{\Lambda})^{-2})n^{-2}tr(Z_{n}^{{}^{\prime}}EE^{{}^{\prime}}Z_{n})\\\
&=\lambda_{max}((\hat{\Lambda})^{-2})n^{-2}\mid\mid
E^{{}^{\prime}}Z_{n}\mid\mid^{2}\\\ &=O_{p}(\frac{Q}{n}),\end{split}$ (30)

where the last equality follows the fact that
$\mid\lambda_{max}((\hat{\Lambda})^{-1})\mid=O_{p}(1)$ and $n^{-2}\mid\mid
E^{{}^{\prime}}Z_{n}\mid\mid^{2}=O_{p}(\frac{Q}{n})$ by Markov’s inequality.

Finally, we compute the convergence order of $B_{3}$. Let
$M=E(WY\mid\textbf{T})-Z_{n}\beta^{*}$, then

$\displaystyle\begin{split}\mid\mid
B_{3}\mid\mid^{2}&=\mid\mid(n\hat{\Lambda})^{-1}Z_{n}^{{}^{\prime}}M\mid\mid^{2}\\\
&=n^{-2}tr(M^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-2}Z_{n}^{{}^{\prime}}M)\\\
&=n^{-2}tr((\hat{\Lambda})^{-1}Z_{n}^{{}^{\prime}}MM^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1})\\\
&=n^{-2}tr((\hat{\Lambda})^{-1/2}Z_{n}^{{}^{\prime}}MM^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1/2}(\hat{\Lambda})^{-1})\\\
&\leq\lambda_{max}((\hat{\Lambda})^{-1})n^{-2}tr((\hat{\Lambda})^{-1/2}Z_{n}^{{}^{\prime}}MM^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1/2})\\\
&=\lambda_{max}((\hat{\Lambda})^{-1})n^{-2}tr(Z_{n}^{{}^{\prime}}MM^{{}^{\prime}}Z_{n}(\hat{\Lambda})^{-1})\\\
&=\lambda_{max}((\hat{\Lambda})^{-1})n^{-1}tr(MM^{{}^{\prime}}Z_{n}(Z_{n}^{{}^{\prime}}Z_{n})^{-1}Z_{n}^{{}^{\prime}})\\\
&=\lambda_{max}((\hat{\Lambda})^{-1})n^{-1}\sum_{i=1}^{n}\\{E(w_{i}Y_{i}\mid\textbf{T}_{i})-B(\textbf{T}_{i}^{{}^{\prime}})\beta^{*}\\}^{2}.\end{split}$
(31)

Since

$\displaystyle
n^{-1}\sum_{i=1}^{n}\\{E(w_{i}Y_{i}\mid\textbf{T}_{i})-B(\textbf{T}_{i}^{{}^{\prime}})\beta^{*}\\}^{2}\leq\text{sup}_{\textbf{t}\in\mathcal{T}}\mid
E(wY\mid\textbf{T}=\textbf{t})-B(\textbf{t})^{{}^{\prime}}\beta^{*}\mid^{2}=O_{p}(Q^{-2\tilde{\alpha}}),$

we can deduce that $\mid\mid B_{3}\mid\mid^{2}\leq
O_{p}(Q^{-2\tilde{\alpha}})$.

By combining Equations (29), (30) and (31), we have
$\mid\mid\hat{\beta}-\beta^{*}\mid\mid=O_{p}(\sqrt{\frac{Q}{n}}+Q^{-\alpha})$.

Next, the proof of the convergence of $\hat{s}(\mathbf{t})$ is provided. Let
$\Lambda=E[B(\textbf{t})^{{}^{\prime}}B(\textbf{t})]$, then we have

$\displaystyle\begin{split}\int_{\mathcal{T}}\mid\hat{s}(\textbf{t})-s(\textbf{t})\mid^{2}dF_{\textbf{T}}(\textbf{t})&=\int_{\mathcal{T}}\mid
B(\textbf{t})^{{}^{\prime}}\hat{\beta}-B(\textbf{t})^{{}^{\prime}}\beta^{*}+B(\textbf{t})^{{}^{\prime}}\beta^{*}-s(\textbf{t})\mid^{2}dF_{\textbf{T}}(\textbf{t})\\\
&\leq
2(\hat{\beta}-\beta^{*})^{{}^{\prime}}\int_{\mathcal{T}}[B(\textbf{t})^{{}^{\prime}}B(\textbf{t})]dF_{\textbf{T}}(\textbf{t})(\hat{\beta}-\beta^{*})+2\int_{\mathcal{T}}\mid
B(\textbf{t})^{{}^{\prime}}\beta^{*}-s(\textbf{t})\mid^{2}dF_{\textbf{T}}(\textbf{t})\\\
&\leq
2\mid\mid\hat{\beta}-\beta^{*}\mid\mid^{2}\lambda_{max}(\Lambda)+2\int_{\mathcal{T}}\mid
B(\textbf{t})^{{}^{\prime}}\beta^{*}-s(\textbf{t})\mid^{2}dF_{\textbf{T}}(\textbf{t})\\\
&=O_{p}(\frac{Q}{n}+Q^{-2\tilde{\alpha}})\end{split}$ (32)

and

$\displaystyle\begin{split}\text{sup}_{\textbf{t}\in\mathcal{T}}\mid\hat{s}(\textbf{t})-s(\textbf{t})\mid&=\text{sup}_{\textbf{t}\in\mathcal{T}}\mid
B(\textbf{t})^{{}^{\prime}}\hat{\beta}-B(\textbf{t})^{{}^{\prime}}\beta^{*}+B(\textbf{t})^{{}^{\prime}}\beta^{*}-s(\textbf{t})\mid\\\
&\leq\text{sup}_{\textbf{t}\in\mathcal{T}}\mid\mid
B(\textbf{t}^{{}^{\prime}})\mid\mid\cdot\mid\mid\hat{\beta}-\beta^{*}\mid\mid+\text{sup}_{\textbf{t}\in\mathcal{T}}\mid
B(\textbf{t})^{{}^{\prime}}\beta^{*}-s(\textbf{t})\mid\\\ &\leq
O_{p}(\sqrt{\frac{Q}{n}}+Q^{-\tilde{\alpha}}\\})+O_{p}(Q^{-\tilde{\alpha}})\\\
&=O_{p}(\sqrt{\frac{Q}{n}}+Q^{-\tilde{\alpha}}\\}),\end{split}$ (33)

Therefore, the proof of Theorem 2 is complete.