An approximate Bayesian inference on propensity score estimation under unit nonresponse

TL;DR

This paper introduces an approximate Bayesian approach for propensity score estimation under unit nonresponse, providing a Bayesian inference method that aligns with frequentist confidence intervals without requiring Taylor linearization.

Contribution

It proposes a novel Bayesian method for nonresponse adjustment that avoids Taylor linearization and can handle over-identified and nonignorable nonresponse cases.

Findings

Simulation studies confirm the method's validity.

The approach aligns Bayesian credible regions with frequentist confidence intervals.

Application to survey data demonstrates practical utility.

Abstract

Nonresponse weighting adjustment using the response propensity score is a popular tool for handling unit nonresponse. Statistical inference after the nonresponse weighting adjustment is complicated because the effect of estimating the propensity model parameter needs to be incorporated. In this paper, we propose an approximate Bayesian approach to handle unit nonresponse with parametric model assumptions on the response probability, but without model assumptions for the outcome variable. The proposed Bayesian method is calibrated to the frequentist inference in that the credible region obtained from the posterior distribution asymptotically matches to the frequentist confidence interval obtained from the Taylor linearization method. Unlike the frequentist approach, however, the proposed method does not involve Taylor linearization. The proposed method can be extended to handle over-identified cases in which there are more estimating equations than the parameters. Besides, the proposed method can also be modified to handle nonignorable nonresponse. Results from two simulation studies confirm the validity of the proposed methods, which are then applied to data from a Korean longitudinal survey.

Tables

Table 1. Table 1: Simulation results: “m” denotes mean function, “c_p” is the coverage probability for the corresponding confidence interval, “CI length” is the average length of the confidence intervals.

Response	m	method	c_p	CI	Response	m	method	c_p	CI
mechanism				length	mechanism				length
R1	$m_1$	PS	0.95	1.83	R3	$m_1$	PS	0.95	1.86
		BPS	0.95	1.84			BPS	0.95	1.87
		OPS	0.95	1.78			OPS	0.95	1.78
		OBPS	0.95	1.78			OBPS	0.94	1.78
	$m_2$	PS	0.94	0.88		$m_2$	PS	0.94	0.89
		BPS	0.94	0.88			BPS	0.94	0.89
		OPS	0.94	0.79			OPS	0.93	0.79
		OBPS	0.94	0.80			OBPS	0.94	0.80
	$m_3$	PS	0.95	1.56		$m_3$	PS	0.94	1.58
		BPS	0.94	1.56			BPS	0.94	1.58
		OPS	0.95	1.53			OPS	0.95	1.53
		OBPS	0.94	1.52			OBPS	0.94	1.52
R2	$m_1$	PS	0.95	1.96	R4	$m_1$	PS	0.95	1.95
		BPS	0.96	2.00			BPS	0.95	1.99
		OPS	0.95	1.83			OPS	0.94	1.82
		OBPS	0.95	1.83			OBPS	0.95	1.82
	$m_2$	PS	0.95	1.16		$m_2$	PS	0.94	1.13
		BPS	0.95	1.16			BPS	0.95	1.13
		OPS	0.94	0.91			OPS	0.93	0.90
		OBPS	0.95	0.97			OBPS	0.95	0.95
	$m_3$	PS	0.95	1.68		$m_3$	PS	0.95	1.67
		BPS	0.95	1.72			BPS	0.95	1.70
		OPS	0.95	1.59			OPS	0.95	1.58
		OBPS	0.95	1.58			OBPS	0.95	1.57

Table 2. Table 2: Simulation results: “m” denotes mean function, “bias” is the estimator subtracting true value, “R_std” is the relative standard error which is relative to the standard error of full sample estimator.

Response	m	method	bias	R_std	Response	m	method	bias	R_std
${ℛ}_1$	$m_1$	CC	-0.16	1.16	${ℛ}_2$	$m_1$	CC	-0.14	1.17
		KC	-0.00	1.10			KC	-0.00	1.09
		FI	-0.00	1.09			FI	-0.00	1.08
		BDA	0.00	1.10			BDA	-0.00	1.09
	$m_2$	CC	-0.18	1.14		$m_2$	CC	-0.14	1.12
		KC	0.00	1.11			KC	-0.00	1.09
		FI	-0.01	1.10			FI	-0.01	1.08
		BDA	-0.00	1.11			BDA	-0.00	1.09
	$m_3$	CC	-1.13	1.15		$m_3$	CC	-0.95	1.13
		KC	-0.00	1.03			KC	0.01	1.02
		FI	-0.01	1.04			FI	0.01	1.03
		BDA	-0.00	1.04			BDA	0.01	1.03

Table 3. Table 3: The coverage probabilities for the proposed method

method	m	res	cp
BDA	m1	${ℛ}_1$	0.95
BDA	m2	${ℛ}_1$	0.94
BDA	m3	${ℛ}_1$	0.95
BDA	m1	${ℛ}_2$	0.95
BDA	m2	${ℛ}_2$	0.94
BDA	m3	${ℛ}_2$	0.95

Equations

$$\frac{1}{n}\sum _{i=1}^{n}U\left(\theta ;{\mathbf{x}}_i,y_i\right)=0$$

$$\displaystyle\delta_{i}=\left\{\begin{array}[]{ll}1&\text{if $y_{i}$ is observed}\\ 0&\text{otherwise.}\end{array}\right.$$

$$Pr({\delta }_i=1\mid {\mathbf{x}}_i,y_i)=\pi \left(\varphi ;{\mathbf{x}}_i,y_i\right)$$

$$\frac{1}{n}\sum _{i=1}^n\frac{{\delta }_i}{\pi (\hat{\varphi };{\mathbf{x}}_i,y_i)}U\left(\theta ;{\mathbf{x}}_i,y_i\right)=0$$

$$Pr\left(\delta =1|\mathbf{x},y\right)=Pr\left(\delta =1|\mathbf{x}\right).$$

$$L(\varphi )=\prod _{i=1}^n{\left\{\pi (\varphi ;{\mathbf{x}}_i)\right\}}^{{\delta }_i}{\left\{1-\pi (\varphi ;{\mathbf{x}}_i)\right\}}^{1-{\delta }_i}.$$

$$l_{obs}(\varphi )=\sum _{i=1}^n{\delta }_i\log \pi (\varphi ;{\mathbf{x}}_i,y_i)+\sum _{i=1}^n(1-{\delta }_i)\log \int \{1-\pi (\varphi ;{\mathbf{x}}_i,y_i)\}{\hat{f}}_1(y\mid {\mathbf{x}}_i)dy,$$

$$\sqrt{n}({\hat{\theta }}_{PS}-\theta )\stackrel{\mathcal{L}}{⟶}N(0,{\sigma }^2)$$

$$U_1\left(\varphi \right)=\frac{1}{n}\sum _{i=1}^n\left\{\frac{{\delta }_i}{\pi \left(\varphi ;{\mathbf{x}}_i\right)}-\frac{1-{\delta }_i}{1-\pi \left(\varphi ;{\mathbf{x}}_i\right)}\right\}\frac{\partial \pi \left(\varphi ;{\mathbf{x}}_i\right)}{\partial \varphi }=:\frac{1}{n}\sum _{i=1}^{n}s\left(\varphi ;{\mathbf{x}}_i,{\delta }_i\right).$$

$$U_2\left(\varphi ,\theta \right)=\frac{1}{n}\sum _{i=1}^n\frac{{\delta }_i}{\pi \left(\varphi ;{\mathbf{x}}_i\right)}U\left(\theta ;{\mathbf{x}}_i,y_i\right),$$

$$U_n(\mathbf{\zeta })=\left(\begin{array}{c}{U}_1\left(\varphi \right)\\ U_2\left(\varphi ,\theta \right)\end{array}\right).$$

$$p(\mathbf{\zeta }\mid \hat{\mathbf{\zeta }})=\frac{g(\hat{\mathbf{\zeta }}\mid \mathbf{\zeta })\pi (\mathbf{\zeta })}{\int g(\hat{\mathbf{\zeta }}\mid \mathbf{\zeta })\pi (\mathbf{\zeta })d\mathbf{\zeta }}$$

$$p(\mathbf{\zeta }\mid U_n)=\frac{g\{U_n(\mathbf{\zeta })\mid \mathbf{\zeta }\}\pi (\mathbf{\zeta })}{\int g\{U_n(\mathbf{\zeta })\mid \mathbf{\zeta }\}\pi (\mathbf{\zeta })d\mathbf{\zeta }},$$

$$\left[\sqrt{n}{U}_n|\mathbf{\zeta }\right]=\left[\sqrt{n}{U}_n|\mathbf{\eta }\right]\stackrel{\mathcal{L}}{\to }N\left(\mathbf{\eta },\Sigma \right),$$

$$\left[\mathbf{\eta }|U_n\right]\propto \left[U_n|\mathbf{\eta }\right]\pi \left(\mathbf{\eta }\right).$$

$$\left[\mathbf{\eta }\mid U_n\right]\sim N\left(0,\Sigma /n\right)$$

$$\hat{\Sigma }=\left(\begin{array}{cc}{n}^{-1}{\sum }_{i=1}^{n}s(\hat{\varphi };{\mathbf{x}}_i{)}^{\otimes 2}& n^{-1}{\sum }_{i=1}^n{\delta }_i{\hat{\pi }}_i^{-1}s(\hat{\varphi };{\mathbf{x}}_i)U^{\prime }(\hat{\theta };{\mathbf{x}}_i,y_i)\\ \text{mbox}symm.& n^{-1}{\sum }_{i=1}^n{\delta }_i{\hat{\pi }}_i^{-2}U(\hat{\theta };{\mathbf{x}}_i,y_i{)}^{\otimes 2}\end{array}\right),$$

$$p(\mathbf{\eta }\mid U_n=0)\stackrel{\mathcal{L}}{⟶}N(0,\hat{\Sigma }/n),$$

$$\sqrt{n}\left(U_n\left(\mathbf{\zeta }\right)-\mathbf{\eta }\left(\mathbf{\zeta }\right)\right)\stackrel{\mathcal{L}}{\to }N\left(0,\Sigma \left(\mathbf{\zeta }\right)\right)$$

$$\sqrt{n}\left(\hat{\mathbf{\zeta }}-{\mathbf{\zeta }}_0\right)\stackrel{\mathcal{L}}{\to }N\left(0,A({\mathbf{\zeta }}_0{)}^{-1}\Sigma \left({\mathbf{\zeta }}_0\right)A^{\prime }({\mathbf{\zeta }}_0{)}^{-1}\right)$$

$$p(\mathbf{\zeta }|\hat{\mathbf{\zeta }})\to {\varphi }_{\hat{\mathbf{\zeta }},Var(\hat{\mathbf{\zeta }})}\left(\mathbf{\zeta }\right)$$

$$p\underset{n\to \infty }{\lim }{\int }_{N_n\left({\mathbf{\zeta }}_0\right)}{\varphi }_{\hat{\mathbf{\zeta }},Var(\hat{\mathbf{\zeta }})}\left(\mathbf{\zeta }\right)d\mathbf{\zeta }=1,$$

$$Pr\{\mathbf{\zeta }:p(\mathbf{\zeta }\mid \hat{\mathbf{\zeta }})\ge k\}=1-\alpha .$$

$$C^{\ast }(\alpha )=\left\{\mathbf{\zeta }:p(\mathbf{\zeta }\mid \hat{\mathbf{\zeta }})\ge k^{\ast }(\alpha )\right\}.$$

$$\sum _{i=1}^n\frac{{\delta }_i}{\pi (\hat{\varphi };{\mathbf{x}}_i)}({\mathbf{x}}_i-{\mathbf{\mu }}_x)=0$$

$$U_3\left(\varphi ,{\mathbf{\mu }}_x\right)$$

$$U_4\left({\mathbf{\mu }}_x\right)$$

$$({\mathbf{\mu }}_x,\varphi ,\theta )\to ({\mathbf{\eta }}_1,{\eta }_2,{\mathbf{\eta }}_3,{\mathbf{\eta }}_4)$$

$$\left(\begin{array}{c}{\mathbf{\eta }}_1\\ {\mathbf{\eta }}_2\\ {\mathbf{\eta }}_3\\ {\mathbf{\eta }}_4\end{array}\right)=E\left\{\left(\begin{array}{c}{U}_1(\varphi )\\ U_2(\theta ,\varphi )\\ U_3\left(\varphi ,{\mathbf{\mu }}_x\right)\\ U_4\left({\mathbf{\mu }}_x\right)\end{array}\right)|{\mathbf{\mu }}_x,\varphi ,\theta \right\}.$$

$$U_n(\mathbf{\psi })={\left(U_1^{\prime }(\varphi ),U_2(\varphi ,\theta ),U_3^{\prime }(\varphi ,{\mathbf{\mu }}_x),U_4^{\prime }({\mathbf{\mu }}_x)\right)}^{\prime }.$$

$$[U_n|\mathbf{\psi }]\sim N(0,\Sigma (\mathbf{\psi })/n)$$

Full text

Key words: Approximate Bayesian computation, Posterior distribution, Missing at random, Nonignorable nonresponse, Nonresponse weighting adjustment.

1 Introduction

Missing data is frequently encountered in many areas of statistics. When the response mechanism is missing at random in the sense of Rubin (1976), one of the popular methods of handling missing data is to build a model for the response probability and use the inverse of the estimated response probability to construct weights for estimating parameters. Such weighting method is often called propensity score weighting and the resulting estimator is called propensity score estimator (Rosenbaum and Rubin, 1983). The propensity score method has been well established in the literature. For examples, see Rosenbaum (1987), Flanders and Greenland (1991), Robins et al. (1994), Robins et al. (1995), and Kim and Kim (2007). However, all the above researches were developed via the frequentist approaches. Variance estimates using a Taylor linearization method or bootstrap are used for making frequentist inference.

In this paper, we are interested in developing Bayesian inference for propensity score estimation. One of the main advantages of Bayesian inference is that all the uncertainty in the estimation process can be built into the Bayesian computation automatically. That is, there is no need to conduct variance estimation separately in the Bayesian inference. While the Bayesian method is widely used in many areas of statistics, the literature on the Bayesian approach of propensity score estimation is sparse. An (2010) proposed a Bayesian propensity score estimator jointly modeling the response mechanism and the outcome variable. However, specifying a correct outcome model is difficult under missing data and incorrect specification may lead to biased inference. McCandless et al. (2009) and Kaplan and Chen (2012) also assumed joint models and obtained Bayesian credible regions in the context of casual inference.

In this paper, our interest is in developing a new Bayesian approach without making any model assumptions on the outcome variable. Since no parametric model assumptions on the outcome variable are used, there is no explicit likelihood function corresponding to $\theta$, the main parameter of interest. This makes it difficult to develop a Bayesian method for propensity score estimation. The challenge thus lies in properly incorporating the uncertainty in the propensity score estimation process into the Bayesian framework.

In this paper, we propose a novel approach featuring approximate Bayesian computation based on the summary statistics (Beaumont et al., 2002). The sampling distribution of summary statistics, which is the estimating equation itself, can be used to replace the likelihood part in deriving the posterior distribution. In the proposed Bayesian method, the credible region obtained from the posterior distribution with a flat prior asymptotically matches the frequentist confidence interval obtained from the Taylor linearization method. The computation for the proposed method is relatively simple and easy to understand.

To guarantee the consistency of estimators, the propensity score method requires the correct specification of the response model. To protect against model misspecification, Robins et al. (1994), Scharfstein et al. (1999), and Bang and Robins (2005) proposed the so-called doubly robust estimation, which requires either the propensity score model or the outcome regression model be correctly specified. To achieve efficiency and robustness, we can add into the proposed Bayesian method additional estimating equations obtained from the auxiliary variables observed throughout the full sample. When we incorporate more equations than the parameters, the proposed Bayesian method is modified to solve the over-identifying situation.

The rest of the paper is organized as follows. In Section 2, we introduce the basic setup of the general propensity score estimation problem. The proposed method is presented in Section 3. The main result and asymptotic theory are discussed in Section 4. In Section 5, we developed a related method by extending our proposed method to incorporate the auxiliary information observed throughout the sample. We also presented how to incorporate data augmentation algorithm to handle nonignorable nonresponse in Section 6. The finite sample performance of the proposed methods is examined in an extensive simulation study in Section 7. An application of the proposed methods to a longitudinal survey is presented in Section 8. Some concluding remarks are made in Section 9.

2 Basic Setup

Suppose that we are interested in estimating $\theta$ defined through $E\left\{U\left(\theta;\bm{X},Y\right)\right\}=0$ for some estimating function $U(\theta;\bm{X},Y)$. Let $\left(\bm{x}_{i},y_{i}\right),i=1,\cdots,n,$ be independently and identically distributed (IID) realizations of random variable $\left(\bm{X},Y\right)$. Under complete data, we can obtain a consistent estimator of $\theta$ by solving

[TABLE]

for $\theta$. We assume that the solution to (1) is unique almost everywhere.

Now, suppose that $\bm{X}$ is always observed and $Y$ is subject to missingness. In this case, we can define the response indicator function for unit $i$ as

[TABLE]

We assume that $\delta_{i}$ are independently generated from a Bernoulli distribution with

[TABLE]

for some parameter vector $\phi$ and $\pi(\cdot)$ is a known function. In the logistic regression model, $\pi(x)=1/\{1+\exp(-x)\}$.

When missing data exist, we cannot apply (1) directly. Instead, using the parametric model for the response probability in (5), we can obtain the propensity score (PS) estimator of $\theta$ by the following two steps:

[Step 1] Compute the maximum likelihood (ML) estimator $\hat{\phi}$ of $\phi$.

[Step 2] Compute the PS estimator of $\theta$ by solving

[TABLE]

for $\theta$.

The computation for the ML estimator of $\phi$ can be greatly simplified if the response mechanism is Missing At Random (MAR) in the sense that

[TABLE]

In this case, the maximum likelihood estimator of $\phi$ can be obtained by finding the maximizer of

[TABLE]

If MAR does not hold, parameter estimation is more complicated. Assuming a parametric model for $f_{1}(y\mid\mathbf{x})=f(y\mid\mathbf{x},\delta=1)$, the ML estimator can be obtained by maximizing

[TABLE]

where $\hat{f}_{1}(y\mid\mathbf{x}_{i})$ is an estimator for $f_{1}(y\mid\mathbf{x}_{i})$. Riddles et al. (2016) proposed an alternative computational tool that avoids computing the above integration using an EM algorithm.

We shall first present our proposed method under the MAR assumption. An extension to Not Missing At Random (NMAR) will be discussed in Section 6. Once the PS estimator $\hat{\theta}_{PS}$ of $\theta$ is obtained from the above two-step procedure, statistical inference for $\theta$ can be made based on the asymptotic normality

[TABLE]

for some $\sigma^{2}>0$, where $\stackrel{{\scriptstyle\mathcal{L}}}{{\longrightarrow}}$ denotes convergence in distribution. See Chapter 5 of Kim and Shao (2013) for a justification for (7).

Under the above setup, we shall introduce the proposed Bayesian approach to estimate the parameter and make inference from the posterior distribution. An advantage of the Bayesian approach is that we can incorporate the uncertainty in estimating $\phi$ into the Bayesian computation automatically.

3 Proposed Method

We now present the proposed Bayesian method in the case of MAR. Under the parametric model assumption (5), the likelihood function for $\phi$ is given in (6). From the likelihood function, we can derive the score function for $\phi$ as

[TABLE]

If we define

[TABLE]

the PS estimator $\hat{\theta}_{PS}$ of $\theta$ can be viewed as the solution to the joint estimating equations: $U_{1}(\phi)=0$ and $U_{2}(\phi,\theta)=0$. Taylor linearization can be used to obtain a consistent variance estimator of $\hat{\theta}_{PS}$. See Chapter 5 of Kim and Shao (2013) for more details.

To introduce the proposed Bayesian inference corresponding to $\hat{\theta}_{PS}$, we first define $\bm{\zeta}=(\theta,\phi)$ and

[TABLE]

Instead of generating the posterior distribution from $p(\bm{\zeta}\mid\mbox{sample})$ directly , we use the posterior distribution $p(\bm{\zeta}\mid\hat{\bm{\zeta}})$ to approximate the posterior distribution $p(\bm{\zeta}\mid\mbox{sample})$, where $\hat{\bm{\zeta}}$ solves $U_{n}(\bm{\zeta})=0$. Thus, we can consider

[TABLE]

as an approximate posterior distribution for $\bm{\zeta}$, where $g(\hat{\bm{\zeta}}\mid\bm{\zeta})$ is the sampling distribution of $\hat{\bm{\zeta}}$ and $\pi(\bm{\zeta})$ is the prior distribution for $\bm{\zeta}$. However, finding the sampling distribution $g(\hat{\bm{\zeta}}\mid\bm{\zeta})$ will involve Taylor linearization.

To consider an alternative computation, instead of generating from $p(\bm{\zeta}\mid\hat{\bm{\zeta}})$ in (10), we use a posterior distribution from

[TABLE]

where $g\{U_{n}({\bm{\zeta}})\mid\bm{\zeta}\}$ is the sampling distribution of $U_{n}(\bm{\zeta})$. To generate samples from (11), we first make a transformation of the parameters, defined as $\bm{\eta}=E(U_{n}\mid\bm{\zeta})$. Thus, $T:\bm{\zeta}\rightarrow\bm{\eta}$ is an one-to-one transformation of the parameter. We will generate $\bm{\eta}^{*}$ from $p(\bm{\eta}\mid U_{n})$ first and then use $\bm{\zeta}^{*}=T^{-1}(\bm{\eta}^{*})$ to obtain the posterior distribution values from (11).

Now, to compute $p(\bm{\eta}\mid U_{n})$, first note that, under some regularity conditions,

[TABLE]

where notation $\left[\cdot\right]$ is used to denote the sampling distribution and $\stackrel{{\scriptstyle\mathcal{L}}}{{\longrightarrow}}$ denotes the convergence in distribution. Writing $\pi(\bm{\eta})$ as a prior distribution of $\bm{\eta}$, the posterior distribution of $\bm{\eta}$ given $U_{n}$ can be expressed as

[TABLE]

If there is no information for the prior, we can use a flat prior for $\bm{\eta}$. The sampling distribution $\left[U_{n}|\bm{\eta}\right]$ serves the role of the likelihood function in the Bayesian inference. Using (12) and a flat prior for $\bm{\eta}$, we obtain

[TABLE]

as the posterior distribution, where a consistent estimator of $\Sigma$ is

[TABLE]

where $\hat{\pi}_{i}=\pi(\hat{\phi};\mathbf{x}_{i})$, $\hat{\phi}$ and $\hat{\theta}$ solve $U_{n}\left(\phi,\theta\right)=\bm{0}$, $\bm{A}^{\otimes 2}=\bm{A}\bm{A}^{\prime}$ and $\bm{A}^{\prime}$ represents the transpose of $\bm{A}$. The details of the derivation are presented in Appendix A. After we obtain the posterior distribution of $\bm{\eta}$, we can use the inverse transformation of $T$ to obtain the posterior distribution of the original parameters. The following algorithm describes how to generate parameters from the posterior distribution of $\bm{\zeta}=(\phi,\theta)$:

[Step 1] Generate $\bm{\eta}^{*}$ from the posterior distribution

[TABLE]

where $\hat{\Sigma}$ is a consistent estimator of $Var(\sqrt{n}U_{n})=\Sigma$ in (12).

[Step 2] Solve $U_{n}\left(\bm{\zeta}\right)=\bm{\eta}^{*}$ for $\bm{\zeta}$ to obtain $\bm{\zeta}^{*}$.

Steps 1–2 can be repeated independently to generate independent samples from the posterior distribution. The samples can be used to obtain the posterior distribution of the induced parameters.

As we have illustrated before, the basic idea is that we use the posterior distribution of $p(\bm{\zeta}|\hat{\bm{\zeta}})$ to approximate the posterior distribution of $p\left(\bm{\zeta}|\text{sample}\right)$. This idea is similar in spirit to the Approximate Bayesian Computation of Soubeyrand and Haon-Lasportes (2015). Note that we do not use Taylor linearization to obtain the posterior distribution of $\theta$. Instead, we use a transformation technique and generate the posterior distribution of $p(\bm{\eta}|U_{n})$ first. After we obtain the posterior distribution of $\bm{\eta}$, we use the inverse transformation $T^{-1}:\bm{\eta}\rightarrow\bm{\zeta}$ to obtain the posterior distribution of the original parameters. The back-transformation plays the role of Taylor linearization in the frequentist approach. See Figure 1 for the illustration of the basic idea. Some asymptotic properties are established in the next section.

4 Asymptotic Properties

To establish the consistency of the parameter estimate and the interval estimate, we assume the following regularity conditions:

[C1] As $n\rightarrow\infty$, $U_{n}\left(\bm{\zeta}\right)\xrightarrow{}\bm{\eta}\left(\bm{\zeta}\right)$ in probability uniformly. That is $\sup_{\bm{\zeta}\in\bm{Z}}\|U_{n}\left(\bm{\zeta}\right)-\bm{\eta}\left(\bm{\zeta}\right)\|\xrightarrow{P}0$, where $\bm{Z}$ is the parameter space.

[C2] The mapping $\bm{\zeta}\mapsto U_{n}\left(\bm{\zeta}\right)$ is continuous and has exactly one zero $\hat{\bm{\zeta}}$ with probability one as $n\rightarrow\infty$.

[C3] Equation $\bm{\eta}\left(\bm{\zeta}\right)=0$ has exactly one root at $\bm{\zeta}=\bm{\zeta}_{0}$.

[C4] There exits a neighbor of $\bm{\zeta}_{0}$, denoted by $N_{n}\left(\bm{\zeta}_{0}\right)$, on which with probability one all $U_{n}\left(\bm{\zeta}\right)$ are continuously differentiable and the Jacobian $\partial U_{n}\left(\bm{\zeta}\right)/\partial\bm{\zeta}$ converge uniformly to a non-stochastic limit which is non-singular. Here, $N_{n}\left(\bm{\zeta}_{0}\right)$ is a ball with center $\bm{\zeta}_{0}$ and radius $r_{n}$, where $r_{n}$ satisfies $r_{n}\xrightarrow{}0$ and $r_{n}\sqrt{n}\xrightarrow{}\infty$. Also, we assume that $\partial^{2}U_{n,j}\left(\bm{\zeta}\right)/\partial\bm{\zeta}\partial\bm{\zeta}^{\prime}$ is finite for each entry for $j=1,2,\cdots,p$ and with probability one as $n\xrightarrow{}\infty$.

[C5] For any $\bm{\zeta}\in N_{n}\left(\bm{\zeta}_{0}\right)$,

[TABLE]

holds for some $\Sigma(\bm{\zeta})=Var\left\{\sqrt{n}U_{n}\left(\bm{\zeta}\right)|\bm{\zeta}\right\}>0$ that is independent of $n$.

As long as the samples satisfy some moment conditions, condition [C1] holds. Condition [C2] and [C3] are used to make sure that the solutions of estimating equation $U_{n}$ and estimating function $\bm{\eta}$ exist and are unique to avoid the model non-identifiability problem. The condition [C4] regulates the derivatives of the estimating equation to make sure that the variance converges. Condition [C5] provides the asymptotic distribution for the estimating equation. Under the above conditions, we can obtain

[TABLE]

where $A(\bm{\zeta})=\partial\bm{\eta}\left(\bm{\zeta}\right)/\partial\bm{\zeta}$.

We now make additional assumptions to establish the posterior consistency and convergence in distribution:

[C6] The prior distribution $\bm{\eta}\mapsto\pi\left(\bm{\eta}\right)$ is positive and Lipschitz continuous over the parameter space.

[C7] For any $\bm{\zeta}\in N_{n}\left(\bm{\zeta}_{0}\right)$, the variance estimator $\hat{\Sigma}\left(\bm{\zeta}\right)$ in (14) satisfies $\hat{\Sigma}\left(\bm{\zeta}\right)=\Sigma\left(\bm{\zeta}\right)\left\{1+o_{p}(1)\right\}$.

[C8] For any $\bm{\zeta}\in N_{n}\left(\bm{\zeta}_{0}\right)$, the mapping $\bm{\zeta}\mapsto\left|\Sigma\left(\bm{\zeta}\right)\right|^{-1}$ is Lipschitz continuous. Also, the mapping $\bm{\zeta}\mapsto x^{\prime}\left\{\Sigma\left(\bm{\zeta}\right)\right\}^{-1}x$ is Lipschitz continuous in the sense that there exists a constant $C\left(x\right)$ satisfying $\left\|x^{\prime}\left\{\Sigma\left(\bm{\zeta}_{1}\right)\right\}^{-1}x-x^{\prime}\left\{\Sigma\left(\bm{\zeta}_{2}\right)\right\}^{-1}x\right\|\leq C\left(x\right)\left\|\bm{\zeta}_{1}-\bm{\zeta}_{2}\right\|$, for any $\bm{\zeta}_{1},\bm{\zeta}_{2}\in N_{n}\left(\bm{\zeta}_{0}\right)$, for all $x\in\mathbb{R}^{p}$, where $p=\dim\left(\bm{Z}\right)$. And $C\left(x\right)$ is also Lipschitz continuous.

[C9] $\bm{\zeta}\mapsto U_{n}\left(\bm{\zeta}\right)$ and $\bm{\zeta}\mapsto\bm{\eta}\left(\bm{\zeta}\right)$ are one to one functions for any $\bm{\zeta}\in N_{n}\left(\bm{\zeta}_{0}\right)$. Also $\bm{\zeta}\mapsto\bm{\eta}\left(\bm{\zeta}\right)$ is Lipschitz continuous.

Condition [C6] is a common assumption for the prior and the flat prior satisfies this condition. Condition [C7] requires the variance estimator to be consistent. Conditions [C8] to [C9] are the sufficient conditions for the posterior distribution to be approximated by the proposed method. All the conditions can be easily satisfied if we assume variance estimate is continuous and has bounded eigenvalues.

Theorem 4.1

Let $\hat{\bm{\zeta}}$ be the solution to $U_{n}(\bm{\zeta})=0$. Under (C1)–(C9), the posterior distribution $p(\bm{\zeta}\mid U_{n}=0)=p(\bm{\zeta}|\hat{\bm{\zeta}})$, generated by the two-step method in Section 3, satisfies

[TABLE]

where $\phi_{\hat{\bm{\zeta}},Var(\hat{\bm{\zeta}})}\left(\cdot\right)$ is the density of the normal distribution with mean $\hat{\bm{\zeta}}$ and variance $Var(\hat{\bm{\zeta}})$.

The proof is shown in Appendix B. Result (17) is a convergence of the posterior distribution to normality and result (18) is the posterior consistency. By (17), the confidence region using the proposed Bayesian method is asymptotically equivalent to the frequentist confidence region based on asymptotic normality of $\hat{\bm{\zeta}}$. Thus, our proposed Bayesian method is calibrated to frequentist inference.

To construct a level $\alpha$ confidence region, let $k^{*}(\alpha)$ be the largest value of $k$ such that

[TABLE]

The level-$\alpha$ Bayesian High Posterior Density (HPD) confidence region (Chen and Shao, 1999) using $k^{*}$ is

[TABLE]

We can show that $\int_{\hat{C}^{*}\left(\alpha\right)}p(\bm{\zeta}\mid\hat{\bm{\zeta}})d\bm{\zeta}\xrightarrow{}1-\alpha$ in probability, where $\hat{C}^{*}(\alpha)$ is the confidence region from Monte Carlo samples, which are generated from the approximate target posterior distribution. See Hyndman (1996).

5 Optimal Estimation

We now extend the proposed method to incorporate additional information from the full sample. Note that the PS estimator applied to $\bm{\mu}_{x}=E(X)$ can be computed as the solution to

[TABLE]

which is not necessarily equal to $\hat{\mu}_{x,n}=n^{-1}\sum_{i=1}^{n}\bm{x}_{i}$. Including this extra information in the propensity score estimation, if done properly, will improve the efficiency of the resulting PS estimator. In the frequentist propensity score method, incorporating such extra information can be implemented by Generalized Method of Moments and it is sometimes called optimal PS estimation. See Cao et al. (2009), Zhou and Kim (2012) and Imai and Ratkovic (2014).

To include such extra information, we may add

[TABLE]

in addition to the original estimating equations based on $U_{1}(\phi)$ and $U_{2}(\phi,\theta)$ in (8) and (9), respectively. Note that we cannot directly apply the proposed two-step method in Section 3 in this case because there are more estimating equations than the parameters and the transformation

[TABLE]

is not one-to-one, where

[TABLE]

To solve this problem, instead of using the two-step method involving generation of $\bm{\eta}^{*}$ first from (14), we consider a direct sampling method that generates $\bm{\psi}^{*}=(\bm{\mu}_{x}^{*},\phi^{*},\theta^{*})$ from the posterior distribution of $\bm{\psi}=(\bm{\mu}_{x},\phi,\theta)$ given the observed data directly. To formally describe the procedure, first define

[TABLE]

Under some regularity conditions, we can obtain

[TABLE]

for sufficiently large $n$, where $\Sigma(\bm{\psi})=Var\left\{\sqrt{n}U_{n}(\bm{\psi})\mid\bm{\psi}\right\}$. Using (19) as the sampling distribution $g(U_{n}|\bm{\psi})$ of $U_{n}$ and using a prior $\pi(\bm{\psi})$ for $\bm{\psi}$, the posterior distribution of $\bm{\psi}$ can be written as

[TABLE]

Note that we can still use the approximate normality of $U_{n}$ to play the role of the likelihood function in the approximate Bayesian analysis. Note that even if the prior distribution is normal, the posterior distribution in (20) is no longer normal.

To obtain the posterior draws from (20), we can use a Monte Carlo method based on a version of Metropolis-Hastings algorithm (e.g. Chib and Greenberg (1995)). The computation details of the Monte Carlo method for generating samples from (20) are presented in Appendix C.

Note that, in generating samples from (20), the number of estimating equations is allowed to be greater than the number of parameters. Therefore, the proposed method is quite flexible in the sense that it can be applied to over-identified situations. Since the point estimator is asymptotically equivalent to the optimal PS estimator, the proposed method can thus be called optimal Bayesian PS (OBPS) method.

6 Nonignorable nonresponse

We now consider an application of the proposed Bayesian method to nonignorable nonresponse. Under the setup of Section 2, we first assume a parametric model for the response mechanism

[TABLE]

where $\pi\left(\cdot\right)$ is known up to $\phi$ and $\bm{x}_{i}=(\bm{x}_{i1},\bm{x}_{i2})$. The auxiliary variable $\bm{x}_{i2}$ is often called the response instrumental variable to avoid the non-identifiable problem in Wang et al. (2014). In addition, we assume a parametric model for the respondents’ outcome model

[TABLE]

for some $\gamma$. Using (21) and (22), we can obtain the following prediction model for the nonrespondents:

[TABLE]

where $O\left(\bm{x}_{1},y;\phi\right)=Pr\left(\delta=0|\bm{x}_{1},y\right)/Pr\left(\delta=1|\bm{x}_{1},y\right)$, $f\left(y|\bm{x},\delta=1\right)$. If $\pi(\phi;\mathbf{x}_{i1},y_{i})$ follows a logistic regression model such as $\pi(\phi;\bm{x}_{i1},y_{i})=\{1+\exp(x_{i1}\phi_{1}+y_{i}\phi_{2})\}^{-1}$ then $O\left(\bm{x}_{1},y;\phi\right)=\exp(-\phi_{2}y)$. See Kim and Yu (2011) for more discussion of the prediction model (23).

If $y_{i}$ were available throughout the sample, we could use

[TABLE]

as the estimating functions for $\bm{\zeta}=(\gamma,\phi,\theta)$, where $s_{1}(\gamma)$ is the score function of $\gamma$ with $s_{1}(\gamma;x_{i},y_{i})=\partial\log f\left(y_{i}|\bm{x}_{i},\delta_{i}=1;\gamma\right)/\partial\gamma$ and $S_{2}\left(\phi\right)$ is the score function of $\phi$. Writing the joint estimating equations as $U_{n}(\bm{\zeta})=\left(S^{\prime}_{1}(\gamma),S^{\prime}_{2}(\phi),U(\theta)\right)^{\prime}$ and $\bm{\eta}=E\left\{U_{n}(\bm{\zeta})\mid\bm{\zeta}\right\}$, the following two-step method can be used to generate the posterior samples of $\bm{\zeta}$.

[Step 1] Generate $\bm{\eta}^{*}$ from the approximate posterior distribution using $p(\bm{\eta}\mid U_{n}(\bm{\zeta}))$. Under a flat prior for $\bm{\eta}$, the posterior distribution of $\bm{\eta}$ can be obtained as a multivariate normal distribution with mean $\bm{0}$ and variance $\Sigma/n$. A consistent estimator of $\Sigma$ is

[TABLE]

where $\hat{\bm{\zeta}}=(\hat{\gamma},\hat{\phi},\hat{\theta})$ is the solution to $U_{n}(\bm{\zeta})=\bm{0}$ under complete response.

[Step 2] The posterior values of $\bm{\zeta}$ can be obtained by solving $U_{n}(\bm{\zeta})=\bm{\eta}^{*}$ for $\bm{\zeta}$.

Now, to implement the proposed Bayesian method under missing data, we can use Data Augmentation (DA) method of Tanner and Wong (1987). The DA algorithm consists of I-step and P-step. In I-step, the imputed values of $y_{i}$ are generated from the prediction model using the current parameter values. In P-step, the posterior values of the parameters are generated from the above two-step method using the current imputed data. To formally describe the proposed method, define $X_{n}=\left\{\bm{x}_{1},\cdots,\bm{x}_{n}\right\}$, $\bm{\delta}_{n}=\left\{\delta_{1},\cdots,\delta_{n}\right\}$ and $Y_{n}=(Y_{\text{obs}},Y_{\text{mis}})$, where $Y_{\text{obs}}$ and $Y_{\text{mis}}$ are the observed and missing part of $Y_{n}=(y_{1},\cdots,y_{n})$, respectively. The proposed DA algorithm can be described as follows:

I-step:

Given current parameter values $\bm{\zeta}^{*}$, generate imputed values $Y^{*}_{\text{mis}}$ from the prediction model (23) evaluated at the current parameter values.

P-step:

Using the current imputed data, apply the above two-step method of generating the parameter values $\bm{\zeta}^{*}$ from $p(\bm{\zeta}\mid U_{n}^{*}(\bm{\zeta}))$, where $U_{n}^{*}(\bm{\zeta})=U_{n}(\bm{\zeta};Y_{\text{obs}},Y^{*}_{\text{mis}})$.

The two steps are iteratively computed until some convergence criterion is satisfied. Once the posterior values of $\bm{\zeta}^{*}$ are obtained, the posterior values of $\theta^{*}$ can be used to perform Bayesian inference for $\theta$. To explain the proposed method further, denote $p_{U}(\bm{\zeta}\mid X_{n},Y_{n},\bm{\delta}_{n})=p(\bm{\zeta}\mid U_{n})$ to emphasize that $U_{n}(\bm{\zeta})$ is a function of $Y_{n}$. The I-step of the proposed method is to generate $Y_{\text{mis}}$ from the posterior predictive distribution of $Y_{\text{mis}}$ by

[TABLE]

where

[TABLE]

is generated from P-step. After convergence, the DA algorithm generates $\zeta$ from the posterior density

[TABLE]

7 Simulation Study

We perform two limited simulation studies to validate our theory and to check the robustness of our proposed methods. In the first simulation, the proposed method is evaluated under ignorable response mechanism. In the second simulation, the proposed method is applied to some nonignorable nonresponse mechasnism.

7.1 Simulation Study One

The first simulation study can be described as a $3\times 4$ factorial design, where the factors are outcome regression model for $E(y\mid\mathbf{x})$ and the response mechanism.

For the outcome regression models, we use $y=m(x_{1},x_{2})+e$ with three different mean functions given by

[TABLE]

where $c_{3}$ is chosen to give the same values for $E(y)$ in different mean functions. The explanatory variables $(x_{1},x_{2})^{T}$ are generated from $N(\bm{\mu},\Sigma_{x})$, with $\bm{\mu}_{x}=\left(2,8\right)^{T}$ and $\Sigma=\mbox{diag}\{4,8\}$. The error distribution is $e\sim N(0,\sqrt{|x_{1}|+1})$.

For the response mechanism, we use four different response mechanisms. In the first response mechanism (R1), the response indicator function $\delta_{i}$ are independently generated from a Bernoulli distribution with probability

[TABLE]

with $(\phi_{0},\phi_{1})=(0.1,0.4)$, which makes the overall response rate approximately equal to 70%. In the second response mechanism (R2), we use the sample logistic regression model with $(\phi_{0},\phi_{1})=(-1.2,0.15)$, which leads to about 30% response rate. In the third response mechanism (R3), the response indicator function $\delta_{i}$ are independently generated from a Bernoulli distribution with probability

[TABLE]

where $\Phi(\cdot)$ is the cumulative distribution function of the standard normal distribution and $(\phi_{0},\phi_{1})=(0,0.28)$, which leads to about 70% response rate. In the fourth response mechanism (R4), we use the same probit model with $(\phi_{0},\phi_{1})=(-0.7,0.1)$ to make the response rate near to 30%.

For each of the $12=3\times 4$ simulation setup, we generate random samples of size $n=500$ independently $B=2,000$ times. From each realized sample, we specify a logistic regression model

[TABLE]

as the response model. Thus, in R3 and R4, the response model is incorrectly specified.

For each Monte Carlo sample, we use the following four methods of inference for $\theta=E(y)$:

PS: Frequentist approach based on Taylor linearization. The point estimator $(\hat{\theta}_{PS},\hat{\phi})$ is computed from

[TABLE]

The confidence intervals are constructed by $\hat{\theta}_{PS}\pm 1.96\sqrt{\hat{V}_{PS}}$, where $\hat{V}_{PS}$ is obtained by the Taylor linearization method.

• 2.

Bayesian PS (BPS): Apply the proposed Bayesian method based on the joint estimating functions

[TABLE]

The estimators for $\phi,\theta$ are obtained by the median of the draws from the approximate posterior distribution. The confidence interval can be constructed by HPD region introduced in Section 4.

• 3.

Optimal PS (OPS): Generalized method of moments using

[TABLE]

in addition to (30) and (31). If we denote $U_{n}(\mu_{x},\phi,\theta)=(U^{\prime}_{1},U_{2},U^{\prime}_{3},U^{\prime}_{4})^{\prime}$, then the OPS estimator is obtained by minimizing $U_{n}^{T}W^{-1}U_{n}$, where $W=Var(U_{n})$. See Section 5.4 of Kim and Shao (2013).

• 4.

OBPS: Optimal Bayesian PS method discussed in Section 5 using the same estimating functions $U_{1}(\phi)$, $U_{2}(\phi,\theta)$, $U_{3}(\phi,\mu_{x})$, and $U_{4}(\mu_{x})$. The point estimators for $\bm{\mu}_{x},\phi,\theta$ are obtained by the median of the draws from the approximate posterior distribution. The confidence intervals can be constructed by the HPD region, introduced in Section 4.

For each of the four methods, 95% confidence intervals for $\theta$ are computed from Monte Carlo samples.

Table 1 presents the simulation results, coverage probabilities and average lengths of confidence intervals (CI), for the four methods. Overall, all the coverage probabilities are approximately 95%, which validates our proposed methods BPS and OBPS. For R1 and R2, we have a correctly specified model for the response mechanism. For R1, which has high response rate 70%, both BPS and OBPS methods provide valid confidence intervals with correct coverage rates. Comparing the average length of confidence intervals, we can see that PS and BPS methods have approximately equal average CI lengths and OPS and OBPS have approximately equal average CI lengths, which confirms the asymptotic equivalence of the two methods. That is, our proposed Bayesian methods are calibrated to the frequentist inference. The same conclusion can be obtained for R2, which has much lower response rates. For different regression mean functions, we find that both OPS and OBPS methods achieve more efficiency gains when the regression model is not linear and the response rate is low. For the probit response mechanism (R3 and R4), BPS and OBPS still provide valid confidence intervals with correct coverages. Thus, the proposed method seems to be robust against model misspecification of the response model.

7.2 Simulation Study Two

In the second simulation study, we consider an extension of the proposed method to nonignorable nonresponse. In the simulation, we generate the covariate variable $x\sim N(0,0.5)$ and use the outcome regression model $y=m(x)+e$ to generate $y$, where $e\sim N(0,1)$. We consider three different mean functions $m(x)$, which are specified as $m_{1}(x)=-1+2x$, $m_{2}(x)=-1.25+2x+0.5x^{2}$ and $m_{3}(x)=-1+8\sin(x)$.

We use two different mechanisms to generate the response indicators. The response indicator function $\delta_{i}$ are independently generated from Bernoulli distribution with the probability for $\delta_{i}=1$ equal to

[TABLE]

where $(\phi_{10},\phi_{11})=(0.8,-0.2),(\phi_{20},\phi_{21})=(0.5,-0.1)$ and $\Phi(\cdot)$ is cumulative distribution function of the standard normal distribution. The overall response rates are approximately around 70%. Thus, we have $3\times 2$ setup for the simulation study.

For each simulation setup, $n=500$ samples are generated independently for 2,000 times. For each Monte Carlo sample, we apply the following methods to estimate $\theta=E(y)$:

Full sample method: Use $\hat{\theta}=\sum_{i=1}^{n}y_{i}/n$, which is computed as a benchmark for the comparison.

• 2.

Complete-Case (CC) method: Estimate $\theta$ by removing nonresponse. That is, $\hat{\theta}_{CC}$ is obtained by solving $\sum_{i=1}^{n}\delta_{i}(y_{i}-\theta)=0$ for $\theta$.

• 3.

Kott and Chang (2010) (KC) method: Assume the response model is

[TABLE]

The KC estimates are obtained by solving

[TABLE]

• 4.

Fractional imputation (FI) method: Use $y|(x,\delta=1)\sim N(\beta_{0}+\beta_{1}x_{i},\sigma^{2})$ and the response mechanism in (35) to obtain the predictive model. The maximum likelihood estimator of $\theta$ is computed by using Fractional Imputation (FI) method in Kim (2011). Set the size of FI is 20. A description of the FI algorithm is described in Appendix D.

• 5.

Bayesian Data Augmentation (BDA) method: Apply the proposed method in Section 6 using the same model for FI method. In the data augmentation algorithm, we choose the burn-in size as 2,000 and after burn-in, iteration size is 2,000.

Thus, in the last two methods, the outcome model is misspecified under $m_{2}$ and $m_{3}$. Under $\mathcal{R}_{2}$, the response mechanism is slightly misspecified.

The simulation results are presented in Table 2 and 3. From Table 2, we can see that the performance of the proposed BDA method is similar to the KC and FI methods. Furthermore, the proposed BDA method can simultaneously construct correct confidence intervals and does not involve Taylor linearization. From Table 3, we can see that the coverage probabilities of the proposed method are around 0.95, which confirms the validity of the proposed BDA method.

8 Application

In this section, we apply the proposed Bayesian propensity score methods to Korea Labor and Income Panel Survey (KLIPS) data. A brief description of the panel survey can be found at http:// www.kli.re.kr/klips/en/about/introduce.jsp. The study variable (y) is the average monthly income for the current year and the auxiliary variable (x) can be demographic variables, such as the age groups and sex. Let $(X_{i},Y_{it})$ be the observations for household $i$ in panel year $t$. The KLIPS has $n=5,013$ households and $T=8$ panel years. We treat the first panel observations as the baseline measurements, and there are no missing data in the first year. In the panel survey, $X_{i}$ are completely observed and $Y_{it}$ are subject to missingness, for $i=1,2,\cdots,n$ and $t=1,2,\cdots,T$. Let $\delta_{it}$ be the response indicator function of $Y_{it}$. Define

[TABLE]

We are interested in estimating the probability of full response

[TABLE]

where $Y_{i,obs}=(Y_{i1},\cdots,Y_{iT})^{\prime}$ represents the observed responses for household $i$. The inverse of the $\pi_{i}$ in (39) can be used as the propensity weight for the penal survey. For monotone missing data, in the sense of $\delta_{it}=1$ implying $\delta_{i,t-1}=1,\cdots,\delta_{i1}=1$, the probability reduces to

[TABLE]

where $\pi_{it}=Pr(\delta_{it}=1|\delta_{i,t-1}=1,X_{i},Y_{i1},\cdots,Y_{i,t-1})$ under MAR assumption.

For arbitrary missing patterns as in KLIPS, we first define $\delta_{it}^{*}=\prod_{k=1}^{t}\delta_{ik}$. Note that $\delta_{it}^{*}=1$ implies that $\delta_{i,t-1}^{*}=1$. Furthermore,

[TABLE]

where $\pi_{i1}=1$ for all samples.

Thus, we can build a parametric model for $\pi_{it}=Pr(\delta_{it}=1|\delta_{i,t-1}^{*}=1,X_{i},Y_{i,t-1})$ and estimate the parameters sequentially. Instead of using the frequentist approach of Zhou and Kim (2012), we apply the BPS method in Section 3 and OBPS method in Section 5 to incorporate the extra information in $X$.

We are interested in estimating the average income for the final year and constructing confidence intervals for the parameters. Assume the response mechanism follows

[TABLE]

which is known up to parameter $\phi_{t}$. Thus, we allow that the response probability at year $t$ depends on the last year income $y_{t-1}$, but not on the current year income. Assume $\delta_{it}$, given $\delta_{i,t-1}^{*}=1,X_{i}$, and $Y_{i,t-1}$, independently follow Bernoulli distribution with probability $\pi(\phi_{t};X_{i},Y_{i,t-1})$ in (40). Therefore, the score function of $\phi_{t}$ is

[TABLE]

Then the joint estimating equations are $U_{n}(\phi_{2},\phi_{3},\cdots,\phi_{T},\theta)=0$, where

[TABLE]

and $\theta=E(Y_{T})$.

The Bayesian propensity score (BPS) method can be described as

Solve $U_{n}(\phi_{2},\phi_{3},\cdots,\phi_{T},\theta)=\mathbf{0}$ to obtain $\hat{\phi}_{2},\cdots,\hat{\phi}_{T}$, and $\hat{\theta}$.

• 2.

Generate $\eta^{*}=(\eta_{1}^{*^{\prime}},\eta_{2}^{*^{\prime}})^{\prime}$ from $N(\mathbf{0},\hat{\Sigma}/n)$, where $\hat{\Sigma}$ is a consistent variance estimator of $\sqrt{n}U_{n}(\phi_{2},\phi_{3},\cdots,\phi_{T},\theta)$.

• 3.

Solve $\left(S^{\prime}(\phi_{2}),\cdots,S^{\prime}(\phi_{T})\right)^{\prime}=\eta_{1}^{*}$ to obtain $\phi_{2}^{*},\cdots,\phi_{T}^{*}$.

• 4.

Compute $\pi_{i}^{*}=\pi(\phi_{2}^{*};X_{i},Y_{i,1})\times\cdots\times\pi(\phi_{T}^{*};X_{i},Y_{i,T-1})$. Solve

[TABLE]

to obtain $\theta^{*}$.

Repeat the above steps independently to generate samples from the posterior distribution of parameters. The variance-covariance matrix $\hat{\Sigma}$ can be derived by

[TABLE]

To improve the efficiency of the point estimator, we also apply OBPS method to the same data. In addition to equations in (45), we add

[TABLE]

where $\mu_{x}$ is the marginal proportion vector for demographical covariates. Therefore, the posterior distribution of $\theta$ can be obtained by applying the proposed algorithm in Section 5.

For a comparison, we also considered a naive method which does not use the propensity model and apply the Bayesian method in the complete cases (CC) only. We apply BPS, OBPS and CC method to $T=2,3,4$. The numerical results are presented below.

From Figure 2, all three methods provide similar estimators for the average income $\theta$.The trend of average income goes up as year $T$ increases. For year $T=2$, all three methods provide similar mean estimates. But the OBPS method is the most efficient. For year $T=3$, we see that the CC method provides lower mean estimate than BPS or OBPS, which is due to the nonresponse bias in the CC method. This phenomenon becomes more obvious for year $T=4$. Also, the lengths of confidence intervals increase as $T$ increases, since the fully observed sample size is decreasing due to panel attrition. The CC method presents smaller values of $\theta$ for $T=4$, which suggests more panel attrition for higher income households. Both BPS and OBPS provide similar mean estimates. But the OBPS method has narrower confidence intervals, which confirms the efficiency of the OBPS method.

9 Concluding Remarks

A new Bayesian inference using PS method is developed using the idea of Approximate Bayesian computation. The proposed method can be widely applicable due to popularity of PS method. The proposed Bayesian approach is calibrated to frequentist inference in the sense that the proposed method provides the same inferential results with its frequentist version asymptotically (Little, 2012). The calibration property holds if the prior distribution for the model parameters is flat. If the prior is informative then the resulting Bayesian inference will be more efficient than frequentist inference thanks to its natural incorporation of the prior information. Thus, the proposed method is applicable when the need of combining information from different sources.

Causal inference, including estimation of average treatment effect from observational studies, can be one promising application area of the PS method (Morgan and Winship, 2014 and Hudgens and Halloran, 2008). Developing tools for causal inference using the Bayesian PS method will be an important extension of this research. Also, Bayesian model selection method (Ishwaran and Rao, 2005) can be naturally applied to this setup. Such extensions will be topics for future research.

Appendix

A. Consistent variance estimator

From the asymptotic distribution in (12), we can write

[TABLE]

To emphasize that $\Sigma$ is a function of $\phi,\theta$, we use $\Sigma=:\Sigma\left(\phi,\theta\right)$. Since the transformation

[TABLE]

is one-to-one, $\Sigma\left(\phi,\theta\right)$ is equivalent to $\Sigma\left(\bm{\eta}_{1},\eta_{2}\right)$, where $\bm{\eta}=(\bm{\eta}^{\prime}_{1},\eta_{2})^{\prime}$. The corresponding density function is

[TABLE]

Since we have assigned a flat prior, in the sense of $\pi\left(\bm{\eta}\right)\propto 1$, we can derive the posterior distribution as

[TABLE]

By the definition of $\bm{\eta}$, $U_{n}$ is the unbiased estimator of $\bm{\eta}$. Thus, we write $\bm{\hat{\eta}}=U_{n}$. To show that

[TABLE]

we first show that $\Sigma(\cdot)$ is continuous, which can be proved by the dominated convergence theorem applied to $\bm{\eta}_{1}(\cdot)$ and $\eta_{2}(\cdot)$. Now, noting that, by asymptotic distribution (12) and Chebyshev’s inequality, we can show that $U_{n}\xrightarrow{P}\bm{\eta}$. Thus, we can obtain $\Sigma(\bm{\hat{\eta}})\xrightarrow{P}U_{n}(\bm{\eta})$. Since $\Sigma$ is positive definite and $x^{-1/2}$ is continuous if $x>0$, $|\Sigma(\bm{\hat{\eta}})|^{-1/2}\xrightarrow{P}|\Sigma(\bm{\eta})|^{-1/2}$. By the continuous mapping theorem,

[TABLE]

Therefore, we can derive the posterior distribution as

[TABLE]

That is $[\bm{\eta}|U_{n}=\mathbf{0}]\sim N(\mathbf{0},\Sigma(U_{n}=\mathbf{0})/n)$, which is equivalent to

[TABLE]

where $(\hat{\phi},\hat{\theta})$ is the solution to $U_{n}=\mathbf{0}$. Furthermore, the consistency of $\hat{\Sigma}=\hat{\Sigma}(\hat{\phi},\hat{\theta})$ in (3) can be proved using the law of large numbers.

B. Proof of Theorem 4.1

Step I

From Condition [C9], we assume that $\bm{\zeta}\mapsto U_{n}(\bm{\zeta})$ and $\bm{\zeta}\mapsto\bm{\eta}(\bm{\zeta})$ are one-to-one functions, for any $\bm{\zeta}\in N(\bm{\zeta}_{0})$. Denote these two mappings as $T_{n}$ and $T$ respectively. Because of their one-to-one property, their inverse mappings exist for $\bm{\zeta}\in N_{n}(\bm{\zeta}_{0})$. Therefore, we can write (15) as

[TABLE]

which leads to

[TABLE]

Thus, by the convergence of $U_{n}$ to $\bm{\eta}$ and using the argument similar to the proof for Lemma 1 in Soubeyrand and Haon-Lasportes (2015), we can show that

[TABLE]

Step II

Note that $U_{n}(\hat{\bm{\zeta}})=0$, thus $T_{n}^{-1}(0)=\hat{\bm{\zeta}}$. From (B.1), we can therefore get the posterior distribution

[TABLE]

Thus, we can write the density $p(\bm{\eta}|U_{n}=0)$ as

[TABLE]

Furthermore, by the consistency of the variance estimator provided in condition [C7], we can obtain $\hat{\Sigma}:=\hat{\Sigma}(\hat{\bm{\zeta}})=\Sigma(\hat{\bm{\zeta}})\left\{1+o_{p}(1)\right\}$. Thus,

[TABLE]

which leads to

[TABLE]

From [C1] and [C5], we have $U_{n}(\bm{\zeta})\xrightarrow{}\bm{\eta}$ in probability and $U_{n}=O_{p}(1/\sqrt{n})$ for $\bm{\zeta}\in N(\bm{\zeta}_{0})$, which leads to $\bm{\eta}=O(1/\sqrt{n})$. Thus,

[TABLE]

in probability and the following follows

[TABLE]

Step III

Let $\bm{\eta}^{*}$ be generated from the asymptotic posterior distribution (B.3) which is a normal distribution with mean 0 and variance $\hat{\Sigma}/n$. Therefore, the $j$-th component $\bm{\zeta}_{j}^{*}$ of $\bm{\zeta}^{*}$ satisfies

[TABLE]

By $E(\bm{Z}^{\prime}\Lambda\bm{Z})=tr(\Lambda\Sigma)+\mu^{\prime}\Lambda\mu$, we derive

[TABLE]

under [C4]. Therefore, we have

[TABLE]

for $j=1,2,\cdots,p$, which establishes

[TABLE]

Step IV

Now, the posterior variance of $\bm{\zeta}^{*}_{j}$:

[TABLE]

The first term is

[TABLE]

For the second term, using

[TABLE]

for $\bm{Z}\sim N(\mu,\Sigma)$, we have

[TABLE]

The covariance of two terms is less than the square root of their variances. We have shown that the variance of the first term is in the order of $O(1/n)$ and the variance of the second term is in the order of $O(1/n^{2})$. So the covariance has the order of $O(n^{-3/2})$.

Similarly, we can derive

[TABLE]

Combining (B.5) and (B.6), we have

[TABLE]

Step V

By Conditions [C1]-[C5], we have

[TABLE]

where $A(\bm{\zeta})=\partial\bm{\eta}(\bm{\zeta})/\partial\bm{\zeta}$. See Theorem 5.21 in Van der Vaart (2000).

Since $T_{n}\xrightarrow{}T$ uniformly by [C1] and both mappings are one-to-one functions, we can state that $T^{-1}_{n}\xrightarrow{}T^{-1}$ uniformly. Thus,

[TABLE]

Also, by [C7], we have ${\hat{\Sigma}}\xrightarrow{P}\Sigma(\hat{\bm{\zeta}})$ and $\hat{\bm{\zeta}}\xrightarrow{P}\bm{\zeta}_{0}$. By the Lipschitz continuity of $\Sigma(\zeta)$ , we can conclude that $\Sigma(\hat{\bm{\zeta}})\xrightarrow{P}\Sigma(\bm{\zeta}_{0})$. Thus, $\hat{\Sigma}\xrightarrow{P}\Sigma(\bm{\zeta}_{0})$

and

[TABLE]

by the continuous mapping theorem. Combining the previous conclusions (B.4), (B.7) and (B.8), we can use Slutsky’s theorem to get

[TABLE]

which proves (17).

Step VI

Let $\alpha\in(0,1)$, and define

[TABLE]

where the $\chi_{p}^{2}(\alpha)$ is the $\alpha$ quantile of Chi-square distribution with $p$ degrees of freedom.

Furthermore, from a property of the Raylei quotient (Horn and Johnson, 1985), there exists a matrix $O$ such that

[TABLE]

where $OO^{T}=\bm{I}_{p}$ and $0<\lambda_{1}\leq\lambda_{2},\cdots,\leq\lambda_{p}$. Thus we obtain

[TABLE]

Also, we can apply the conclusion (B.9) to get

[TABLE]

for all $\bm{\zeta}^{*}\in C_{n,\alpha}$.

Similarly, by the asymptotic normality of the estimator $\hat{\bm{\zeta}}$ and applying the conclusion (B.9),

[TABLE]

Next, from (B.10), we can conclude that

[TABLE]

By the inequality $\|\bm{\zeta}^{*}-\bm{\zeta}_{0}\|\leq\|\hat{\bm{\zeta}}-\bm{\zeta}^{*}\|+\|\hat{\bm{\zeta}}-\bm{\zeta}_{0}\|$, we obtain

[TABLE]

Since we have defined $N_{n}(\bm{\zeta}_{0})$ in a neighborhood with center $\bm{\zeta}_{0}$ and radius $r_{n}$, where $r_{n}$ satisfies $r_{n}\xrightarrow{}0$ and $\sqrt{n}r_{n}\xrightarrow{}\infty$. From (B.11),

[TABLE]

Therefore,

[TABLE]

This is equivalent to

[TABLE]

The above conclusion holds for any $\alpha\in(0,1)$. Thus

[TABLE]

C. Computational Details for the Metropolis-Hastings Algorithm

Implementing the optimal Bayesian propensity score (OBPS) method is done through the following algorithm.

Choose the initial value for $\bm{\psi}$ and denote it as $\bm{\psi}_{0}$.

• 2.

For iteration $t$, given the current parameter value $\bm{\psi}_{t}$, generate $\Delta\bm{\psi}$ from $N(\mathbf{0},V)$, where $V$ is a tunning parameter obtainable by the data-driven method discussed below. Let the candidate value be $\bm{\psi}^{*}=\bm{\psi}_{t}+\Delta\bm{\psi}$.

• 3.

Compute the acceptance probability

[TABLE]

• 4.

Generate $u$ from Uniform $(0,1)$ distribution. If $u<\alpha$, accept the candidate $\bm{\psi}_{t+1}=\bm{\psi}^{*}$. Otherwise let $\bm{\psi}_{t+1}=\bm{\psi}_{t}$.

• 5.

For burning in period, discard the values from the first $B$ iterations. Then collect $M$ values. These $M$ values can be treated as values generated from the target posterior distribution.

For the choice of the initial value for $\bm{\psi}$, we can use the solution to

[TABLE]

In Metropolis-Hastings algorithm, the value of $V$ for the random walk will directly affect the speed of convergence of the Markov chain and the acceptance rate. We recommend a data-driven method to set $V$. A data-driven choice of $V$ can be obtained from the posterior variance of the Monte Carlo samples from $p\{\bm{\mu}_{x},\phi,\theta|(U_{1},U_{2},U_{4})=\mathbf{0}\}$.

To compute the acceptance probability, we need to compute the ratio

[TABLE]

which can be approximated by

[TABLE]

where

[TABLE]

and $(\hat{\bm{\mu}}_{x},\hat{\phi},\hat{\theta})$ are the consistent estimators.

D. Fractional imputation algorithm in simulation study two

Let $s(\phi;\delta_{i},\mathbf{x}_{1},y)$ be the score function for the response model. In additional to the response model, we also assume $f(y\mid x,\delta=1;\gamma)$. To solve the observed score function of $\phi$, that is

[TABLE]

where the conditional expectation is with respect to the prediction model in (23). The estimate $\hat{\gamma}$ of $\gamma$ can be obtained by using the observed data to fit the model $f(y\mid x,\delta=1;\gamma)$. To compute the solution $\hat{\phi}$ to (D.1), EM algorithm using fractional imputation (Kim, 2011) can be used. The algorithm is described as followings:

• E-step:

For each unit $i$ with $\delta_{i}=0$, generate $y_{ij}^{*}$ from $f_{1}(y\mid x_{i},\delta_{i}=1;\hat{\gamma})$ for $j=1,2,\cdots,b$. Given the current value of $\phi_{1}$, compute the fractional weights of $y_{ij}$ as

[TABLE]

subject to $\sum_{j=1}^{b}w_{ij}^{*}=1$.

• M-step:

Update $\phi$ by solving

[TABLE]

Repeat E-step and M-step iteratively until convergence. After convergence, the final estimator of $\theta=E(Y)$ is constructed by

[TABLE]