© Wuhan University 2024

0 Introduction

Many empirical issues in social sciences depend on evaluating treatment or policy effects. The estimation of causal effects is a hot issue in statistics and econometrics. Half of the 2021 Nobel Prize in Economics was awarded to Angrist and Imbens to recognize their contributions to the causal inference methodology^[1,2]. Because the average treatment effect is commonly used to study causal effects in Ref. [3], we consider the average treatment effect to study causal effects in this article.

Randomized Controlled Trails (RCTs) have been considered as the "golden" method to evaluate the average treatment effect, where the treatment assignment is randomized. Hence, the average treatment effect can be obtained by comparing the difference in responses between treatment and control groups. However, it is often infeasible to conduct RCTs due to the ethics, morality, technology, or cost in practice. For example, it is unethical to randomize pregnant women to smoke for studying its effects on neonatal weight, and we should not randomly expose people to air pollution to evaluate the impact of PM2.5 on health.

With the rapid development of science and technology, observational studies are widely used to study treatment effects in which the subjects are observed and recorded in the natural state. Compared with RCTs, observational studies can provide richer information about markets, humans, and their behavioral characteristics. However, the treatment assignment in observational studies is far from random. Due to the lack of randomization, systematic differences in the distribution of baseline covariates between treated and controlled groups usually exist. If these differences cannot be appropriately adjusted, the estimator of treatment effects will be seriously biased. The propensity score (PS), proposed by Rosenbaum and Rubin^[4], has been widely used to adjust for previously mentioned differences, which is defined as the conditional probability of receiving a particular treatment given the baseline covariates. Subjects with identical PS will have the same distribution of the baseline covariates regardless of which groups they come from.

There are usually four propensity score-based methods to adjust for heterogeneity in observational studies, including propensity score matching, propensity score stratification, propensity score inverse probability weighting, and propensity score covariate adjustment^[5-9]. Propensity score covariate adjustment and partial least squares can estimate treatment effects in survival analysis^[10,11]. Unfortunately, the propensity score is usually unknown and needs to be evaluated in observational studies. In many practical applications, PS is estimated by the Logistic or Probit models. However, the misspecification of the propensity score model will bring serious bias to the estimator of treatment effects. Lee^[12] proposed a simple least squares estimator based on propensity score residuals to reduce the bias of the estimator of treatment effects, in which the propensity score was estimated by the Probit model. Inspired by Lee's method^[12], we will study treatment effects with propensity score residuals and propose a more robust approach based on the generalized single-index model to estimate the propensity score.

The rest of the paper is organized as follows. We introduce the statistical inference procedures in Section 1. In Section 2, we evaluate the finite-sample performance of the proposed method through simulation studies. Section 3 applies the proposed method to evaluate the policy of "Sunshine Running" in the university. Conclusions are shown in Section 4.

1 Estimating Procedures of Average Treatment Effect

We adopted the counterfactual framework to study the average treatment effect^[13,14]. In the counterfactual framework, two potential outcomes exist under the treatment and control, respectively. The foundation problem of causal inference is that only one of two potential outcomes can be observed for the same individual^[15]. The average treatment effect is defined as the expectation of the difference between two potential outcomes.

To simplify the expression, we adopt the following notation: the binary variable D denotes the treatment assignment ($D=\mathrm{1}$ for treatment, $D=\mathrm{0}$ for control) and let X represent a p-vector of the baseline covariates;$Y(\mathrm{1})$,$Y(\mathrm{0})$ are potential outcomes under the active treatment $(D=\mathrm{1})$ and control $(D=\mathrm{0})$, respectively. Hence, the observed outcome is $Y=DY(\mathrm{1})+(\mathrm{1}-D)Y(\mathrm{0})$. The observed data $\{(Y_i,D_i,{\mathit{X}}_i),i=\mathrm{1},...,n\}$ are assumed to be the independent copies of $(Y,D,\mathit{X})$. In this paper, we consider the following average treatment effect (ATE) to study the causal effect:

$\tau =E[Y(\mathrm{1})]-E[Y(\mathrm{0})]$(1)

The conditional average treatment effect (CATE) is:

$\tau (\mathit{X})=E[Y(\mathrm{1})-Y(\mathrm{0})|\mathit{X}]=E[Y\left(\mathrm{1}\right)|\mathit{X}]-E[Y(\mathrm{0})|\mathit{X}]$(2)

According to the properties of the conditional expectation, we can have the following:

$\tau =E[\tau (\mathit{X})]$(3)

1.1 Identifiability of Average Treatment Effect

In order to identify the average treatment effect, we have to make the following assumptions.

Remarks: The above three assumptions are commonly used in the causal inference literature. Assumption 2 means that all the covariates that affect both the outcome and treatment assignment are measured; Assumption 3 ensures overlap exists in the distributions of covariates in treatment and control groups.

1.2 Models

There usually exist two types of models in the study of treatment effects: the outcome model, which links the response with the treatment assignment and baseline covariates, and the treatment assignment model, which defines the conditional probability of receiving a particular treatment given the baseline covariates. We generate the outcome through the following partial linear model:

$Y_i=\alpha D_i+g\left({\mathit{X}}_{\mathit{i}}\right)+{\epsilon }_i,i=\mathrm{1},\cdots ,n$(4)

where the parameter $\alpha$ denotes the average treatment effect, $g(\mathit{X})$ denotes the unknown smoothed function, and $\epsilon$ is the random error with $E[\epsilon |D,\mathit{X}]=\mathrm{0}$. We consider the following generalized single-index model for the propensity score:

$P(D=\mathrm{1}|\mathit{X})=\frac{\mathrm{e}\mathrm{x}\mathrm{p}(h({\mathit{\gamma }}_{\mathrm{0}}^{\mathrm{T}}\mathit{X}))}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(h({\mathit{\gamma }}_{\mathrm{0}}^{\mathrm{T}}\mathit{X}))}$(5)

where $h(\cdot )$ is an unknown smoothed function and ${\mathit{\gamma }}_{\mathrm{0}}$ denotes the unknown parameter. Furthermore, in order to identify model (5), ${\mathit{\gamma }}_{\mathrm{0}}$ needs to satisfy $\Vert {\mathit{\gamma }}_{\mathrm{0}}\Vert =\mathrm{1}$ with ${\mathit{\gamma }}_{\mathrm{0,1}}\ge \mathrm{0}$. The proposed model (5) includes the commonly used Logistic and Probit models as special cases.

1.3 Estimation Method of Average Treatment Effect

Inspired by Lee's method^[12], we use the simple least squares estimator based on propensity score residuals to estimate the parameter $\alpha$. Let $\pi (\mathit{X})$ and $\widehat{\pi }(\mathit{X})$ denote the true and estimated propensity score and parameters ${\widehat{\mathit{\gamma }}}_n$, ${\widehat{\alpha }}_n$ are estimators of $\mathit{\gamma }$, $\alpha$, respectively. Firstly, we estimate the propensity score by the maximum likelihood estimate (MLE). Since the function $h(\cdot )$ in model (5) is unknown, we use the sieve method to approximate the unknown function $h(\cdot )$. Define ${\mathit{H}}_k(\mathit{x})=(h_{\mathrm{0}}(\mathit{x}),h_{\mathrm{1}}(\mathit{x}),$$\cdots ,h_{k-\mathrm{1}}{(\mathit{x}))}^{\mathrm{T}}$to denote the k-dimensional orthogonal basis ${\mathit{C}}_k=(c_{\mathrm{0}},c_{\mathrm{1}},...,c_{k-\mathrm{1}}{)}^{\mathrm{T}}$ denotes k-dimensional coefficients (A detailed description of the construction of orthogonal basis and the related properties can be found in Chen^[16], Dong et al^[17]). For model (5), we can use the following maximum likelihood function:

$L_n(\mathit{\gamma },{\mathit{C}}_k)=\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}\{D_i({\mathit{H}}_k({\mathit{\gamma }}^T{\mathit{X}}_i{)}^T{\mathit{C}}_k)-log(\mathrm{1}+exp({\mathit{H}}_k({\mathit{\gamma }}^T{\mathit{X}}_i{)}^T{\mathit{C}}_k))\}$(6)

Then, $({\widehat{\mathit{\gamma }}}_n,{\widehat{\mathit{C}}}_k)$, the estimator of $(\mathit{\gamma },{\mathit{C}}_k)$, can be obtained by:

$({\widehat{\mathit{\gamma }}}_n,{\widehat{\mathit{C}}}_k)=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{\gamma },{\mathit{C}}_k}{\mathrm{m}\mathrm{i}\mathrm{n}}[-L_n(\mathit{\gamma },{\mathit{C}}_k)+\lambda ({\Vert \mathit{\gamma }\Vert }_{\mathrm{2}}^{\mathrm{2}}-\mathrm{1})]$(7)

where the positive integer value $k$ is a regulation parameter, $\lambda$ is a Lagrange multiplier, and ${\Vert \mathit{\gamma }\Vert }_{\mathrm{2}}^{\mathrm{2}}$ denotes the square of the Euclidean norm of $\mathit{\gamma }$. Once obtained $({\widehat{\mathit{\gamma }}}_n,{\widehat{\mathit{C}}}_k)$, the propensity score can be estimated as:

$\widehat{\pi }(\mathit{X})=\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\widehat{\mathit{C}}}_k)}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\widehat{\mathit{C}}}_k)}$(8)

According to the simple least squares estimator based on the propensity score residuals, we can obtain the estimator of the average treatment effect:

${\widehat{\alpha }}_n=\frac{{\displaystyle \sum _{i=\mathrm{1}}^n}(D_i-\widehat{\pi }({\mathit{X}}_i))(Y_i-{\displaystyle \sum _{j=\mathrm{0}}^q}({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^j)}{{\displaystyle \sum _{i=\mathrm{1}}^n}(D_i-\widehat{\pi }({\mathit{X}}_i{))}^{\mathrm{2}}}$(9)

where q denotes the polynomial order, which is usually taken as 2 in practical applications. Furthermore, by similar arguments as Lee^[12], we establish the asymptotic properties of ${\widehat{\alpha }}_n$, presented in Theorem 1.

To simply notations, we set

$\pi (\mathit{X})=\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^T{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}$(10)

$\mathrm{\Delta }=D-\pi (\mathit{X})$(11)

$V=Y-{\displaystyle \sum _{j=\mathrm{0}}^q}({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^j-\alpha \mathrm{\Delta }$(12)

and its estimations are:

$\widehat{\pi }({\mathit{X}}_i)=\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}{\widehat{\mathit{C}}}_k)}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}{\widehat{\mathit{C}}}_k)}$(13)

${\widehat{\mathrm{\Delta }}}_i=D_i-\widehat{\pi }({\mathit{X}}_i)$(14)

${\widehat{V}}_i=Y_i-{\displaystyle \sum _{j=\mathrm{0}}^q}({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^j-{\widehat{\alpha }}_n{\widehat{\mathrm{\Delta }}}_i$(15)

In order to establish the asymptotic properties of the proposed estimators ${\widehat{\alpha }}_n$, we need the following lemmas.

The moment conditions:

$\begin{array}{l}E[\{Y-{\displaystyle \sum _{j=\mathrm{0}}^q}({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^j-\alpha (D-\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)})\}\\ \cdot \{D-\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}\}]=\mathrm{0}\end{array}$(16)

for ${\widehat{\alpha }}_n$, with $\mathit{\gamma }$ and ${\mathit{C}}_k$ replaced by $\mathit{g}$ and $\mathit{c}$, is

$E[\{Y-{\displaystyle \sum _{j=\mathrm{0}}^q}({\mathit{g}}^{\mathrm{T}}{\mathit{X})}^j-\alpha (D-\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{g}}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}\mathit{c})}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{g}}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}\mathit{c})})\}$

$\cdot \{D-\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{g}}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}\mathit{c})}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{g}}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}\mathit{c})}\}]=\mathrm{0}$(17)

The effect of ${\widehat{\mathit{\gamma }}}_n-\mathit{\gamma }$ and ${\widehat{\mathit{C}}}_k-{\mathit{C}}_k$ on the asymptotic distribution of ${\widehat{\alpha }}_n-\alpha$ would be zero if the derivatives of the left-hand side were zero at $\mathit{g}=\mathit{\gamma }$ and $\mathit{c}={\mathit{C}}_k$. The derivatives with respect to $\mathit{g}$ at $\mathit{g}=\mathit{\gamma }$ and $\mathit{c}={\mathit{C}}_k$ are

$\mathrm{0}-E[V\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}{(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k{))}^{\mathrm{2}}}{\dot{\mathit{H}}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k{\mathit{X}}^{\mathrm{T}}]$(18)

where ${\dot{\mathit{H}}}_k(\cdot )$ denotes the derivative of ${\mathit{H}}_k(\cdot )$. The first term is zero because $E[(D-\pi (\mathit{X}))|\mathit{X}]=\mathrm{0}$, we denote:

$q_{\alpha \gamma }=-E[V\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}{(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k{))}^{\mathrm{2}}}{\dot{\mathit{H}}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k{\mathit{X}}^{\mathrm{T}}]$(19)

and $q_{\alpha \gamma }$ can be estimated consistently by:

${\widehat{q}}_{\alpha \gamma }=-\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{\widehat{V}}_i\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}{\widehat{\mathit{C}}}_k)}{(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}{\widehat{\mathit{C}}}_k{))}^{\mathrm{2}}}{\dot{\mathit{H}}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}{\widehat{\mathit{C}}}_k{\mathit{X}}_i^{\mathrm{T}}$(20)

The derivatives with respect to $\mathit{c}$ at $\mathit{g}=\mathit{\gamma }$ and $\mathit{c}={\mathit{C}}_k$ are

$\mathrm{0}-E[V\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}{(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k{))}^{\mathrm{2}}}{\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}]$(21)

The first term is zero because $E[(D-\pi (\mathit{X}))|\mathit{X}]=\mathrm{0}$, we denote

$q_{\alpha C_k}=-E[V\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k)}{(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}{\mathit{C}}_k{))}^{\mathrm{2}}}{\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X})}^{\mathrm{T}}]$(22)

and $q_{\alpha C_k}$ can be estimated consistently by:

${\widehat{q}}_{\alpha C_k}=-\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}{\widehat{V}}_i\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}{\widehat{\mathit{C}}}_k)}{(\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}{\widehat{\mathit{C}}}_k{))}^{\mathrm{2}}}{\mathit{H}}_k({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}$(23)

We denote the log-likelihood function for model (5) is

$F_n\left(\mathit{\gamma },{\mathit{C}}_k\right)={\displaystyle \sum _{i=\mathrm{1}}^n}\{D_i({\mathit{H}}_k({\mathit{\gamma }}^T{\mathit{X}}_i{)}^T{\mathit{C}}_k)-log(\mathrm{1}+exp({\mathit{H}}_k({\mathit{\gamma }}^T{\mathit{X}}_i{)}^T{\mathit{C}}_k))\}$(24)

then we have the score function with respect to $\mathit{\gamma }$

$\begin{array}{l}{\mathit{s}}_{\gamma i}=\frac{\partial F\left(\mathit{\gamma },{\mathit{C}}_k\right)}{\partial \mathit{\gamma }}\\ =D_i{\dot{\mathit{H}}}_k{\left({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X}}_i\right)}^{\mathrm{T}}{\mathit{C}}_k{\mathit{X}}_i^{\mathrm{T}}-\pi ({\mathit{X}}_i){\dot{\mathit{H}}}_k{\left({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X}}_i\right)}^{\mathrm{T}}{\mathit{C}}_k{\mathit{X}}_i^{\mathrm{T}}\end{array}$(25)

By the definition of the influence function, we have the influence function with respect to $\mathit{\gamma }$,

${\mathit{\eta }}_{\gamma i}=E^{-\mathrm{1}}[{\mathit{s}}_{\gamma }{\mathit{s}}_{\gamma }^{\mathrm{T}}]{\mathit{s}}_{\gamma i}$(26)

then

$\sqrt[]{n}({\widehat{\mathit{\gamma }}}_n-\mathit{\gamma })=\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\mathit{\eta }}_{\gamma i}+o_p(\mathrm{1})$(27)

Similarly, we also have the score function with respect to ${\mathit{C}}_k$,

${\mathit{s}}_{{\mathit{C}}_{k}i}=\frac{\partial F(\mathit{\gamma },{\mathit{C}}_k)}{\partial {\mathit{C}}_k}=D_i{\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}-\pi ({\mathit{X}}_i){\mathit{H}}_k({\mathit{\gamma }}^{\mathrm{T}}{\mathit{X}}_i{)}^{\mathrm{T}}$(28)

By the definition of the influence function, we have the influence function with respect to ${\mathit{C}}_k$,

${\mathit{\eta }}_{C_{k}i}=E^{-\mathrm{1}}[{\mathit{s}}_{C_k}{\mathit{s}}_{C_k}^{\mathrm{T}}]{\mathit{s}}_{C_{k}i}$(29)

then

$\sqrt[]{n}({\widehat{\mathit{C}}}_k-{\mathit{C}}_k)=\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\mathit{\eta }}_{C_{k}i}+o_p(\mathrm{1})$(30)

$\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}\{D_i-\widehat{\pi }({\mathit{X}}_i)\}\{Y_i-{\displaystyle \sum _{j=\mathrm{0}}^q}({\widehat{\mathit{\gamma }}}_n^T{\mathit{X}}_i{)}^j-{\widehat{\alpha }}_n(D_i-\widehat{\pi }({\mathit{X}}_i)\}=\mathrm{0}$(31)

then apply the mean value theorem to ${\widehat{\alpha }}_n$ around $\alpha$ to get

$\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}[\{D_i-\widehat{\pi }({\mathit{X}}_i)\}\{Y_i-{\displaystyle \sum _{j=\mathrm{0}}^q}({\widehat{\mathit{\gamma }}}_n^T{\mathit{X}}_i{)}^j-\alpha (D_i-\widehat{\pi }({\mathit{X}}_i)\}-\{D_i-\widehat{\pi }({\mathit{X}}_i{)\}}^{\mathrm{2}}\{{\widehat{\alpha }}_n-\alpha \}]=\mathrm{0}$(32)

So, the term $\sqrt[]{n}({\widehat{\alpha }}_n-\alpha )$ can be written as

$\begin{array}{l}\sqrt[]{n}({\widehat{\alpha }}_n-\alpha )\\ =\frac{\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}(D_i-\widehat{\pi }({\mathit{X}}_i))(Y_i-{\displaystyle \sum _{j=\mathrm{0}}^q}({\widehat{\mathit{\gamma }}}_n^{\mathrm{T}}{\mathit{X}}_i{)}^j-\alpha (D_i-\widehat{\pi }({\mathit{X}}_i)))}{\frac{\mathrm{1}}{n}{\displaystyle \sum _{i=\mathrm{1}}^n}(D_i-\widehat{\pi }({\mathit{X}}_i{))}^{\mathrm{2}}}\end{array}$(33)

Due to ${\widehat{\mathit{\gamma }}}_n$ and ${\widehat{\mathit{C}}}_k$ , we should go one step further applying the mean value theorem to ${\widehat{\mathit{\gamma }}}_n$ and ${\widehat{\mathit{C}}}_k$ as follows. Expand the numerator of (33) to $({\widehat{\mathit{\gamma }}}_n,{\widehat{\mathit{C}}}_k)$ around $(\mathit{\gamma },{\mathit{C}}_k)$ to get

$\begin{array}{l}\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\widehat{\mathrm{\Delta }}}_i{\widehat{V}}_i=\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\mathrm{\Delta }}_i{V}_i+{\widehat{q}}_{\alpha \gamma }\sqrt[]{n}({\widehat{\mathit{\gamma }}}_n-\mathit{\gamma })\\ +{\widehat{q}}_{\alpha C_k}\sqrt[]{n}({\widehat{\mathit{C}}}_k-{\mathit{C}}_k)\end{array}$(34)

Invoking the uniform LLN (Law of Large Numbers), the denominator of (33), ${\widehat{q}}_{\alpha \gamma }$ and ${\widehat{q}}_{\alpha C_k}$ can be replaced with $E[(D-\pi {(\mathit{X}))}^{\mathrm{2}}]$ ,$q_{\alpha \gamma }$ and $q_{\alpha C_k}$ , respectively. $q_{\alpha \gamma }$ ,${\widehat{q}}_{\alpha \gamma }$ ,$q_{\alpha C_k}$ and ${\widehat{q}}_{\alpha C_k}$ are defined by (19), (20), (22) and (23) in Lemma 1.

Then we have,

$\begin{array}{l}\sqrt[]{n}({\widehat{\alpha }}_n-\alpha )\\ =\frac{\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}{\mathrm{\Delta }}_i{V}_i+q_{\alpha \gamma }\sqrt[]{n}({\widehat{\mathit{\gamma }}}_n-\mathit{\gamma })+q_{\alpha C_k}\sqrt[]{n}({\widehat{\mathit{C}}}_k-{\mathit{C}}_k)}{E[(D-\pi {(\mathit{X}))}^{\mathrm{2}}]}+o_p(\mathrm{1})\end{array}$(35)

By (27) and (30) in Lemma 2, we have

$\sqrt[]{n}({\widehat{\alpha }}_n-\alpha )=\frac{\frac{\mathrm{1}}{\sqrt[]{n}}{\displaystyle \sum _{i=\mathrm{1}}^n}\{{\mathrm{\Delta }}_i{V}_i+q_{\alpha \gamma }{\mathit{\eta }}_{\gamma i}+q_{\alpha C_k}{\mathit{\eta }}_{C_{k}i}\}}{E[(D-\pi {(\mathit{X}))}^{\mathrm{2}}]}+o_p(\mathrm{1})$(36)

Hence,

$\begin{array}{l}\sqrt[]{n}({\widehat{\alpha }}_n-\alpha )\stackrel{d}{\to }N(\mathrm{0},E^{-\mathrm{1}}[(D-\pi {(\mathit{X}))}^{\mathrm{2}}]U{E}^{-\mathrm{1}}\\ \cdot [(D-\pi {(\mathit{X}))}^{\mathrm{2}}])\end{array}$(37)

where

$U=E[\{\mathrm{\Delta }V+q_{\alpha \gamma }{\mathit{\eta }}_{\gamma }+q_{\alpha C_k}{\mathit{\eta }}_{C_k}\}\{\mathrm{\Delta }V+q_{\alpha \gamma }{\mathit{\eta }}_{\gamma }+q_{\alpha C_k}{\mathit{\eta }}_{C_k}{\}}^{\mathrm{T}}]$(38)

The proof is completed.

2 Simulation Studies

In this section, extensive simulation studies are conducted to evaluate the finite-sample performance of the proposed method. We consider different scenarios to simulate the observational studies in reality.

2.1 Simulation 1: Idealized Scenario

The outcome is generated by the following linear regression model:

$Y=\alpha D+{\beta }_{\mathrm{1}}{X}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{X}_{\mathrm{2}}+{\beta }_{\mathrm{3}}{X}_{\mathrm{3}}+\epsilon$(39)

where,$D$denotes the treatment assignment, the covariates $X_{\mathrm{1}}$, $X_{\mathrm{2}}$ and $X_{\mathrm{3}}$ are independent of each other, $X_{\mathrm{1}}$ and $X_{\mathrm{2}}$ follow the standard normal distribution,$X_{\mathrm{3}}$ follows Bernoulli distribution with successful probability of 0.4, the error term $\epsilon$ follows the standard normal distribution, $({\beta }_{\mathrm{1}},{\beta }_{\mathrm{2}},{\beta }_{\mathrm{3}})=(\mathrm{0.5,0},\mathrm{1})$ and $\alpha$ is set to 0 or 0.5. The propensity score is generated by the Logistic model:

$P(D=\mathrm{1}|X_{\mathrm{1}},X_{\mathrm{2}},X_{\mathrm{3}})=\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}})}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}})}$(40)

where $({\gamma }_{\mathrm{1}},{\gamma }_{\mathrm{2}},{\gamma }_{\mathrm{3}})=(\mathrm{0.8,0},-\mathrm{0.6})$, then the proportion of subjects in the treatment group is about 45%. We compare the proposed estimator ${\widehat{\alpha }}_P$ with two estimators:${\widehat{\alpha }}_R$, the estimator based on the covariate adjustment with propensity score^[18,19]; ${\widehat{\alpha }}_L$, and the estimator based on propensity score residuals^[12].

We generate 500 simulated data sets with the total sample size $n$ being 150 or 300. The sample mean, sample standard deviation and sample root mean square error of 500 estimators are given in the columns "Mean", "SD" and "RMSE", respectively. Furthermore, the column "ESD" shows a mean of the estimated standard deviation, and "CP" gives the nominal 95% confidence interval coverage rate using the estimated standard deviation. The bootstrap method is adopted to obtain the estimated variance, and the number of bootstraps is 100. The simulation results are summarized in Table 1.

From Table 1, we obtain that the three estimators ${\widehat{\alpha }}_R$, ${\widehat{\alpha }}_L$ and ${\widehat{\alpha }}_P$ are all approximately unbiased. The corresponding estimated standard error agrees well with the sampling standard error, and the coverage probability of a 95% confidence interval is around the nominal level of the three estimators.

2.2 Simulation 2: Single-Index Model Scenario

The outcome is generated by the following partial single-index model:

$Y=\alpha D+\mathrm{e}\mathrm{x}\mathrm{p}({\beta }_{\mathrm{1}}{X}_{\mathrm{1}}+{\beta }_{\mathrm{2}}{X}_{\mathrm{2}}+{\beta }_{\mathrm{3}}{X}_{\mathrm{3}})+\epsilon$(41)

The treatment assignment model is:

$\begin{array}{l}P(D=\mathrm{1}|X_{\mathrm{1}},X_{\mathrm{2}},X_{\mathrm{3}})\\ =\frac{\mathrm{e}\mathrm{x}\mathrm{p}(({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}}{)}^{\mathrm{2}}+\mathrm{4}({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}}))}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}}{)}^{\mathrm{2}}+\mathrm{4}({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}}))}\end{array}$(42)

The parameter settings are the same as in Simulation 1. The simulation results are summarized in Table 2.

From Table 2, the three estimators ${\widehat{\alpha }}_R$, ${\widehat{\alpha }}_L$ and ${\widehat{\alpha }}_P$ are all approximately unbiased. When the sample size is 150, we can know that the proposed estimator ${\widehat{\alpha }}_P$ has the smallest sample standard deviation and sample root mean square error. The coverage proportions of the proposed estimator ${\widehat{\alpha }}_P$ is apparently closer to the nominal 95% than that of the estimators ${\widehat{\alpha }}_R$ and ${\widehat{\alpha }}_L$. When the sample size is increased to 300, all indicators show that the proposed estimator ${\widehat{\alpha }}_P$ is the most effective.

2.3 Simulation 3: More Complex Scenario

The outcome is generated by the following model:

$Y=\alpha D+({\beta }_{\mathrm{1}}{X}_{\mathrm{1}}+{\beta }_{\mathrm{1}}{X}_{\mathrm{2}}+{\beta }_{\mathrm{1}}{X}_{\mathrm{3}}{)}^{\mathrm{2}}+\epsilon$(43)

The treatment assignment model is:

$P(D=\mathrm{1}|X_{\mathrm{1}},X_{\mathrm{2}},X_{\mathrm{3}})=\frac{\mathrm{e}\mathrm{x}\mathrm{p}({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}}+X_{\mathrm{1}}^{\mathrm{2}})}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}({\gamma }_{\mathrm{1}}{X}_{\mathrm{1}}+{\gamma }_{\mathrm{2}}{X}_{\mathrm{2}}+{\gamma }_{\mathrm{3}}{X}_{\mathrm{3}}+X_{\mathrm{1}}^{\mathrm{2}})}$(44)

The parameter settings are the same as in Simulation 1. The simulation results are summarized in Table 3.

From Table 3, when the treatment assignment model is no longer the single-index model, the three estimators ${\widehat{\alpha }}_R$, ${\widehat{\alpha }}_L$ and ${\widehat{\alpha }}_P$ are all biased. However, the biases of the estimators ${\widehat{\alpha }}_R$ and ${\widehat{\alpha }}_L$ are relatively large, while the bias of the proposed estimator ${\widehat{\alpha }}_P$ is relatively small. The average of the estimated standard deviation of ${\widehat{\alpha }}_P$ is close to the sample standard deviation. The coverage proportions of the proposed estimator ${\widehat{\alpha }}_P$ is much more relative to the nominal level than the estimators ${\widehat{\alpha }}_R$ and ${\widehat{\alpha }}_L$. Therefore, from the simulation results of three different scenarios, the performance of our proposed estimator ${\widehat{\alpha }}_P$ is optimal.

3 Real Data Analysis

Nowadays, due to unhealthy work and rest habits, many college students have the problem of poor physical quality. In order to continuously improve the overall level of physical health of college students, relevant government documents put forward mandatory requirements for the physical health level of college students; many universities also organize "Sunshine Running" to improve the physical health of students^[20,21]. In this section, we apply the proposed method to evaluate the effectiveness of the "Sunshine Running" policy offered by the university to improve college students' physical test scores. The college where we collected the data implemented the "Sunshine Running" policy among the freshmen admitted in 2018 but did not implement the policy in other grades. We collected physical test data from the class of grad 2017 at this university in 2018 as the control group and from the class of grad 2018 in 2019 as the treatment group. In order to make the collected data independent of each other, we randomly selected one male and one female in each major of the university to record their physical test scores, gender, age, and Body Mass Index (BMI, weight divided by height squared). We collected 82 students in the treatment group and 69 students in the control group, for a total of 151 students.

To simplify the expression, we use the following notation: Y denotes the physical test score, D represents the treatment assignment, Sex denotes gender (male=1, female=2), Age denotes age, and BMI denotes Body Mass Index. We consider the following outcome model:

$Y=\alpha D+g(\mathrm{S}\mathrm{e}\mathrm{x},\mathrm{A}\mathrm{g}\mathrm{e},\mathrm{B}\mathrm{M}\mathrm{I})+\epsilon$(45)

where the error term $\epsilon$ follows the standard normal distribution, $g(\cdot ,\cdot ,\cdot )$ is an unknown smoothed function.

We consider the following generalized single-index model to generate the treatment assignment:

$\begin{array}{l}P(D=\mathrm{1}|\mathrm{S}\mathrm{e}\mathrm{x},\mathrm{A}\mathrm{g}\mathrm{e},\mathrm{B}\mathrm{M}\mathrm{I})\\ =\frac{\mathrm{e}\mathrm{x}\mathrm{p}(h({\gamma }_{\mathrm{1}}\mathrm{S}\mathrm{e}\mathrm{x}+{\gamma }_{\mathrm{2}}\mathrm{A}\mathrm{g}\mathrm{e}+{\gamma }_{\mathrm{3}}\mathrm{B}\mathrm{M}\mathrm{I}))}{\mathrm{1}+\mathrm{e}\mathrm{x}\mathrm{p}(h({\gamma }_{\mathrm{1}}\mathrm{S}\mathrm{e}\mathrm{x}+{\gamma }_{\mathrm{2}}\mathrm{A}\mathrm{g}\mathrm{e}+{\gamma }_{\mathrm{3}}\mathrm{B}\mathrm{M}\mathrm{I}))}\end{array}$(46)

where $h(\cdot )$ is an unknown smoothed function, ${\gamma }_{\mathrm{1}},{\gamma }_{\mathrm{2}},{\gamma }_{\mathrm{3}}$ are unknown parameters.

Likewise, we considered three estimators ${\widehat{\alpha }}_R$,${\widehat{\alpha }}_L$ and ${\widehat{\alpha }}_P$. We considered them in the simulation studies to study the effect of the "Sunshine Running" policy. We use the bootstrap to derive the variance of the three estimators, and the number of bootstraps is 100. The estimators, the standard deviations, and p-value of ${\widehat{\alpha }}_R$, ${\widehat{\alpha }}_L$ and ${\widehat{\alpha }}_P$ are summarized in Table 4.

From Table 4, we can conclude that all the estimators confirm that the "Sunshine Running" policy has a significant impact on improving college students' physical test scores. From the standard deviations of the three estimators, the proposed estimator ${\widehat{\alpha }}_P$ is the most effective, and the estimator ${\widehat{\alpha }}_R$ is the worst.

4 Conclusion

Estimating treatment effects or evaluating policy effects is an essential issue in empirical science, which can provide us with the basis for quantitative decision-making. With the development of science and technology, a large number of high-quality observational data have been collected, and how to use high-quality and information-rich observational data to estimate treatment effects is one of the hot issues in current research. The propensity score has been widely used to adjust for heterogeneity in observational studies. However, the misspecification of the propensity score will cause serious bias. In this paper, we proposed the generalized single-index model to estimate the propensity score and used the propensity score residuals to reduce the estimation bias. From the results of simulation studies and real data analysis, our proposed method is more effective than other competitive methods. Next, we will consider more robust methods to estimate the treatment effect, such as deep learning, in the future work.

References

All Tables