Models for the Propensity Score that Contemplate the Positivity Assumption and their Application to Missing Data and Causality

TL;DR

This paper proposes models for propensity scores that explicitly incorporate a lower bound to address the positivity assumption, improving the stability of inverse probability weighting estimators in cases with continuous covariates.

Contribution

It introduces a new modeling approach for propensity scores that ensures a lower bound, reconciling the positivity assumption with generalized linear models.

Findings

Provides a model with explicit lower bounds for propensity scores.

Enhances the stability of IPW estimators when positivity does not hold.

Addresses issues with continuous covariates in propensity score modeling.

Abstract

Generalized linear models are often assumed to fit propensity scores, which are used to compute inverse probability weighted (IPW) estimators. In order to derive the asymptotic properties of IPW estimators, the propensity score is supposed to be bounded away from cero. This condition is known in the literature as strict positivity (or positivity assumption) and, in practice, when it does not hold, IPW estimators are very unstable and have a large variability. Although strict positivity is often assumed, it is not upheld when some of the covariates are continuous. In this work, we attempt to conciliate between the strict positivity condition and the theory of generalized linear models by incorporating an extra parameter, which results in an explicit lower bound for the propensity scores.

Tables

Table 1. Table 2: Estimates for the ATE of smoking in children, their standard errors and confidence intervals.

	Estimator	Standard Error	Confidence Interval
${\hat{\tau }}_O$	-0.1538	0.2188	(-0.5826 , 0.2750)
${\hat{\tau }}_P$	-0.2663	0.1952	(-0.6727 , 0.1402)
${\hat{\tau }}_{LD}$	-0.1911	0.2074	(-0.5737 , 0.1914)
${\hat{\tau }}_{PLD}$	-0.2637	0.2086	(-0.6726 , 0.1452)

Equations

$$P(A=1\mid \mathbf{X},Y)=P(A=1\mid \mathbf{X})=\pi (X).$$

$$\hat{\mu }(\hat{\pi })={\mathbb{P}}_n\left\{\frac{AY}{\hat{\pi }(\text{X})}\right\},$$

$$\pi (\text{X})=\varphi ({\mathbf{\beta }}_0^T\text{X}),$$

$$\pi (\text{X})=(1-{\epsilon }_0)\varphi ({\mathbf{\beta }}_0^T\text{X})+{\epsilon }_0,$$

$$\pi (\mathbf{X},\mathbf{\theta })=(1-\epsilon )\varphi ({\mathbf{\beta }}^T\mathbf{X})+\epsilon .$$

$${\hat{\mu }}_P={\mathbb{P}}_n\left\{\frac{AY}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\},$$

$$\inf \left\{z\in \mathbb{R}:P(|{\mathbf{\beta }}^T\mathbf{X}-z|<r)>0,for\; allr>0\right\}=-\infty .$$

$$(Y^{(0)},Y^{(1)})\coprod T|\mathbf{X}.$$

$$E\left(Y^{(1)}\right)=E\left\{\frac{TY}{\pi (\mathbf{X})}\right\},E\left(Y^{(0)}\right)=E\left\{\frac{(1-T)Y}{1-\pi (\mathbf{X})}\right\},$$

$$\hat{\tau }(\hat{\pi })={\mathbb{P}}_n\left\{\frac{TY}{\hat{\pi }(\text{X})}\right\}-{\mathbb{P}}_n\left\{\frac{(1-T)Y}{1-\hat{\pi }(\text{X})}\right\},$$

$$\pi (\text{X},\mathbf{\theta })=(1-\delta -\epsilon )\varphi ({\mathbf{\beta }}^T\text{X})+\epsilon$$

$${\hat{\tau }}_P={\mathbb{P}}_n\left\{\frac{TY}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}-{\mathbb{P}}_n\left\{\frac{(1-T)Y}{1-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\},$$

$${\hat{\tau }}_{LDP}={\left({\mathbb{P}}_n\left[\frac{T}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\left\{1-\frac{C_1}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}\right]\right)}^{-1}{\mathbb{P}}_n\left[\frac{TY}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\left\{1-\frac{C_1}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}\right]$$

$$-{\left({\mathbb{P}}_n\left[\frac{1-T}{1-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\left\{1-\frac{C_0}{1-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}\right]\right)}^{-1}{\mathbb{P}}_n\left[\frac{(1-T)Y}{1-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\left\{1-\frac{C_0}{1-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}\right]$$

$$C_1=\frac{{\mathbb{P}}_n\left\{\frac{T-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}}{{\mathbb{P}}_n\left[{\left\{\frac{T-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}{\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}}^2\right]},C_0=-\frac{{\mathbb{P}}_n\left\{\frac{T-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}{1-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}}{{\mathbb{P}}_n\left[{\left\{\frac{T-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}{1-\pi (\text{X},{\hat{\mathbf{\theta }}}_n)}\right\}}^2\right]},$$

$$\sup \left\{z\in \mathbb{R}:P(|{\mathbf{\beta }}^T\mathbf{X}-z|<r)>0,for\; allr>0\right\}=+\infty .$$

$$\pi (\text{X})=\pi (\text{X},{\mathbf{\theta }}_0)={\epsilon }_0+\left(1-{\delta }_0-{\epsilon }_0\right)\left\{1+exp({\mathbf{\beta }}_0^T\text{X})\right\}{,}^{-1}$$

$$Y={\nu }_0+{\nu }_1{X}_1+{\nu }_2{X}_2+{\nu }_3{X}_3+{\nu }_{4}T+{\xi }_1{V}_1+{\xi }_2{V}_2+{\xi }_3{V}_3+Z,$$

$$\Sigma_{0}=\Sigma_{1}=A=\left(\begin{array}[]{cccc}1&0.5&-0.5&-0.5\\ 0.5&1&-0.5&-0.5\\ -0.5&-0.5&1&0.5\\ -0.5&-0.5&0.5&1\end{array}\right)$$

$$MSE(\hat{\tau })=\frac{1}{Nrep}\sum _{i=1}^{Nrep}{\left(\hat{\tau }-\tau \right)}^2.$$

$$\widehat{SE}=\sqrt{\frac{1}{Nboot-1}\sum _{i=1}^{Nboot}{\left({\hat{\tau }}_i-\overline{\hat{\tau }}\right)}^2}$$

$$\left[\hat{\tau }-1.96\widehat{SE},\hat{\tau }+1.96\widehat{SE}\right]$$

Full text

Keywords: observational studies, positivity, inverse probability weighting, average treatment effect, missing data.

1 Introduction

In the last 20 years, inverse probability weighted estimators (IPW) have attracted considerable attention in the statistical community. Among other things, they are used for estimating a population mean $E(Y)$ from an incomplete data set under the missing at random (MAR) assumption. MAR establishes that the variable of interest $Y$ and the response indicator $A$ are conditionally independent given an always observed vector $\mathbf{X}$ of covariates. See References [1, 2, 3, 4]. In the causal framework IPW estimators are used for estimating the average effect of a binary treatment on a scalar outcome under the assumption of no unmeasured confounders in observational studies. See References [5, 6, 7, 8, 9, 10].

In the missing data setting, the propensity score is defined as the response probability given the vector of covariates $\mathbf{X}$ whereas, in the causal context, the propensity score is the conditional probability of treatment assignment given a set of measured baseline covariates. IPW estimators are essentially weighted means of observed responses in which the weights are determined as the inverse of the estimated propensity scores. The aim of these weights is to compensate for the missing responses.

One of the conditions required to derive the asymptotic properties of IPW estimators is the strict positivity condition (also known as the positiviy assumption). It states, in the missing data setting, that the propensity score is bounded away from zero and, in the causal context, that the propensity score is bounded away from zero and one (see References [2, 4, 7, 8, 10]). Besides theoretical issues, the violation of the strict positivity condition causes the estimates to be very unstable and to have a large variability. See also References [11, 12, 13].

Most users of IPW procedures assume models for the propensity score that are usually incompatible with the positivity condition. For instance, this is the case of the generalized linear models (GLM), when some of the covariates are continuous. Despite this incompatibility, GLM are the most popular models considered in the field.

In this work we present a slight modification to the GLM originally postulated for the propensity score that is compatible with the strict positivity condition. To explore plausible applications of the new model, we revise a real example in the causal context. We consider the data collected by Tager et. al [14, 15], who investigate the effects of cigarette smoking on children’s pulmonary function. They study the effect of parental smoking on the pulmonary function of their children as well as the effect of direct smoking by the children themselves. Besides, in his book, Rosner [16] performs another analysis of the mentioned data. The author considers the forced expiratory volume (FEV) of a group of children as a response variable as well as their height, age, sex and a binary variable indicating whether or not they smoke. Kahn [17] studies this data using linear regression. In this work we study the effect of smoking in the forced expiratory volume of children by estimating the average treatment effect. This is an example of an observational study in which, under the assumption of no unmeasured confounders, the average treatment effect can be estimated using IPW procedures. The validity of no unmeasured confounders is beyond the scope of this work, and thus, the estimates here presented for the ATE are included only to illustrate the use of the new model for the propensity score. Since the proportion of smokers among the younger children is very small, the corresponding propensity scores take on values near zero and therefore the classical IPW estimators may have a large variance. In situations such as this, we expect that the models we propose in this work will provide more stable estimators than the classical IPW procedures.

The paper is structured as follows: in Section 2 we review methods based on propensity scores for estimating a population mean from incomplete data and we present our model for the propensity score in that missing data context. In Section 3 we review methods for estimating average treatment effects that use weighting by the inverse of the probability of treatment and present our model in the causality context. In Section 4 we describe the computational method we used to estimate the parameters indexing the new model. In Section 5 we study the performance of our model through its impact in IPW estimators by means of a Monte Carlo study. As an illustration, in Section 6, we apply these methods to estimate the effect of smoking in children.

2 Missing Data

Assume that we are interested in estimating the mean value of a scalar response $Y$, based on a sample composed by an always observed vector $\textbf{X}\in\mathbb{R}^{p}$ of covariates while the response of interest is missing by happenstance for some subjects. We will assume that data are missing at random (MAR) [18]. This means that the missing mechanism is not related to the response of interest and it is only related to X, the observed vector of covariates. Let $A$ be a binary variable indicating whether $Y$ is observed or not, namely, $A=1$ if $Y$ is observed and $A=0$ if $Y$ is missing. MAR establishes that

[TABLE]

$\pi(\boldsymbol{X})$ is known in the literature as the propensity score [19]. Up to integrability conditions, under the MAR assumption, $\mu=E(Y)$ can be represented in terms of the distribution of the observed data as $E(Y)=E\left\{AY/\pi(\textbf{X})\right\}$. This representation invites us to estimate $\mu$ by

[TABLE]

where $\hat{\pi}(\textbf{X})$ is a consistent estimator of $\pi(\textbf{X})$, and $\mathbb{P}_{n}$ is the empirical mean operator, i.e. $\mathbb{P}_{n}V=n^{-1}\sum_{i=1}^{n}V_{i}$ [20]. Observed responses corresponding to low values of the estimated propensity score are highly weighted since they should compensate for the high missing rate associated to such a level of covariates. For more details see [21, 1]. Different ways of estimating $\pi(\textbf{X})$ will lead to different estimators for $\mu$; parametric and nonparametric estimators of the propensity score have been considered. Robins, Rotnitzky and Zhao [1] presented a class of estimators for $\mu$, which containes $\hat{\mu}(\hat{\pi})$. These estimators are asymptotically normal when $\pi$ follows a parametric model. Little and An [3] proposed to estimate the propensity by fitting a spline to the logistic regression of the missing-data indicator $A$ on $\boldsymbol{X}$. Kang and Shafer [4] presented $\hat{\mu}(\hat{\pi})$, where $\pi$ follows a linear logistic regression model. More generally, parametric models are often assumed for the propensity score and, in particular, generalized linear models (GLM) play a prominent role among them. These models postulate that

[TABLE]

where $\phi$ is a strictly increasing cumulative distribution function and $\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$. Taking $\phi(u)=1/\{1+exp(-u)\}$ results in the linear logistic model, which is one of the most popular choices in the literature.

Most of the asymptotic results established for $\hat{\mu}(\hat{\pi})$, defined in (2), require the strict positivity condition, which establishes that $P\left\{\pi(\textbf{X})\geq\varepsilon\right\}=1$, for some $\varepsilon>0$ (see References [2, 4]). In practice, when this condition does not hold, values of $\hat{\pi}(\boldsymbol{X}_{i})$ close to zero may arise causing $\hat{\mu}(\hat{\pi})$ to be a very unstable estimator, having a large variability ( see also Reference [11]).

Unfortunately, generalized linear models prevent the validity of the strict positivity condition, except when $\boldsymbol{\beta}_{0}^{T}\textbf{X}$ is bounded from below. Thus, strict positivity is typically violated as far as continuous covariates are included in the vector $\boldsymbol{X}$.

In this work, we attempt to conciliate between the strict positivity condition and the most popular parametric models used for the propensity score. To do so, we slightly perturb the original parametric model postulated by the practitioner by incorporating an explicit lower bound for the propensity score. Namely, model (3) is replaced by

[TABLE]

where $\varepsilon_{0}\in(0,1)$ and $\boldsymbol{\beta}_{0}\in\mathbb{R}^{p}$. This model contemplates the validity of the strict positivity condition and thus we call it strictly positive propensity score (SPPS) model for the missing mechanism. Let $\boldsymbol{\theta}=(\varepsilon,\boldsymbol{\beta}^{T})^{T}$, with $\varepsilon\in(0,1)$ and $\boldsymbol{\beta}\in\mathbb{R}^{p}$. Define

[TABLE]

Assume that for some $\boldsymbol{\theta}_{0}\in\mathbb{R}^{p+1}$, $\pi(\textbf{X})=\pi(\textbf{X},\boldsymbol{\theta}_{0})$. We propose to estimate the mean of $Y$ by

[TABLE]

where $\hat{\boldsymbol{\theta}}_{n}$ is the maximum likelihood estimator of $\boldsymbol{\theta}_{0}$ under model (4).

In order to identify the parameter $\boldsymbol{\theta}_{0}$ from the distribution of $(\textbf{X},A)$, the following assumptions are required:

1. A.1

The function $\phi$ is injective. 2. A.2

The distribution of the observed covariates is not concentrated at a hyperplane. 3. A.3

For all $\boldsymbol{\beta}\in\mathbb{R}^{p}$ such that $A\,|\,\mathbf{X}\sim A\,|\,\boldsymbol{\beta}^{T}\mathbf{X}$, the support of $\boldsymbol{\beta}^{T}\mathbf{X}$ is unbounded from below:

[TABLE]

While the first two assumptions are the classic requirements to obtain the identifiability of the parameter given in a GLM model as (3), the last assumption guarantees the identifiability of the parameters involved in the model introduced in this work.

Assume that A.1, A.2 and A.3 hold. Under regularity conditions, it can be proved that $\hat{\boldsymbol{\theta}}_{n}$ is a consistent estimator of $\boldsymbol{\theta}_{0}$, and therefore $\hat{\mu}_{P}\,$ is a consistent estimator of $\mu$.

3 Average Treatment effect

Causal Inference is a second field where inverse probability weighted procedures play a crucial role. Consider, for instance, a dichotomous treatment variable $T$, where $T=1$ represents an active treatment and $T=0$ means that a control is assigned. The potential outcomes framework, introduced by Neyman [22] and Rubin [23], is used to quantify the effect of the treatment on some response of interest, whenever this difference is different from zero. To do so, two potential outcomes (or counterfactual variables) $Y^{(0)}$ and $Y^{(1)}$ are defined to represent the outcome variable of interest that would be seen if an individual were to receive the treatment and the control, respectively. We are interested in estimating the average treatment effect (ATE), defined as the difference between the mean values of the potential outcomes: $\tau=E(Y^{(1)})-E(Y^{(0)})$. $E(Y^{(1)})$ (respectively $E(Y^{(0)})$) represents the hypothetical mean of the response for the population of individuals where all of them are assumed to receive treatment (respectively control). So, the difference between these means may be considered a resultant of the treatment, meaning that it has a causal effect on the response of interest, whenever this difference is different from zero.

The potential outcomes $Y^{(1)}$ and $Y^{(0)}$ constitute an artificial contraption that allows us to conceptualize what we mean by causality. Only one of these variables is observed in each individual and it is related to the observed response through the consistency assumption which establishes (see Reference [24]) that, if an individual follows the treatment ($T=1$), then the potential outcome $Y^{(1)}$ is precisely his observed outcome. It also establishes the same thing for the control case. Therefore, under the consistency assumption, the observed outcome $Y$ is related to the counterfactual variables through the identity $Y=TY^{(1)}+(1-T)Y^{(0)}$. Thus, estimating the average treatment effect is a missing data problem, since one of the counterfactuals ($Y^{(1)}$ or $Y^{(0)}$) is missing for each individual. To identify $\tau$ we assume that a vector $\boldsymbol{X}$ with all possible counfounders is observed at each subject. This means that potential outcomes $Y^{(1)}$ and $Y^{(0)}$ are conditionally independent of the treatment exposure $T$ given $\boldsymbol{X}$:

[TABLE]

This assumption is known in the literature as strongly ignorable treatment assignment or no unmeasured confounders (see Reference [19]). Assuming (6), we have that

[TABLE]

where now the propensity score $\pi(\boldsymbol{X})$ is defined by $\pi(\boldsymbol{X})=P(T=1\mid\boldsymbol{X})$. These representations of $E\left(Y^{(1)}\right)$ and $E\left(Y^{(0)}\right)$ immediately suggest the estimator for the ATE proposed by Rosenbaum [5], given by

[TABLE]

where $\hat{\pi}(\textbf{X})$ is a consistent estimator of $\pi(\textbf{X})$. Lunceford and Davidian [7], following the general framework of [1] and the theory of M-estimation (see Reference [25]), presented large-sample theoretical properties of this estimator when $\pi$ follows a linear logistic model. Yao, Sun and Wang [9] postulate a generalized linear model for the propensity score, as in (3), and estimate $\pi(\boldsymbol{X})$ with $\phi\left(\hat{\boldsymbol{\beta}}_{n}^{T}\boldsymbol{X}\right)$, where $\hat{\boldsymbol{\beta}}_{n}$ denotes the MLE under model (3). Non parametric versions of $(\ref{est.gen})$ are presented and analysed by Hirano, Imbens and Ridder [6], who propose splines in order to estimate $\pi(\textbf{X})$. Their estimator achieves the semiparametric efficiency bound established by Hahn [26]. Both $\hat{\pi}(\boldsymbol{X})^{-1}$ and $\{1-\hat{\pi}(\boldsymbol{X})\}^{-1}$ are now involved in the estimator presented in (8). Thus, in this scenario, the strict positivity assumption is restated in terms of a lower and an upper bound for the propensity score, which is now assumed to be bounded away from zero and one, in the sense that $P\left\{\varepsilon\leq\pi(\boldsymbol{X})\leq 1-\delta\right\}=1$, for some $\varepsilon$ and $\delta$ greater than zero. This assumption, also known as (existence of) overlap in the covariate distributions, is usually assumed to derive the asymptotic properties of $\hat{\tau}(\hat{\pi})$, for the different estimators $\hat{\pi}$ of $\pi$ considered in the literature. See References [2, 6, 7, 8]. Besides, the lack of overlap leads to an erratic behavior of $\hat{\tau}(\hat{\pi})$, making the precise estimation of $\tau$ difficult. To deal with this issue, some authors propose trimmed versions of (8) (see References [12, 8, 13, 10]). As already mentioned in Section 2, generalized linear models for the propensity score are typically incompatible with the strict positivity assumption. Thus, using the same ideas developed in the missing data setting, for $\boldsymbol{\theta}=(\varepsilon,\delta,\boldsymbol{\beta}^{T})^{T}\in\boldsymbol{\theta}\subseteq\mathbb{R}^{p+2}$, where $\varepsilon,\delta\in(0,1)$ with $\varepsilon+\delta<1$ and $\boldsymbol{\beta}\in\mathbb{R}^{p},$ we propose the following model

[TABLE]

and assume that $\pi(\textbf{X})=\pi(\textbf{X},\boldsymbol{\theta}_{0})$, for some $\boldsymbol{\theta}_{0}\in\boldsymbol{\theta}.$ Model (9) intends to preserve as much as possible from the original family postulated for the propensity ($\phi(\boldsymbol{\beta}^{T}\textbf{X})$); the proposed modification contemplates the strict positivity condition by the inclusion of $\varepsilon$ and $\delta$. Model (9) will be called strictly positive propensity score (SPPS) for treatment assignment, or simply SPPS model if the context is clear. So, the estimator for the ATE that we consider under model (9) for the propensity score is given by

[TABLE]

where $\hat{\boldsymbol{\theta}}_{n}$ is the maximum likelihood estimator of $\boldsymbol{\theta}_{0}$ under model (9).

Lunceford and Davidian [7] propose a modified IPW estimator for $\tau$, in which the weights in $\hat{\tau}(\hat{\pi})$ are redefined in such a way that they not do not take on extremely large values. The authors show that this modification gives rise to a more stable procedure. We apply the same modification to the weights computed assuming the SPPS model and obtain the following estimator for the ATE, which combines the proposal of Lunceford and Davidian [7] and the one proposed in this work.

[TABLE]

[TABLE]

where

[TABLE]

and $\hat{\boldsymbol{\theta}}_{n}$ is the maximum likelihood estimator of $\boldsymbol{\theta}_{0}$ under model (9).

In this setting, in order to identify the parameter $\boldsymbol{\theta}_{0}$ and in addition to A.1, A.2 and A.3 we need the following assumption.

1. A.4

For all $\boldsymbol{\beta}\in\mathbb{R}^{p}$ such that $A\,|\,X\sim A\,|\,\boldsymbol{\beta}^{T}\mathbf{X},$ the support of $\boldsymbol{\beta}^{T}\mathbf{X}$ is unbounded from above:

[TABLE]

Under regularity conditions and assumptions A.1, A.2, A.3 and A.4 , it can be proved that $\boldsymbol{\theta}_{0}$ is identifiable and consistently estimated by $\hat{\boldsymbol{\theta}}_{n}$, and so, $\hat{\tau}_{P}$ is a consistent estimator of $\tau$.

4 Computational Method for Fitting the SPPS Model

In this section we describe the computational method used for the maximum likelihood estimation of the $p+1$ parameters in the SPPS model for the missing mechanism and the $p+2$ parameters in the SPPS model for treatment assignment. We describe in detail the method we propose for estimation in the treatment assignment case, giving in between brackets the necessary modifications for the missing data setting. Given $(\mathbf{X}_{1},T_{1})\dots(\mathbf{X}_{n},T_{n})$ a random sample from model (9) described in Section 3, we estimate $\varepsilon,\delta$ and $\boldsymbol{\beta}$ iteratively with the following procedure: consider that $T_{i}|\mathbf{X}_{i}$ follows a Bernoulli distribution with $P(T_{i}=1|\mathbf{X}_{i})=\pi(\textbf{X}_{i},\boldsymbol{\theta}_{0})$ where $\pi$ is given by equation (9). Estimate iteratively $\boldsymbol{\beta}$ assuming $\varepsilon$ and $\delta$ are known, $\varepsilon$ assuming $\boldsymbol{\beta}$ and $\delta$ are known, and $\delta$ assuming $\boldsymbol{\beta}$ and $\varepsilon$ are known.

1. Step 1

Substitute [math] for $\varepsilon$ and $\delta$ in (9) and compute an initial estimate of $\boldsymbol{\beta}$ by maximum likelihood. Call this estimate $\hat{\boldsymbol{\beta}}_{0}$. 2. Step 2

Compute the fitted values of the model considered in Step 1, $\hat{\pi_{i}}=\phi\left(\hat{\boldsymbol{\beta}}_{0}^{T}\mathbf{X}_{i}\right)$. 3. Step 3

Compute initial estimates of $\varepsilon$ and $\delta$ by $\hat{\varepsilon}_{1}=min\left(\hat{\pi}_{i}\right)$ and $\hat{\delta}_{1}=1-max\left(\hat{\pi}_{i}\right)$. (In the missing data case set $\hat{\delta}_{1}=0$). 4. Step 4

Substitute $\hat{\varepsilon}_{1}$ for $\varepsilon$ and $\hat{\delta}_{1}$ for $\delta$ in (9). Then recompute an estimate of $\boldsymbol{\beta}$ by maximum likelihood. Call this estimate $\hat{\boldsymbol{\beta}_{1}}$. 5. Step 5

Substitute $\hat{\boldsymbol{\beta}}_{1}$ for $\boldsymbol{\beta}$ and $\hat{\delta}_{1}$ for $\delta$ in (9) and estimate $\varepsilon$ by maximum likelihood. This estimate is the new $\hat{\varepsilon_{1}}$. 6. Step 6

Substitute $\hat{\boldsymbol{\beta}}_{1}$ for $\boldsymbol{\beta}$ and $\hat{\varepsilon}_{1}$ for $\varepsilon$ in (9) and estimate $\delta$ by maximum likelihood. This estimate is the new $\hat{\delta_{1}}$. (Skip this step in the missing data case). 7. Step 7

Repeat steps 4 to 6 until convergence (4 and 5 in the missing data case). Return $\hat{\boldsymbol{\pi}}=(\hat{\pi}_{1},\dots,\hat{\pi}_{n})$.

This iterative method converges rapidly in the case of the missing data setting. In the causality setting it also converges rapidly if the real values of $\varepsilon$ and $\delta$ are such that $\varepsilon+\delta$ is not too large. If this is not the case, the modified link function in the GLM in Step 4 can be very flat, growing very slowly from $\hat{\varepsilon}_{1}$ to $1-\hat{\delta}_{1}$, which are very similar. When this happens the GLM method in R usually does not converge, except for very large samples. To solve this problem, in the causality setting, we change Step 3 by the following

1. Step 3 ’

Compute initial estimators of $\varepsilon$ and $\delta$ by $\hat{\varepsilon}_{1}=min\left(\hat{\pi}_{i}\right)$ and $\hat{\delta}_{1}=1-max\left(\hat{\pi}_{i}\right)$. If $\hat{\varepsilon}_{1}+\hat{\delta}_{1}>0.6$, return $\hat{\boldsymbol{\pi}}=(\hat{\pi}_{1},\dots,\hat{\pi}_{n})$, otherwise proceed to step 4.

We remark that Step 3’ implies that, if after fitting a regular GLM we have that the range of the fitted values is not very large, namely smaller than $0.4$, then we keep the regular GLM; there is no need for the SPPS model.

The code that defines the necessary functions for computing the ATE and IPW estimators is available under request.

5 Monte Carlo Study

In this section we report the results of a Monte Carlo study we performed in order to assess the advantages of considering the proposed model over the traditional approach. In practice, several covariates will be available for modeling the propensity score. In order to investigate performance in a realistic setting, Lunceford and Dadivian [7], carried out simulations involving continuous and discrete variables, some of them associated with both treatment exposure and potential response and others associated with potential responses but not treatment exposure. We generate variables as they did in one of their proposed scenarios. We consider covariates $\textbf{X}=(X_{1},X_{2},X_{3},V_{1},V_{2},V_{3}),\,$ a binary variable $T\,$ and an outcome $Y\,$ such that the variable $\,T\,$ follows a Bernoulli distribution with

[TABLE]

where $\boldsymbol{\theta}_{0}=(\varepsilon_{0},\delta_{0},\boldsymbol{\beta_{0}})$ and $\boldsymbol{\beta_{0}}=(0,0.6,-0.6,0.6,0,0,0)$.

We remark that, when $\varepsilon_{0}=\delta_{0}=0$ this SPPS model reduces to the usual linear logistic model. Different settings of $\varepsilon_{0}$ and $\delta_{0}$ will be chosen. The outcome is generated by

[TABLE]

where $Z$ is a normal standard variable independent of $\mathbf{X}$ and $T$, $\left(\nu_{0},\nu_{1},\nu_{2},\nu_{3},\nu_{4}\right)=\left(0,-1,1,-1,2\right)^{T}$ and $\xi=(-1,1,1).$ The joint distribution of $\mathbf{X}$ is specified by taking $X_{3}$ with Bernoulli(0.2) distribution, and then generating $V_{3}$ as Bernoulli with $P(V_{3}=1|X_{3})=0.75X_{3}+0.25(1-X_{3}).$ Conditional on $X_{3},$ the vector $\left(X_{1},V_{1},X_{2},V_{2}\right)$ is generated as multivariate normal $\mathcal{N}(\rho_{X_{3}},\Sigma_{X_{3}}),$ where $\rho_{0}=\left(1,1,-1,-1\right),\,\rho_{1}=\left(-1,-1,1,1\right),$ and

[TABLE]

We generate $Nrep=1000$ samples of size $n=1000$ following the described model, under which the real value of the average treatment effect is $\tau=2$. For each sample we compute four estimators, namely $\hat{\tau}_{O}$, $\hat{\tau}_{P}$, $\hat{\tau}_{LD}$ and $\hat{\tau}_{PLD}$. $\hat{\tau}_{P}$ and $\hat{\tau}_{PLD}$ are defined in (10 ) and (11) respectively; both of them use the SPPS model (9) introduced in this work for the propensity score. On the other hand, $\hat{\tau}_{O}$ is the original IPW estimator, defined by modeling the propensity score with the linear logistic model. That is to say, $\hat{\tau}_{O}=\hat{\tau}(\hat{\pi})$ given in (8) assuming that $\pi(\textbf{X})=\phi(\boldsymbol{\beta}_{1}^{T}\textbf{X})$, with $\phi(u)=1/(1+exp(-u))$ for some $\boldsymbol{\beta}_{1}\in\mathbb{R}^{7}.$ Finally, $\hat{\tau}_{LD}$ denotes the estimator proposed in Lunceford and Davidian (2004). This estimator is defined as in equation (8) but using a linear logistic regression model for the propensity score. For each estimator $\hat{\tau}$, we compute an empirical mean square error with the following formula

[TABLE]

The results of the Monte Carlo study are reported in Table 1. In these simulations we can see that, when the true model is the linear logistic model, both our proposed estimators give a smaller MSE than the corresponding classical ones. On the other hand, if the samples are generated following the SPPS model with $\epsilon_{0}$ and $\delta_{0}$ both greater than zero, then $\hat{\tau}_{PLD}$ remains very small, as compared to the other proposals, followed sometimes by $\hat{\tau}_{P}$ and sometimes by $\hat{\tau}_{LD}$. We remark that in almost all the situations considered our proposed estimators give better results than the corresponding classical ones. In the remaining situations, namely when samples are generated with large values of $\epsilon_{0}$ and $\delta_{0}$, the mean square errors are similar and all the estimators give good results. When $\epsilon_{0}+\delta_{0}$ is greater than $0.6$, all four estimators are equal or almost equal.

6 Example: Children’s FEV Data

This data set contains measurements of the forced expiratory volume (FEV) of 654 children and teenagers aged 3 to 19 years, together with their height, age, sex and and a binary variable indicating whether or not they smoke. It is basically the data set considered in [16], which has been included in the R package covreg [27] with slight modifications. In order to estimate the average smoking effect in the FEV in this population, we consider only children aged 9 or older since there are not any smokers among the younger children in this data set.

We wish to emphasize that this example is included to show the use of the proposed model in a causal inference context. The confiability of results depend on the no unmeasured counfounders assumption, and the discussion on its validity given the observed covariates is beyond the scope of this work. This example is therefore considered only as an illustration of the application of the proposed methods.

All four estimates yield a negative average treatment effect but differ in the size of this effect. As a means to determine the significance of these differences, we computed $95\%$ normal confidence intervals were we estimated the standard error based on 1000 bootstrap samples. To obtain the bootstrap samples, we first break the data set in two, smokers (65 subjects) and non-smokers (374 subjects). Then we resample separately 65 observations from the smokers group and 374 observations from the non smokers group. Let $\hat{\tau}_{i}$ be the estimator of the ATE based on the $ith$ sample and $\overline{\hat{\tau}}$ the mean of $\hat{\tau}_{1},\dots\hat{\tau}_{Nboot}$. The bootstrap estimator of the standard error of each estimator $\hat{\tau}$ is computed by

[TABLE]

and the $95\%$ bootstrap confidence interval based on $\hat{\tau}$ is computed as

[TABLE]

As reported in Table 2, the obtained intervals show that the ATE is not significant. The conclusion is the same for all the estimators considered.

In Table 2 we report the estimates, together with their corresponding standard errors and confidence interval for the ATE. We remark that the standard error of our proposed estimator $\hat{\tau}_{P}$ is slightly smaller than the standard errors of the classical ones.

Acknowledgement This research was partly supported by grants 20020150200110BA and 20020130100279BA from Universidad de Buenos Aires.