Propensity score weighting for causal inference with clustered data

TL;DR

This paper introduces a calibration technique for propensity score weighting that accounts for unmeasured cluster-specific confounders in observational studies with clustered data, improving causal effect estimation.

Contribution

It proposes a robust, model-free propensity score weighting estimator that handles unmeasured cluster effects and can be combined with sampling weights for complex survey data.

Findings

Estimator outperforms competitors in simulations

Effective in estimating effects of school BMI screening

Handles unmeasured cluster confounders robustly

Abstract

Propensity score weighting is a tool for causal inference to adjust for measured confounders in observational studies. In practice, data often present complex structures, such as clustering, which make propensity score modeling and estimation challenging. In addition, for clustered data, there may be unmeasured cluster-specific variables that are related to both the treatment assignment and the outcome. When such unmeasured cluster-specific confounders exist and are omitted in the propensity score model, the subsequent propensity score adjustment may be biased. In this article, we propose a calibration technique for propensity score estimation under the latent ignorable treatment assignment mechanism, i.e., the treatment-outcome relationship is unconfounded given the observed covariates and the latent cluster effects. We then provide a consistent propensity score weighting estimator of the average treatment effect when the propensity score and outcome follow generalized linear mixed effects models. The proposed propensity score weighting estimator is attractive, because it does not require specification of functional forms of the propensity score and outcome models, and therefore is robust to model misspecification. The proposed weighting method can be combined with sampling weights for an integrated solution to handle confounding and sampling designs for causal inference with clustered survey data. In simulation studies, we show that the proposed estimator is superior to other competitors. We estimate the effect of School Body Mass Index Screening on prevalence of overweight and obesity for elementary schools in Pennsylvania.

Figures

Tables

Table 1. Table 1: Simulation results: bias, variance ( × 10 − 3 absent superscript 10 3 \times 10^{-3} ) and coverage ( % percent \% ) of 95 % percent 95 95\% confidence intervals based on 1 , 000 1 000 1,000 Monte Carlo samples; the outcome is linear and logistic linear mixed effects model and the propensity score is logistic, probit or complementary log-log (C-loglog).

	$(m,n_e)=(50,50)$				$(m,n_e)=(100,30)$				$(m,n_e)=(30,100)$
Method	bias	var	cvg		bias	var	cvg		bias	var	cvg
Scenario 1: Linear outcome & Logistic propensity score
${\hat{\tau }}_{\mathrm{simp}}$	-0.37	22	27.4		-0.38	12	8.7		-0.38	35	42.3
${\hat{\tau }}_{\mathrm{fix}}$	-0.01	36	95.6		0.00	21	95.6		-0.01	42	95.2
${\hat{\tau }}_{\mathrm{ran}}$	0.14	26	90.2		0.21	14	64.6		0.07	37	94.7
${\hat{\tau }}_{\mathrm{cal}}$	0.01	26	94.5		0.02	11	95.1		0.00	33	95.6
Scenario 2: Linear outcome & Probit propensity score
${\hat{\tau }}_{\mathrm{simp}}$	-0.29	16	34.4		-0.08	9	2.3		-0.22	30	65.6
${\hat{\tau }}_{\mathrm{fix}}$	0.08	35	90.3		-0.10	19	4.5		0.12	69	90.4
${\hat{\tau }}_{\mathrm{ran}}$	0.24	28	73.9		-0.07	16	29.9		0.21	60	85.5
${\hat{\tau }}_{\mathrm{cal}}$	0.01	22	94.9		0.01	11	95.4		0.00	33	94.6
Scenario 3: Linear outcome & C-loglog propensity score
${\hat{\tau }}_{\mathrm{simp}}$	-0.21	20	62.0		-0.21	10	41.2		-0.22	30	65.6
${\hat{\tau }}_{\mathrm{fix}}$	0.12	48	88.8		0.12	36	82.7		0.12	69	90.4
${\hat{\tau }}_{\mathrm{ran}}$	0.29	38	69.1		0.36	22	32.5		0.21	60	85.5
${\hat{\tau }}_{\mathrm{cal}}$	0.00	21	95.3		0.00	10	95.1		0.00	33	94.6
Scenario 4: Logistic outcome & Logistic propensity score
${\hat{\tau }}_{\mathrm{simp}}$	-0.11	100	9.1		-0.11	540	0.5		-0.11	160	20.5
${\hat{\tau }}_{\mathrm{fix}}$	-0.11	44	0.3		-0.11	38	0.1		-0.11	39	0.1
${\hat{\tau }}_{\mathrm{ran}}$	-0.09	33	1.3		-0.08	21	0.5		-0.10	34	0.3
${\hat{\tau }}_{\mathrm{cal}}$	0.01	74	96.3		0.01	55	95.2		0.01	74	95.9
Scenario 5: Logistic outcome & Probit propensity score
${\hat{\tau }}_{\mathrm{simp}}$	-0.08	58	13.1		-0.08	34	2.3		-0.08	81	25.3
${\hat{\tau }}_{\mathrm{fix}}$	-0.10	93	6.9		-0.10	85	4.5		-0.10	73	3.8
${\hat{\tau }}_{\mathrm{ran}}$	-0.08	67	23.0		-0.07	48	29.9		-0.09	61	8.3
${\hat{\tau }}_{\mathrm{cal}}$	0.01	89	94.7		0.01	65	95.4		0.01	84	95.0
Scenario 6: Logistic outcome & C-loglog propensity score
${\hat{\tau }}_{\mathrm{simp}}$	-0.06	0.3	3.2		-0.06	0.2	1.0		-0.06	0.2	3.7
${\hat{\tau }}_{\mathrm{fix}}$	-0.05	0.5	44.6		-0.05	0.5	43.6		-0.05	0.5	43.0
${\hat{\tau }}_{\mathrm{ran}}$	-0.03	0.5	95.4		-0.03	0.4	97.3		-0.03	0.4	92.8
${\hat{\tau }}_{\mathrm{cal}}$	-0.01	0.7	95.5		0.00	0.6	95.8		-0.01	0.7	95.2

Table 2. Table 2: Balance Check

		simple	fixed	random	calibration
	Cluster 1	1.68	-0.22	0.68	0.20
	Cluster 2	1.21	0.10	-0.41	0.10
$X_1$	Cluster 3	1.75	-0.02	0.99	0.02
	Cluster 4	0.86	-0.04	-1.05	0.02
	Cluster 5	-0.36	0.37	-1.39	0.33
	Whole Pop	1.28	-0.02	-0.02	0
	Cluster 1	0.48	0.02	0.30	0.03
	Cluster 2	0.43	0.13	-0.01	0.14
$X_2$	Cluster 3	0.73	0.01	0.46	0.02
	Cluster 4	0.18	-0.08	-0.34	-0.07
	Cluster 5	-0.57	-0.39	-1.53	-0.44
	Whole Pop	0.39	-0.003	-0.001	0

Table 3. Table 3: Results: estimate, variance estimate (ve) based on 500 500 500 bootstrap replicates, and 95 % percent 95 95\% confidence interval (c.i.)

	estimate	ve	$95\%$ c.i.
simple	8.78	2.11	(5.94, 11.63)
fixed	0.47	0.44	(-0.83, 1.77)
random	0.52	0.44	(-0.77, 1.82)
calibration	0.53	0.39	(-0.71, 1.76)

Equations

$$\begin{gathered} \prod _{i=1}^m\prod _{j=1}^{n_i}{\left\{f_1(Y_{ij}\mid X_{ij},U_i)\mathrm{pr}(A_{ij}=1\mid X_{ij},U_i)\right\}}^{A_{ij}} \\ \times {\left[f_0(Y_{ij}\mid X_{ij},U_i)\{1-\mathrm{pr}(A_{ij}=1\mid X_{ij},U_i)\}\right]}^{1-A_{ij}}, \end{gathered}$$

$$g\{{\mu }_{ij}(a)\}=X_{ij}^{\mathrm{T}}{\beta }_a+U_i,$$

$$\mathrm{pr}(A_{ij}=1\mid X_{ij},U_i)=h(X_{ij}^{\mathrm{T}}\gamma +U_i)\equiv e(X_{ij},U_i;\eta ),$$

$${\hat{\tau }}_{\mathrm{IPTW}}(\nu ,\eta )=\frac{1}{n}\sum _{i=1}^m\sum _{j=1}^{n_i}\left\{\frac{A_{ij}{Y}_{ij}}{e(X_{ij},U_i;\eta )}-\frac{(1-A_{ij})Y_{ij}}{1-e(X_{ij},U_i;\eta )}\right\}.$$

$$E\left\{\frac{A}{e(X,U)}\left(\begin{array}[]{c}X\\ U\end{array}\right)\right\}=E\left[E\left\{\frac{A}{e(X,U)}\mid X,U\right\}\left(\begin{array}[]{c}X\\ U\end{array}\right)\right]=E\left\{\left(\begin{array}[]{c}X\\ U\end{array}\right)\right\},$$

$$E\left\{\frac{1-A}{1-e(X,U)}\left(\begin{array}[]{c}X\\ U\end{array}\right)\right\}=E\left[E\left\{\frac{1-A}{1-e(X,U)}\mid X,U\right\}\left(\begin{array}[]{c}X\\ U\end{array}\right)\right]=E\left\{\left(\begin{array}[]{c}X\\ U\end{array}\right)\right\}.$$

$$\sum _{i=1}^m\sum _{j=1}^{n_i}\frac{A_{ij}}{{\hat{e}}_{ij}}{X}_{ij}=\sum _{i=1}^m\sum _{j=1}^{n_i}\frac{1-A_{ij}}{1-{\hat{e}}_{ij}}{X}_{ij}=\sum _{i=1}^m\sum _{j=1}^{n_i}{X}_{ij},$$

$$\sum _{j=1}^{n_i}\frac{A_{ij}}{{\hat{e}}_{ij}}=\sum _{j=1}^{n_i}\frac{1-A_{ij}}{1-{\hat{e}}_{ij}}=\sum _{j=1}^{n_i}1=n_i,(i=1,\dots ,m).$$

$$\sum _{i=1}^m\sum _{j=1}^{n_i}G({\alpha }_{ij},d_{ij})=\sum _{i=1}^m\sum _{j=1}^{n_i}{\alpha }_{ij}\log \frac{{\alpha }_{ij}}{d_{ij}},$$

$$\begin{gathered} {\alpha }_{ij}({\lambda }_1,{\lambda }_2)=\frac{n_i{A}_{ij}{d}_{ij}\exp ({\lambda }_1{X}_{ij}{A}_{ij})}{{\sum }_{j=1}^{n_i}{A}_{ij}{d}_{ij}\exp ({\lambda }_1{X}_{ij}{A}_{ij})} \\ +\frac{n_i(1-A_{ij})d_{ij}\exp \{{\lambda }_2{X}_{ij}(1-A_{ij})\}}{{\sum }_{j=1}^{n_i}(1-A_{ij})d_{ij}\exp \{{\lambda }_2{X}_{ij}(1-A_{ij})\}}, \end{gathered}$$

$$\hat{Q}(\lambda_{1},\lambda_{2})=\left(\begin{array}[]{c}\hat{Q}_{1}(\lambda_{1},\lambda_{2})\\ \hat{Q}_{2}(\lambda_{1},\lambda_{2})\end{array}\right)\\ =\left(\begin{array}[]{c}n^{-1}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\left\{A_{ij}\alpha_{ij}(\lambda_{1},\lambda_{2})-1\right\}X_{ij}\\ n^{-1}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\left\{(1-A_{ij})\alpha_{ij}(\lambda_{1},\lambda_{2})-1\right\}X_{ij}\end{array}\right)=0.$$

$${\hat{\tau }}_{\mathrm{IPTW}}=\frac{1}{n}\sum _{i=1}^m\sum _{j=1}^{n_i}\left\{\frac{A_{ij}{Y}_{ij}}{{\hat{e}}_{ij}}-\frac{(1-A_{ij})Y_{ij}}{1-{\hat{e}}_{ij}}\right\}.$$

$${\mu }_1(U_i)=\frac{\int q_1(x,U_i)E\{Y_{ij}(1)\mid x,U_i\}f(x)dx}{\int q_1(x,U_i)f(x)dx},q_1(X_{ij},U_i)=E\left(\frac{A_{ij}}{{\hat{e}}_{ij}}-1\mid X_{ij},U_i\right),$$

$$\displaystyle E\left[\frac{1}{n}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\left(\frac{A_{ij}}{\hat{e}_{ij}}-1\right)\left\{Y_{ij}(1)-\mbox{$\mu$}_{1}(U_{i})\right\}\right]$$

$$E\left\{\frac{1}{n}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\left(\frac{A_{ij}}{\hat{e}_{ij}}-1\right)Y_{ij}(1)\right\}=E\left\{\frac{1}{n}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\left(\frac{A_{ij}}{\hat{e}_{ij}}-1\right)\mbox{$\mu$}_{1}(U_{i})\right\}=0,$$

$$\begin{gathered} E\left(\frac{1}{n}\sum _{i=1}^m\sum _{j=1}^{n_i}\frac{A_{ij}}{{\hat{e}}_{ij}}{Y}_{ij}\right)=E\left\{\frac{1}{n}\sum _{i=1}^m\sum _{j=1}^{n_i}\frac{A_{ij}}{{\hat{e}}_{ij}}{Y}_{ij}(1)\right\} \\ =E\left\{\frac{1}{n}\sum _{i=1}^m\sum _{j=1}^{n_i}{Y}_{ij}(1)\right\}=E\{Y(1)\}. \end{gathered}$$

$$E\left(\frac{1}{n}\sum _{i=1}^m\sum _{j=1}^{n_i}\frac{1-A_{ij}}{1-{\hat{e}}_{ij}}{Y}_{ij}\right)=E\{Y(0)\}.$$

$$V_1^{-1/2}({\hat{\tau }}_{\mathrm{IPTW}}-\tau )\to \mathcal{N}(0,1),$$

$$B_1$$

$$B_2$$

$${\hat{B}}_1$$

$${\hat{B}}_2$$

$$\hat{V}({\hat{\tau }}_{\mathrm{IPTW}})=\frac{1}{n}\left(\frac{1}{m-1}\sum _{i=1}^m({\hat{\tau }}_i-{\hat{\tau }}_{\mathrm{IPTW}}{)}^2+\frac{1}{m}\sum _{i=1}^m{\hat{V}}_i\right).$$

$${\hat{\tau }}_{\mathrm{IPTW}}=\frac{1}{N}\sum _{i\in S_I}\sum _{j=1}^{n_i}{\omega }_{ij}\left\{\frac{A_{ij}{Y}_{ij}}{e(X_{ij},U_i)}-\frac{(1-A_{ij})Y_{ij}}{1-e(X_{ij},U_i)}\right\}.$$

$$E({\hat{\tau }}_{\mathrm{IPTW}})$$

$$\sum _{i\in S_I}\sum _{j=1}^{n_i}{\omega }_{ij}\frac{A_{ij}}{{\hat{e}}_{ij}}{X}_{ij}$$

$$\sum _{j=1}^{n_i}{\omega }_{ij}\frac{A_{ij}}{{\hat{e}}_{ij}}$$

$$\begin{gathered} {\alpha }_{ij}({\lambda }_1,{\lambda }_2)=\frac{N_i{A}_{ij}{d}_{ij}\exp ({\lambda }_1{X}_{ij}{A}_{ij})}{{\sum }_{j=1}^{n_i}{\omega }_{ij}{A}_{ij}{d}_{ij}\exp ({\lambda }_1{X}_{ij}{A}_{ij})} \\ +\frac{N_i(1-A_{ij})d_{ij}\exp \{{\lambda }_2{X}_{ij}(1-A_{ij})\}}{{\sum }_{j=1}^{n_i}{\omega }_{ij}(1-A_{ij})d_{ij}\exp \{{\lambda }_2{X}_{ij}(1-A_{ij})\}}, \end{gathered}$$

$$\hat{Q}(\lambda_{1},\lambda_{2})=\left(\begin{array}[]{c}\hat{Q}_{1}(\lambda_{1},\lambda_{2})\\ \hat{Q}_{2}(\lambda_{1},\lambda_{2})\end{array}\right)\\ =\left(\begin{array}[]{c}N^{-1}\sum_{i\in S_{I}}\sum_{j=1}^{n_{i}}\omega_{ij}\left\{A_{ij}\alpha_{ij}(\lambda_{1},\lambda_{2})-1\right\}X_{ij}\\ N^{-1}\sum_{i\in S_{I}}\sum_{j=1}^{n_{i}}\omega_{ij}\left\{(1-A_{ij})\alpha_{ij}(\lambda_{1},\lambda_{2})-1\right\}X_{ij}\end{array}\right)=0.$$

$${\hat{\tau }}_{\mathrm{IPTW}}=\frac{1}{N}\sum _{i\in S_I}\sum _{j=1}^{n_i}{\omega }_{ij}\left\{\frac{A_{ij}{Y}_{ij}}{{\hat{e}}_{ij}}-\frac{(1-A_{ij})Y_{ij}}{1-{\hat{e}}_{ij}}\right\}.$$

$$V_2^{-1}({\hat{\tau }}_{\mathrm{IPTW}}-\tau )\to \mathcal{N}(0,1),$$

$$B_1$$

Full text

Keywords: Calibration; Inverse probability weighting; Mixed effects model; Survey sampling; Unmeasured confounding.

1 Introduction

Observational studies are often used to infer causal effects in medical and social science studies. In observational studies, there often is confounding by indication: some covariates are predictors of both the treatment and outcome. One implication is that the covariate distributions differ between the treatment groups. Under the assumption of ignorable treatment assignment and that all confounders are observed, the causal effect of treatments can be obtained by comparing the outcomes for units from different treatment groups, adjusting for the observed confounders. Rosenbaum and Rubin (1983) further claimed the central role of the propensity score, and showed that adjusting for the propensity score removes confounding bias. An extensive literature thereafter proposed a number of estimators based on the propensity score. In particular, propensity score weighting can be used to create a weighted population where the covariate distributions are balanced between the treatment groups, on average. Therefore, the comparison between the weighted outcomes has a causal interpretation; see Imbens and Rubin (2015) for a textbook discussion.

Propensity score weighting has been mainly developed and applied in settings with independently and identically distributed (i.i.d.) data. However, in many research area, data often present complex structures, such as clustering. For a motivating example, we examine the 2007–2010 body mass index (BMI) surveillance data from Pennsylvania Department of Health to estimate the effect of School Body Mass Index Screening (SBMIS) on the annual overweight and obesity prevalence in elementary schools in Pennsylvania. The data set includes $493$ school districts in Pennsylvania, which are clustered by two factors: location (rural, suburban, and urban), and population density (low, median, and high). In this data set, $63\%$ of schools implemented SBMIS, and the percentages of schools implemented SBMIS across the clusters are from $45\%$ to $70\%$, indicating cluster-level heterogeneity of treatment. Moreover, even if we collect a rich set of unit-level covariates, there may be unobserved cluster effects that are related to both the treatment and outcome. In our motivating example, we have school-level covariates including the baseline prevalence of overweight and obesity and percentage of reduced and free lunch. However, certain key health factors, such as accessibility to and quality of care, socioeconomic   and environmental variables, which can be very different across clusters, are preceivably important factors for children’s obesity rate and implementing prevention screening strategy. Unfortunately, these cluster-specific confounders are not available. When such unmeasured confounders exist and are omitted in the propensity score model, the subsequent analysis may be biased. Although important for empirical practice, much less work has been done for causal inference with clustered data subject to unmeasured cluster-specific confounding.

The goal of this article is to develop a novel propensity score weighting method for causal inference with clustered data in the presence of unmeasured cluster-level confounders. Natural models for the propensity score and outcome are generalized linear mixed effects models, where cluster-level confounding is captured via cluster random effect terms. Prior to our work, Li et al. (2013) investigated the performance of the propensity score weighting estimators under generalized linear fixed/mixed effects models for the propensity score and outcome. However, their approach requires correct specification of functional forms of the propensity score and outcome models. In this article, we provide a robust construction of inverse propensity score weights under the latent ignorable treatment assignment mechanism, i.e., the treatment-outcome relationship is unconfounded given the observed covariates and the latent cluster effects. The key insight is based on the central role of the propensity score in balancing the covariate distributions between the treatment groups. For propensity score estimation, we then adopt the calibration technique and impose balancing constraints for moments of the observed and latent cluster-specific confounders between the treatment groups. These constraints eliminate confounding biases. Under certain regularity conditions, the propensity score weighting estimator of the average treatment effect is consistent when the propensity score and outcome follow generalized linear mixed effects models. The proposed propensity score weighting estimator does not require correct specification of functional forms of the treatment and outcome models, and therefore is robust to model misspecification.

2 Basic Setup

2.1 Observed data structure

To fix the ideas, we first focus on clustered data. The extension to clustered survey data will be addressed in Section 5. Suppose that the sample consists of $m$ clusters, and cluster $i$ includes $n_{i}$ units. Denote the sample size by $n=\sum_{i=1}^{m}n_{i}$. For unit $j$ in cluster $i$, we observe a $p$-dimensional vector of pre-treatment variables $X_{ij}$, which may include the observed individual and cluster characteristics, a binary treatment variable $A_{ij}\in\{0,1\}$, with [math] and $1$ being the labels of control and active treatments, respectively, and lastly an outcome variable $Y_{ij}$.

2.2 Models and assumptions

We use the potential outcomes framework (Rubin, 1974). Assume that each unit has two potential outcomes: $Y_{ij}(0)$, the outcome that would be realized, possibly contrary to the fact, had the unit received the control treatment, and $Y_{ij}(1)$, the outcome that would be realized, possibly contrary to the fact, had the unit received the active treatment. This notation implicitly assumes that there is no interference between units and no versions of each treatment (Rubin, 1978). The observed outcome is $Y_{ij}=Y_{ij}(A_{ij})$ (Rubin, 1974). Suppose that $\{A_{ij},X_{ij},Y_{ij}(0),Y_{ij}(1):i=1,\ldots,m;j=1,\ldots n_{i}\}$ follows an infinite super-population distribution. Our goal is to estimate the average treatment effect, $\tau=E\{Y_{ij}(1)-Y_{ij}(0)\}$, where the expectation is taken with respect to the super-population distribution. In the binary case, $\tau$ is called the causal risk difference.

The fundamental problem is that not all potential outcomes can be observed for each unit in the sample; only one potential outcome, the outcome corresponding to the treatment the unit actually followed, can be observed. With unstructured i.i.d. data, Rubin (1974) described the following assumption for identifying average treatment effect.

Assumption 1** (Ignorability)**

\{Y_{ij}(0),Y_{ij}(1)\}\mathbin{\mathchoice{\hbox to0.0pt{\hbox{\set@color\displaystyle\perp}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color\displaystyle\perp}}{\hbox to0.0pt{\hbox{\set@color\textstyle\perp}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color\textstyle\perp}}{\hbox to0.0pt{\hbox{\set@color\scriptstyle\perp}\hss}\kern 2.36812pt{}\kern 2.36812pt\hbox{\set@color\scriptstyle\perp}}{\hbox to0.0pt{\hbox{\set@color\scriptscriptstyle\perp}\hss}\kern 1.63437pt{}\kern 1.63437pt\hbox{\set@color\scriptscriptstyle\perp}}}A_{ij}\mid X_{ij}.**

Assumption 1 indicates that all confounders are observed. For clustered data, confounding may vary across clusters and related to some cluster-specific variables that are not always observable. In these cases, Assumption 1 does not hold. We assume a cluster-specific latent effect $U_{i}$ that summarizes the effect of unobserved cluster-level confounders, and consider the following modified ignorability assumption.

Assumption 2** (Latent ignorability)**

\{Y_{ij}(0),Y_{ij}(1)\}\mathbin{\mathchoice{\hbox to0.0pt{\hbox{\set@color\displaystyle\perp}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color\displaystyle\perp}}{\hbox to0.0pt{\hbox{\set@color\textstyle\perp}\hss}\kern 3.46875pt{}\kern 3.46875pt\hbox{\set@color\textstyle\perp}}{\hbox to0.0pt{\hbox{\set@color\scriptstyle\perp}\hss}\kern 2.36812pt{}\kern 2.36812pt\hbox{\set@color\scriptstyle\perp}}{\hbox to0.0pt{\hbox{\set@color\scriptscriptstyle\perp}\hss}\kern 1.63437pt{}\kern 1.63437pt\hbox{\set@color\scriptscriptstyle\perp}}}A_{ij}\mid X_{ij},U_{i}.**

Under Assumption 2, $\mathrm{pr}\{A_{ij}=1\mid X_{ij},U_{i},Y_{ij}(0),Y_{ij}(1)\}=\mathrm{pr}(A_{ij}=1\mid X_{ij},U_{i}),$ which is the propensity score. We assume the numbers of units received the active and control treatments are nonzero in each cluster; otherwise, there exist units in some clusters for which we can not estimate the average treatment effect without extrapolation assumptions.

Assumption 3** (Positivity)**

There exist constants $\underline{e}$ and $\bar{e}$ such that, with probability $1$, $0<\underline{e}<\mathrm{pr}(A_{ij}=1\mid X_{ij},U_{i})<\bar{e}<1$.

Under Assumption 2, we can write the conditional distribution of $\{(A_{ij},Y_{ij}):i=1,\ldots,m;j=1,\ldots n_{i}\}$ given $\{(X_{ij},U_{i}):i=1,\ldots,m;j=1,\ldots n_{i}\}$ as

[TABLE]

where $f_{a}(\cdot\mid X_{ij},U_{i})$ is a conditional distribution of $Y_{ij}(a)$ given $(X_{ij},U_{i})$, for $a=0,1$.

There are two different model specifications regarding the cluster-level effects. The fixed effects model treats $U_{i}$ as fixed but unknown parameters across clusters. In this fixed-effects approach, the treatment assignment mechanism is an ignorable process, which complies with Assumption 1 given that $X_{ij}$ stacks all observed confounders and cluster-specific dummy variables. On the other hand, the random effects model treats $U_{i}$ as random and i.i.d. drawn from a distribution. The difference between the two modeling strategies has been addressed in both statistics and econometrics literature; see, e.g., Baltagi (1995) and Wooldridge (2002). Briefly, there are two main considerations: one on statistical consideration and the other on logic consideration. First, if the number of clusters is relatively large, the estimator from the fixed effects approach becomes inconsistent in propensity score estimation (Wallace and Hussain, 1969). In this case, the random effects approach is preferred. Second, the fixed effects approach does not make distributional assumptions of the cluster-specific effects; whereas, the random effects approach assumes that $U_{i}$ is random and i.i.d. drawn from a distribution. To justify this random effects assumption, we can use the exchangeability criterion of Chamberlain (1984). If the $U_{i}$’s can be randomly permuted to ensure exchangeability; a version of de Finetti’s theorem implies then that they are i.i.d. as draws from an appropriately defined distribution. In the case of evaluating the causal effect of SBMIS on children’s obesity, this assumption implies that unobserved cluster characteristics that influence both implementing SBMIS and children’s obesity are not correlated with school characteristics that are included in the models.

In this article, we assume that the $U_{i}$’s are random variables and independently follow a certain distribution. Figure 1 provides a directed acyclic graph (Pearl, 2009) that implies the dependence of variables under our assumptions in cluster $i$.

We now posit generalized linear mixed effects models for $f_{a}(Y_{ij}\mid X_{ij},U_{i})$ and $\mathrm{pr}(A_{ij}=1\mid X_{ij},U_{i})$. To be specific, we suppose that $Y_{ij}(a)$ follows a generalized linear mixed effects model with a random effect $U_{i}$ as

[TABLE]

where $\mu_{ij}(a)=E\{Y_{ij}(a)\mid X_{ij},U_{i}\}$, $g(\cdot)$ is the link function, and $\beta_{a}$ is a $p$-dimensional vector of fixed effects of $X_{ij}$. Similarly, we assume that $A_{ij}$ given $(X_{ij},U_{i})$ follows a generalized linear mixed effects model as

[TABLE]

where $h(\cdot)$ is the inverse link function, and $\gamma$ is a $q$-dimensional vector of parameters.

2.3 Inverse probability of treatment weighting estimator

To estimate the average treatment effect $\tau$, let $\nu=(U_{1},\ldots,U_{m})$ denote the vector of random effects. The inverse propensity score or probability of treatment weighting (IPTW) estimator of $\tau$ can be expressed as

[TABLE]

Under Assumptions 2 and 3, if the propensity score is known, it is straightforward to verify that $\hat{\tau}_{\mathrm{IPTW}}(\nu,\eta)$ is unbiased for $\tau$. Moreover, if it is unknown but depends only on fixed parameters, $\hat{\tau}_{\mathrm{IPTW}}(\nu,\eta)$ with the consistently estimated propensity score is asymptotically unbiased for $\tau$. The challenge with clustered data is that $\hat{\tau}_{\mathrm{IPTW}}(\nu,\eta)$ depends on a growing number of unobserved random effects with the number of clusters. To resolve this issue, there are several options:

(i)

Weight based on predicted random effects; i.e., based on the generalized linear mixed effects model in (2), predict the propensity score as $e(X_{ij},\hat{U}_{i}^{\mathrm{ran}};\hat{\eta})$, where $\hat{U}_{i}^{\mathrm{ran}}$ is the mode of a predictive distribution for $U_{i}$ given the observed $A_{ij}$ and $X_{ij}$, and $\hat{\eta}$ is the maximum likelihood estimator of $\eta$.

(ii)

Weight based on estimated fixed effects; i.e., treat the $U_{i}$’s in model (2) as fixed effects, and estimate the propensity score as $e(X_{ij},\hat{U}_{i}^{\mathrm{fix}};\hat{\eta})$, where $\hat{U}_{i}^{\mathrm{fix}}$ and $\hat{\eta}$ are maximum likelihood estimators.

Let $\hat{\tau}_{\mathrm{IPTW}}(\nu,\eta)$ in (3) be denoted as $\hat{\tau}_{\mathrm{ran}}$ or $\hat{\tau}_{\mathrm{fix}}$ when the propensity score is predicted under option (i) or estimated under option (ii), respectively. The two approaches suffer several drawbacks. Firstly, to obtain $\hat{\tau}_{\mathrm{ran}}$ often involves numerical integration, and therefore can be computationally heavy. Secondly, the predicted value of the propensity score does not guarantee the balance of covariate distributions between the treatment groups, due to the shrinkage of random effects toward zero. Lastly, $\hat{\tau}_{\mathrm{fix}}$ does not yield a consistent estimator for $\tau$ for small $n_{i}$ (Skinner et al., 2011).

3 Proposed methodology

3.1 The new IPTW estimator

To motivate our estimation of the propensity score, we note

[TABLE]

and

[TABLE]

(4) and (5) clarify the central role of the propensity score in balancing the covariate distributions between the treatment groups in the super-population. For simplicity of presentation, let $e_{ij}$ be the propensity score for unit $j$ in cluster $i$. We consider the propensity score estimate $\hat{e}_{ij}$ that satisfies the empirical version of (4) and (5):

[TABLE]

[TABLE]

To obtain the propensity score estimate that achieves (6) and (7), we use the calibration technique in the following steps:

Step  1.

Obtain an initial propensity score estimate $\hat{e}_{ij}^{0}$ under some working propensity score model, e.g. a logistic linear model fitted to $(A_{ij},X_{ij})$. This in turn provides an initial set of inverse propensity score weights, $\mathbb{W}^{0}=\{d_{ij};i=1,\ldots,m,j=1,\ldots,n_{i}\}$, where $d_{ij}=1/e_{ij}^{0}$ if $A_{ij}=1$ and $d_{ij}=1/(1-e_{ij}^{0})$ if $A_{ij}=0$.

Step  2.

Modify the initial set of weights $\mathbb{W}^{0}$ to a new set of weights $\mathbb{W}=\{\alpha_{ij};i=1,\ldots,m,j=1,\ldots,n_{i}\}$ by minimizing the Kullback-Leibler distance (Kullback and Leibler, 1951) of $\mathbb{W}^{0}$ and $\mathbb{W}$:

[TABLE]

subject to (6) and (7). By Lagrange Multiplier, the minimizer of (8) subject to (6) and (7) is

[TABLE]

where $(\lambda_{1},\lambda_{2})^{\mathrm{\scriptscriptstyle T}}$ is the solution to the following equation

[TABLE]

Step  3.

Obtain the propensity score estimate as $\hat{e}_{ij}=\alpha_{ij}(\hat{\lambda}_{1},\hat{\lambda}_{2})^{-A_{ij}}\{1-\alpha_{ij}(\hat{\lambda}_{1},\hat{\lambda}_{2})\}^{-1+A_{ij}}$.

Finally, our proposed IPTW estimator is

[TABLE]

Remark 1

Calibration has been used in many scenarios. In survey sampling, calibration is widely used to integrate auxiliary information in estimation or handle nonresponse; see, e.g., Wu and Sitter (2001), Chen et al. (2002), Kott (2006), Chang and Kott (2008) and Kim et al. (2016). In causal inference, calibration has been used such as in Constrained Empirical Likelihood (Qin and Zhang, 2007), Entropy Balancing (Hainmueller, 2012), Inverse Probability Tilting (Graham et al., 2012), and Covariate Balance Propensity Score of Imai and Ratkovic (2014). Chan et al. (2015) showed that estimation of average treatment effects by empirical balancing calibration weighting can achieve global efficiency. However, all these works were developed in settings with i.i.d. variables and they assumed that there are no unmeasured confounders. Our article is the first to use calibration for causal inference with unmeasured cluster-specific confounders.

Remark 2

In Step 2 of the calibration algorithm, different distance functions, other than the Kullback-Leibler distance, can be considered. For example, if we choose $G(\alpha_{ij},d_{ij})=d_{ij}(\alpha_{ij}/d_{ij}-1)^{2}$, then the minimum distance estimation leads to generalized regression estimation (Park and Fuller, 2012) of the $\alpha_{ij}$’s. If we choose $G(\alpha_{ij},d_{ij})=-d_{ij}\log(\alpha_{ij}/d_{ij})$, then it leads to empirical likelihood estimation (Newey and Smith, 2004). We use the Kullback–Leibler distance function, which leads to exponential tilting estimation (Kitamura and Stutzer, 1997, Imbens et al., 1998, Schennach et al., 2007). The advantage of using the Kullback-Leibler distance is that the resulting weights are always non-negative. Also, with Kullback-Leibler distance, the calibration constraint (7) can be built into a closed form expression for the weights, and thus avoiding solving a large number of equations. This reduces the computation burden greatly, when there is a large number of clusters.

4 Main results

To discuss the asymptotic properties of the proposed estimator, we assume that the number of clusters increases with $n$, i.e., $m\rightarrow\infty$, as $n\rightarrow\infty$ , and that the cluster sample sizes satisfy the condition that $\sup_{1\leq i\leq m}n_{i}=O(n^{1/2})$. Under this asymptotic framework, the number of clusters increases but some of the cluster sample sizes may remain small. We also impose certain regularity conditions specified in the Appendix. Denote $A\cong B$ as $A=B+o_{p}(1)$, where the reference distribution is the super-population model.

To show the consistency of $\hat{\tau}_{\mathrm{IPTW}}$, we first introduce a cluster-specific mean function:

[TABLE]

where $f(x)$ is the density of $X$. The key then is to note

[TABLE]

as $m\rightarrow\infty$, which follows from the definition of $\mu_{1}(U_{i})$ in (12). (13) implies that

[TABLE]

where zero follows from the constraint (7). Under Assumption 2 and (14), it follows

[TABLE]

Similarly, we establish

[TABLE]

Combining (15) and (16), we have $E(\hat{\tau}_{\mathrm{IPTW}})\cong E\{Y(1)-Y(0)\}=\tau,$ which yields the asymptotic unbiasedness of $\hat{\tau}_{\mathrm{IPTW}}$. Under certain regularity conditions as in the Appendix, we then have $\mathrm{plim}_{n\rightarrow\infty}\hat{\tau}_{\mathrm{IPTW}}=\tau$.

It is important to note that to the logistic model is only a working model. The proposed estimator $\hat{\tau}_{\mathrm{IPTW}}$ does not require specification of this working model to be true. (6) and (7) play the key role for the unbiasedness of $\hat{\tau}_{\mathrm{IPTW}}$. Therefore, $\hat{\tau}_{\mathrm{\mathrm{IPTW}}}$ is robust to specification of this working propensity score model.

Theorem 1

Suppose that Assumptions 2 and 3, and the regularity conditions specified in the Appendix hold, and that the outcome and propensity score follow the generalized linear mixed effects models in (1) and (2). Suppose further that the number of clusters $m$ and the cluster sample sizes $n_{i}$, for $i=1,\ldots m$, satisfy the condition that $m\rightarrow\infty$, as $n\rightarrow\infty$, and $\sup_{1\leq i\leq m}n_{i}=O(n^{1/2})$. Then, the proposed propensity score weighting estimator in (11), subject to constraints (6) and (7), satisfies

[TABLE]

in distribution, as $n\rightarrow\infty$, where $V_{1}=\mathrm{var}\left(n^{-1}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\tau_{ij}\right),$ with $\tau_{ij}=\{\alpha_{ij}(\lambda_{1}^{*},\lambda_{2}^{*})A_{ij}(Y_{ij}-B_{1}^{\mathrm{\scriptscriptstyle T}}X_{ij})+B_{1}^{\mathrm{\scriptscriptstyle T}}X_{ij}\}-\{\alpha_{ij}(\lambda_{1}^{*},\lambda_{2}^{*})(1-A_{ij})(Y_{ij}-B_{2}^{\mathrm{\scriptscriptstyle T}}X_{ij})+B_{2}^{\mathrm{\scriptscriptstyle T}}X_{ij}\}$,

[TABLE]

and $(\lambda_{1}^{*},\lambda_{2}^{*})^{\mathrm{\scriptscriptstyle T}}$ satisfies $E\{\hat{Q}(\lambda_{1}^{*},\lambda_{2}^{*})\}=0$ with $\hat{Q}(\lambda_{1},\lambda_{2})$ defined in (10).

The proof is relegated to the Appendix. We now discuss variance estimation. Let $\hat{\tau}_{ij}=\alpha_{ij}(\hat{\lambda}_{1},\hat{\lambda}_{2})\{A_{ij}(Y_{ij}-\hat{B}_{1}^{\mathrm{\scriptscriptstyle T}}X_{ij})-(1-A_{ij})(Y_{ij}-\hat{B}_{2}^{\mathrm{\scriptscriptstyle T}}X_{ij})\}+(\hat{B}_{1}-\hat{B}_{2})^{\mathrm{\scriptscriptstyle T}}X_{ij}$, where

[TABLE]

Let $\hat{\tau}_{i}=n_{i}^{-1}\sum_{j=1}^{n_{i}}\hat{\tau}_{ij}$ and $\hat{V}_{i}=(n_{i}-1)^{-1}\sum_{j=1}^{n_{i}}(\hat{\tau}_{ij}-\hat{\tau}_{i})^{2}.$ The variance estimator can be constructed as

[TABLE]

5 Extension to clustered survey data

Clustered data often arise in survey sampling. In complex surveys, the challenge is to take design information or design weights into account when developing propensity score methods for causal inference. In this section, we extend the proposed propensity score weighting estimator to clustered survey data. Consider a finite population $\mathcal{F}_{N}$ with $M$ clusters and $N_{i}$ units in the $i$th cluster, where $N=\sum_{i=1}^{M}N_{i}$ denotes the population size. We assume that in the finite population, $\{A_{ij},X_{ij},Y_{ij}(0),Y_{ij}(1)\}$ follows the superpopulation model $\xi$ as described in Section 2. We are interested in estimating the average treatment effect $\tau=E\{Y(1)-Y(0)\}$.

The sample is selected according to a two-stage cluster sampling design. Specifically, at the first stage, cluster $i$ is sampled with the first inclusion probability $\pi_{i}$, $i\in S_{I}$, where $S_{I}$ is the index set for the sampled clusters. Let $\pi_{ij}=\mathrm{pr}(i,j\in S_{I})$ be the second inclusion probability for clusters $i$ and $j$ being sampled. At the second stage, given that cluster $i$ was selected at the first stage, unit $j$ is sampled with conditional probability $\pi_{j|i}$, $j=1,\ldots,n_{i}$. The final sample size is $n=\sum_{i\in S_{I}}n_{i}$. The design weight for unit $j$ in cluster $i$ be $\omega_{ij}=(\pi_{i}\pi_{j|i})^{-1}$, which reflects the number of units for cluster $i$ in the finite population this unit $j$ represents. We assume that the design weights are positive and known throughout the sample. Also, let $\pi_{kl|i}$ be the second inclusion probability for units $k$ and $l$ being sampled given that cluster $i$ was selected. The second inclusion probabilities, the $\pi_{ij}$ and $\pi_{kl|i}$’s, are often used for variance estimation.

For clustered survey data, if the propensity score $e(X_{ij},U_{i})$ is known, we can express the IPTW estimator of $\tau$ as

[TABLE]

Let $E_{\xi}$ and $E_{p}$ denote expectation under the super-population model and the sampling design, respectively. It is easy to verify that

[TABLE]

In practice, since the propensity score $e(X_{ij},U_{i})$ is often unknown, (17) is not infeasible. To estimate the propensity score, we now require the propensity score estimate $\hat{e}_{ij}$ satisfy the following design-weighted moment constraints

[TABLE]

These moment constraints (18) and (19) are the sample version of (4) and (5), respectively.

To obtain the propensity score estimate that achieves (18) and (19), we use the calibration technique in the following steps:

Step  1.

Obtain an initial propensity score estimate $\hat{e}_{ij}^{0}$ under some working propensity score model, e.g. a logistic linear model fitted to $(A_{ij},X_{ij})$, each unit weighted by the design weight $\omega_{ij}$. This in turn provides an initial set of inverse propensity score weights, $\mathbb{W}^{0}=\{d_{ij};i=1,\ldots,m,j=1,\ldots,n_{i}\}$, where $d_{ij}=1/e_{ij}^{0}$ if $A_{ij}=1$ and $d_{ij}=1/(1-e_{ij}^{0})$ if $A_{ij}=0$.

Step  2.

Modify the initial set of weights $\mathbb{W}^{0}$ to a new set of weights $\mathbb{W}=\{\alpha_{ij};i=1,\ldots,m,j=1,\ldots,n_{i}\}$ by minimizing $\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\omega_{ij}\alpha_{ij}\log(\alpha_{ij}/d_{ij}),$ subject to (18) and (19). By Lagrange Multiplier, $\alpha_{ij}$ can be obtained as

[TABLE]

where $(\lambda_{1},\lambda_{2})^{\mathrm{\scriptscriptstyle T}}$ is the solution to the following equation

[TABLE]

Step  3.

Obtain the propensity score estimate as $\hat{e}_{ij}=\alpha_{ij}(\hat{\lambda}_{1},\hat{\lambda}_{2})^{-A_{ij}}\{1-\alpha_{ij}(\hat{\lambda}_{1},\hat{\lambda}_{2})\}^{-1+A_{ij}}$.

Finally, our proposed IPTW estimator is

[TABLE]

In the above procedure, the design wights are used in both the propensity score estimates and the weighting estimator.

We now consider the asymptotic property of $\hat{\tau}_{\mathrm{IPTW}}$ in (21). We use an asymptotic framework, where the sample size $n$ indexes a sequence of finite populations and samples (Fuller, 2009; Section $1.3$), such that the population size $N$ increases with $n$, but the cluster sample sizes $N_{i}$ may remain small. In addition, we have the following regularity conditions for the sampling mechanism.

Assumption 4

(i) The first-order inclusion probability $\pi_{i}\pi_{j|i}$ is positive and uniformly bounded in the sense that there exist positive constants $C_{1}$ and $C_{2}$ that do not depend on $N$, such that $C_{1}\leq\pi_{i}\pi_{j|i}Nn^{-1}\leq C_{2},$ for any $i$ and $j$; (ii) the sequence of Hotvitz-Thompson estimators $\hat{Y}_{\mathrm{HT}}=N^{-1}\sum_{i\in S_{I}}\sum_{j=1}^{n_{i}}\omega_{ij}y_{i}$ satisfies $\mathrm{var}_{p}(\hat{Y}_{\mathrm{HT}})=O\left(n^{-1}\right)$ and $\left\{\mathrm{var}_{p}(\hat{Y}_{\mathrm{HT}})\right\}^{-1/2}(\hat{Y}_{\mathrm{HT}}-\bar{Y})\mid\mathcal{F}_{N}\rightarrow\mathcal{N}(0,1),$ in distribution, as $n\rightarrow\infty$, where $\bar{Y}=N^{-1}\sum_{i=1}^{M}\sum_{j=1}^{N_{i}}y_{i}$ is the population mean of $Y$, and the reference distribution is the randomization distribution generated by the sampling mechanism.

Sufficient conditions for the asymptotic normality of the Hotvitz-Thompson estimators are discussed in Chapter 1 of Fuller (2009).

Theorem 2

Suppose that Assumptions 2 and 3, and the regularity conditions specified in the Appendix hold, and that the outcome and propensity score follow generalized linear mixed effects models in (1) and (2). Suppose further that the sequence of finite populations and samples satisfy Assumption 4. Then, the proposed propensity score weighting estimator in (21), subject to constraints (18) and (19), satisfies

[TABLE]

in distribution, as $n\rightarrow\infty$, where $V_{2}=\mathrm{var}\left(N^{-1}\sum_{i\in S_{I}}\sum_{j=1}^{n_{i}}\omega_{ij}\tau_{ij}\right),$ with $\tau_{ij}=\{\alpha_{ij}(\lambda_{1}^{*},\lambda_{2}^{*})A_{ij}(Y_{ij}-B_{1}^{\mathrm{\scriptscriptstyle T}}X_{ij})+B_{1}^{\mathrm{\scriptscriptstyle T}}X_{ij}\}-\{\alpha_{ij}(\lambda_{1}^{*},\lambda_{2}^{*})(1-A_{ij})(Y_{ij}-B_{2}^{\mathrm{\scriptscriptstyle T}}X_{ij})+B_{2}^{\mathrm{\scriptscriptstyle T}}X_{ij}\}$,

[TABLE]

and $(\lambda_{1}^{*},\lambda_{2}^{*})^{\mathrm{\scriptscriptstyle T}}$ satisfies $E\{\hat{Q}(\lambda_{1}^{*},\lambda_{2}^{*})\}=0$ with $\hat{Q}(\lambda_{1},\lambda_{2})$ defined in (20).

For variance estimation of $\hat{\tau}_{\mathrm{IPTW}}$, let $\hat{\tau}_{ij}=\alpha_{ij}(\hat{\lambda}_{1},\hat{\lambda}_{2})\{A_{ij}(Y_{ij}-\hat{B}_{1}^{\mathrm{\scriptscriptstyle T}}X_{ij})-(1-A_{ij})(Y_{ij}-\hat{B}_{2}^{\mathrm{\scriptscriptstyle T}}X_{ij})\}+(\hat{B}_{1}-\hat{B}_{2})^{\mathrm{\scriptscriptstyle T}}X_{ij}$, where

[TABLE]

Let $\hat{\tau}_{i}=\sum_{j=1}^{n_{i}}\pi_{j|i}^{-1}\hat{\tau}_{ij}$ and

[TABLE]

The variance estimator can be constructed as

[TABLE]

6 A simulation study

We conduct simulation studies to evaluate the finite-sample performance of the proposed estimator. We first generate finite populations and then select a sample from each finite population using a two-stage cluster sampling design.

In the first setting, we specify the number of clusters in the population to be $M=10,000$, and the size of the $i$th cluster size $N_{i}$ to be the integer part of $500\times\exp(2+U_{i})/\{1+\exp(2+U_{i})\}$, where $U_{i}\sim\mathcal{N}(0,1)$. The cluster sizes range from $100$ to $500$. The potential outcomes are generated according to linear mixed effects models, $Y_{ij}(0)=X_{ij}+U_{i}+e_{ij}$ and $Y_{ij}(1)=X_{ij}+\tau+\tau U_{i}+e_{ij}$, where $\tau=2$, $X_{ij}\sim\mathcal{N}(0,1)$, $e_{ij}\sim\mathcal{N}(0,1)$, $U_{i}$, $X_{ij}$, and $e_{ij}$ are independent, for $i=1,\ldots,M$, $j=1,\ldots,N_{i}$. The parameter of interest is $\tau$. We consider three propensity score models, $\mathrm{pr}(A_{ij}=1\mid X_{ij};U_{i})=h(\gamma_{0}+\gamma_{1}U_{i}+X_{ij})$, with $h(\cdot)$ being the inverse logit, probit and complementary log-log link function, for generating $A_{ij}$. We set $(\gamma_{0},\gamma_{1})$ to be $(-0.5,1)$, $(-0.25,0.5)$ and $(-0.5,0.1)$ for the above three propensity score models, respectively. The observed outcome is $Y_{ij}=A_{ij}Y_{ij}(1)+(1-A_{ij})Y_{ij}(0)$. From each realized population, $m$ clusters are sampled by Probability-Proportional-to-Size (PPS) sampling with the measure of size $N_{i}$. So the first-order inclusion probability of selecting cluster $i$ is equal to $\pi_{i}=mN_{i}/\sum_{i=1}^{m}N_{i}$, which implicitly depends on the unobserved random effect. Once the clusters are sampled, the $n_{i}$ units in the $i$th selected cluster are sampled by Poison sampling with the corresponding first-order inclusion probability $\pi_{j|i}=n_{e}z_{ij}/(\sum_{j=1}^{M_{i}}z_{ij})$, where $z_{ij}=0.5$ if $e_{ij}<0$ and 1 if $e_{ij}>0$. With this sampling design, the units with $e_{ij}>0$ are sampled with a chance twice as big as the units with $e_{ij}<0$. We consider three combinations of $m$ and $n_{e}$: (i) $(m,n_{e})=(50,50)$; (ii) $(m,n_{e})=(100,30)$, representing a large number of small clusters; and (iii) $(m,n_{e})=(30,100)$, representing a small number of large clusters.

In the second setting, all data-generating mechanisms are the same as in the first setting, except that the potential outcomes are generated according to logistic linear mixed effects models, $Y_{ij}(0)\sim\mathrm{Bernoulli}(p_{ij}^{0})$ with $\mathrm{logit}(p_{ij}^{0})=X_{ij}+U_{i}$ and $Y_{ij}(1)\sim\mathrm{Bernoulli}(p_{ij}^{1})$ with $\mathrm{logit}(p_{ij}^{1})=X_{ij}+\tau+\tau u_{i}$. Moreover, in the 2-stage cluster sampling, $\pi_{j|i}=n_{e}z_{ij}/(\sum_{j=1}^{M_{i}}z_{ij})$, where $z_{ij}=0.5$ if $Y_{ij}=0$ and 1 if $Y_{ij}=1$. With this sampling design, the units with $Y_{ij}=1$ are sampled with a chance twice as big as the units with $Y_{ij}=0$.

We compare four estimators for $\tau$: (i) $\hat{\tau}_{\mathrm{simp}}$, the simple design-weighted estimator without propensity score adjustment; (ii) $\hat{\tau}_{\mathrm{fix}}$, the weighting estimator in (3) with the propensity score estimated by a logistic linear fixed effects model with a cluster-level main effect; (iii) $\hat{\tau}_{\mathrm{ran}}$, the weighting estimator in (3) with the propensity score estimated by a logistic linear mixed effects model where the cluster effect is random; (iv) $\hat{\tau}_{\mathrm{\mathrm{IPTW}}}$, the proposed estimator with calibrations (18) and (19). Table 1 shows biases, variances and coverages for $95\%$ confidence intervals from $1,000$ simulated data sets. The simple estimator shows large biases across difference scenarios, even adjusting for sampling design. This suggests that the covariate distributions are different between the treatment groups in the finite population, contributing to the bias. $\hat{\tau}_{\mathrm{fix}}$ works well under Scenario 1 with the linear mixed effects model for the outcome and the logistic linear mixed effects model for the propensity score; however, its performance is not satisfactory in other scenarios. Moreover, $\hat{\tau}_{\mathrm{fix}}$ shows the largest variance among the four estimators in most of scenarios. This is because for a moderate or large number of clusters, there are too many free parameters, and hence the propensity score estimates may not be stable. For $\hat{\tau}_{\mathrm{ran}}$, we assume that the cluster effect is random, which reduces the number of free parameters. As a result, $\hat{\tau}_{\mathrm{ran}}$ shows less variability than $\hat{\tau}_{\mathrm{fix}}$. Nonetheless, both $\hat{\tau}_{\mathrm{fix}}$ and $\hat{\tau}_{\mathrm{ran}}$ can not control the bias well. The proposed estimator shows small biases and good empirical coverages across all scenarios. Notably, to compute $\hat{\tau}_{\mathrm{IPTW}},$ we used a working model, a logistic linear model, to provide an initial set of weights. When the true propensity score is probit or complementary log-log model, $\hat{\tau}_{\mathrm{\mathrm{IPTW}}}$ still has small biases. This suggests that our proposed estimator is robust to the working model in our simulation settings.

7 An Application

We examine the 2007–2010 BMI surveillance data from Pennsylvania Department of Health to investigate the effect of School Body Mass Index Screening (SBMIS) on the annual overweight and obesity prevalence in elementary schools in Pennsylvania. Early studies have shown that SBMIS has been associated with increased parental awareness of child weight (Harris et al., 2009; Ebbeling et al., 2012). However, there have been mixed findings about the effect of screening on reducing prevalence of overweight and obesity (Harris et al., 2009; Thompson and Card-Higginson, 2009). The data set includes $493$ school districts in Pennsylvania. The baseline is the school year 2007. The schools are clustered by two factors: location (rural, suburban, and urban), and population density (low, median, and high). This results in five clusters: rural-low, rural-median, rural-high, suburban-high, and urban-high. Let $A=1$ if the school implemented SBMIS, and $A=0$ if the school did not. In this data set, $63\%$ of schools implemented SBMIS, and the percentages of schools implemented SBMIS across the clusters range from $45\%$ to $70\%$, indicating cluster-level heterogeneity of treatment. The outcome variable $Y$ is the annual overweight and obesity prevalence for each school district in the school year 2010. The prevalence is calculated by dividing the number with BMI$>85$th by the total number of students screened for each school district. For each school, we obtain school characteristics including the baseline prevalence of overweight and obesity ($X_{1}$), and percentage of reduced and free lunch ($X_{2}$).

For a direct comparison, the average difference of the prevalence of overweight and obesity for schools that implemented SBMIS and those that did not is $8.78\%$. This unadjusted difference in the prevalence of overweight and obesity ignores differences in schools and clusters. To take the cluster-level heterogeneity of treatment into account, we consider three propensity score models: (i) a logistic linear fixed effects model with linear predictors including $X_{1}$, $X_{2}$, and a fixed intercept for each cluster; (ii) a logistic linear mixed effects model with linear predictors including fixed effects $X_{1}$, $X_{2}$, and a random effect for each cluster; (iii) the proposed calibrated propensity score. Using the estimated propensity score, we estimate $\tau=E\{Y(1)-Y(0)\}$ by the weighting method.

Table 2 displays the standardized differences of means for $X_{1}$ and $X_{2}$ between the treated and control groups for each cluster and the whole population, standardized by the standard errors in the whole population. Without any adjustment, there are large differences in means for $X_{1}$ and $X_{2}$. For this specific data set, the three methods for modeling and estimating the propensity score are similar in balancing the covariate distributions between the treated and control groups. All three propensity score weighting methods improve the balance for $X_{1}$ and $X_{2}$. Table 3 displays point estimates and variance estimates based on $500$ bootstrap replicates. The simple estimator shows that the screening has a significant effect in reducing the prevalence of overweight and obesity. However, this may be due to confounders. After adjusting for the confounders, the screening does not have a significant effect. Given the different sets of assumptions for the different methods, this conclusion is reassuring.

8 Discussion

The IPTW estimator is not efficient in general. Semiparametric efficiency bounds for estimating the average treatment effects in the setting with i.i.d. random variables were derived by Hahn (1998). He showed that the efficient influence function for the average treatment effect depends on both the propensity score and the outcome model. An important implication is that combining the propensity score model and the outcome regression model can improve efficiency of the IPTW estimator. For clustered data, since the data are correlated through the random cluster effects, the efficiency theory established for the i.i.d. data is not applicable. It remains an interesting avenue for future research to develop the semiparametric efficiency theory for clustered data.

In this article, we assumed that there is no interference between units. This setup is not uncommon. In our application, the treatment was implemented school-wise. The potential outcomes for one school are likely to be unaffected by the treatments implemented at other schools, and therefore the assumption of no interference is likely to hold. However, in other settings, this assumption may not hold. A classical example is given in infectious diseases (Ross, 1916; Hudgens and Halloran, 2008), where whether one person becomes infected depends on who else in the population is vaccinated. Extension of our calibration estimation to take the interference structure into account in these settings is also an interesting topic for future research.

Acknowledgments

We would like to thank Peng Ding and Fan Li for insightful and fruitful discussions.

Appendix

Appendix A. Regularity conditions

Condition A1

(i) $\hat{\tau}_{\mathrm{IPTW}}(\lambda_{1},\lambda_{2})\rightarrow\tau$ in probability uniformly in a compact set $\mathcal{B}$ as $n\rightarrow\infty$; (ii) $\hat{Q}(\lambda_{1},\lambda_{2})\rightarrow Q(\lambda_{1},\lambda_{2})=E\{\hat{Q}(\lambda_{1},\lambda_{2})\}$ in probability uniformly in $\mathcal{B}$ as $n\rightarrow\infty$, and there exists a unique $(\lambda_{1}^{*},\lambda_{2}^{*})\in\mathcal{B}$ such that $Q(\lambda_{1}^{*},\lambda_{2}^{*})=0$; (iii) $\partial\hat{\tau}_{\mathrm{IPTW}}(\lambda_{1},\lambda_{2})/\partial(\lambda_{1},\lambda_{2})^{\mathrm{\scriptscriptstyle T}}$ and $\partial\hat{Q}(\lambda_{1},\lambda_{2})/\partial(\lambda_{1},\lambda_{2})^{\mathrm{\scriptscriptstyle T}}$ are continuous at $(\lambda_{1},\lambda_{2})$ in $\mathcal{B}$ almost surely; (iv)$E||X_{ij}||^{3}<\infty$, $E|Y_{ij}(0)|^{3}<\infty$, and $E|Y_{ij}(1)|^{3}<\infty$; (v) the matrix $E\left\{\partial\hat{Q}(\lambda_{1}^{*},\lambda_{2}^{*})/\partial(\lambda_{1},\lambda_{2})^{\mathrm{\scriptscriptstyle T}}\right\}$ is invertible.

Conditions A1 (i)–(iii) ensure that $\mathrm{plim}_{n\rightarrow\infty}(\hat{\lambda}_{1},\hat{\lambda}_{2})=(\lambda_{1}^{*},\lambda_{2}^{*})$ and

[TABLE]

which is similar to Corollary II.2 of Andersen and Gill (1982). Condition A1 (iv) is a moment condition for the central limit theorem.

Appendix B. Proof of Theorem 1

Write $\hat{\tau}_{\mathrm{IPTW}}(\lambda_{1},\lambda_{2})=n^{-1}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\alpha_{ij}(\lambda_{1},\lambda_{2})Y_{ij},$ where $\alpha_{ij}(\lambda_{1},\lambda_{2})$ is defined in (9). The proposed estimator is $\hat{\tau}_{\mathrm{IPTW}}(\hat{\lambda}_{1},\hat{\lambda}_{2})$, where $(\hat{\lambda}_{1},\hat{\lambda}_{2})$ satisfies $\hat{Q}(\hat{\lambda}_{1},\hat{\lambda}_{2})=0$, where $\hat{Q}(\lambda_{1},\lambda_{2})$ is defined in (10). Let $(\lambda_{1}^{*},\lambda_{2}^{*})$ satisfy $E\{Q(\lambda_{1}^{*},\lambda_{2}^{*})\}=0$.

Under Conditions A1, using the standard linearization technique, we obtain

[TABLE]

where

[TABLE]

Therefore, $\mathrm{var}(\hat{\tau}_{\mathrm{IPTW}})=\mathrm{var}(n^{-1}\sum_{i=1}^{m}\sum_{j=1}^{n_{i}}\tau_{ij})$, denoted as $V_{1}$.

To establish the asymptotic normality of $\hat{\tau}_{\mathrm{IPTW}}$, we use the central limit theory for dependent variables (Hoeffding et al., 1948; Serfling, 1968). Let $\mathrm{var}(\tau_{ij})=\sigma_{\tau}^{2}$ and $\mathrm{cov}(\tau_{ij},\tau_{ik})=\nu_{\tau}$ for $j\neq k$. Arrange the $\tau_{ij}$’s in a $n$-length sequence $\{\tau_{11},\ldots,\tau_{1n_{1}},\tau_{21},\ldots,\tau_{mn_{m}}\}$. To simplify the notation, let the $k$th random variable in this sequence be denoted by $\tau_{k}$, for $k=1,\ldots,n$. We now consider such sequences $\{\tau_{k}:k=1,\ldots,n\}$ are indexed by $n$. By Condition A1 (iv), the absolute central moments $E|\tau_{k}-E(\tau_{k})|^{3}$ are bounded uniformly in $k$. Moreover, by the assumption of $\sup_{1\leq i\leq m}n_{i}=O(n^{1/2})$, we then have $\mathrm{var}(\sum_{k=a+1}^{a+n}\tau_{k})\sim nA^{2}$, uniformly in $a$, as $n\rightarrow\infty,$ where $A^{2}$ is a positive constant. Following Serfling (1968), these are typical criterion for verifying the Lindeberg condition (Loève, 1960), and therefore $V_{1}^{-1/2}(\hat{\tau}_{\mathrm{IPTW}}-\tau)\rightarrow\mathcal{N}(0,1)$, in distribution, as $n\rightarrow\infty$.