Non-response in sample survey is a common source of non-sampling error that appears when part of the data to be collected is not observed. In the case of missing at random (MAR) in which non-response occurs randomly, many appropriate statistical methods have been developed. On the other hand, in the case of non-ignorable non-response, there are relatively few studies on this subject. Non-ignorable non-response is known to cause bias and so accurate bias estimation is the key to properly handle the non-response. The PSA estimator, propensity-score-adjusted estimator, defined by

is widely used to reduce non-response bias. Here R_i is 1 if response and R_i is 0 if non-response and s is an index set of sample. Also, w_i is the sample weight and p̂_i is the estimated response probability. Kim and Riddles (2012) developed some asymptotic theories of the PSA estimator and suggested minimum variance of the estimator. Kim and Yu (2011) studied a mean estimation method for non-ignorable non-response case in which paper the missing value of the variable of interest is non-parametrically calculated using a kernel estimator. Riddles et al. (2016) proposed a propensity score adjustment method for non-ignorable non-response and showed that the proposed method is more efficient than the calibration-weighting method in Chang and Kott (2008) and Kott and Chang (2010).

In order to use the PSA estimator, it is necessary to estimate the response probability. Estimating response probabilities relies heavily on the use of model. Iannacchione et al. (1991) used a logistic regression model for the response probability estimation and this logistic regression model is widely used. Also Da Silva and Opsomer (2006) and Da Silva and Opsomer (2009) considered nonparametric methods to obtain the response probability. Bethlehem (2020) presented an approximate bias of the sample mean for non-ignorable non-response. Various response probability estimation methods have been addressed in that paper.

In general, the final sample weight considering non-response is used for estimation defined by

Therefore, the PSA estimator using ${\widehat{w}}_i^F=w_i/{\widehat{p}}_i$ is one of the non-response adjusted weight estimators. Chung and Shin (2017) suggested a non-response adjusted weight estimation method using a post-stratification which is addressed in Bethlehem (2020). This method reduces the effect of bias and gives a good performance in terms of the MSE.

In this study, we propose a new response probability estimation method that combines the response probability estimation method using logistic regression model with the post-stratification method suggested by Chung and Shin (2017).

The composition of this paper is as follows. In Section 2, the existing methods and a proposed method for the response probability estimation are explained. Section 3 describes the bias corrected PSA estimators using the response probability estimates obtained in Section 2. Section 4 confirms the superiority of the proposed method through simulation studies. There is a conclusion in Section 5.

2.1. Modeling method

In order to properly handle non-ignorable non-response, the response probability of each data must be appropriately estimated. As mentioned in Bethlehem (2020), various methods can be used for the estimation. The propensity score using the logistic regression model is the most commonly used. Let R_i be the indicator variable representing the response and p_i be the response probability. Also let the auxiliary variable x_i be a p-dimensional vector and always observable. Usually for non-ignorable non-response, it is assumed that p_i is a function of y_i and x_i. However, in this study, we assume that the response probability p_i is a function of y_i for simplicity.

Then, we can write p_i = P(Rⁱ = 1|y_i, x_i) = P(R_i = 1|y_i) = g₂(y_i; φ). Practically, we do not have enough information to estimate p_i and mostly the function is not known. To handle this situation, Kim and Yu (2011) and Kim and Riddles (2012) used some follow-up samples. However, this is not the usual case.

Mostly, the variable of interest is related to the auxiliary variables and therefore many studies consider a super-population model. Let the super-population model be y_i = g₁(x_i; β) + ε_i and assume

where θ is a vector of some unknown parameters and η_i is an error related to ε_i. Since this study considers non-ignorable non-response, it may not be appropriate to use a function of the auxiliary variable x_i. However, in reality, in order to estimate the response probability, we have no choice but to use the available auxiliary variables. Considering this, equation (2.1) may be useful. However, bias may occur in the result using the response probability obtained in equation (2.1).

The approximate estimate of p_i can be obtained using the available auxiliary variable x_i. The most common response probability estimation method is to use the parametric logistic regression model. Of course, nonparametric methods can be used as Da Silva and Opsomer (2006) and Da Silva and Opsomer (2009). Da Silva and Opsomer (2009) showed that local polynomial regression results in better practical and theoretical properties for the ignorable case. Recently Bethlehem (2020) studied the non-ignorable non-response and used the logistic regression modeling. Since the logistic regression model is commonly used and simple to implement, we use the logistic regression model in this paper. The logistic regression model is defined by

Using this model, we can obtain the response probability estimate ${\widehat{p}}_i^P$. Now let ${\widehat{w}}_i^P=w_i/{\widehat{p}}_i^P$. Then the final adjusted sample weight is obtained by following,

In Riddles et al. (2016), (2.2) is called naive estimator for non-ignorable non-response.

2.2. Post-stratification method

Several methods estimating response probability have been developed. One of them is the post-stratification method. As mentioned in Bethlehem (2020), post-stratification is a well-known and frequently used weighting technique. Usually categorical variables are used. Using these variables, population is divided into a number of non-overlapping subpopulations, called strata. Of course, continuous variables can be used and Bethlehem (2020) used the estimated response probability to construct subpopulations. Although (2.1) is an approximate result, naturally we can use the available auxiliary variable in order to construct strata. Chung and Shin (2017) and Min and Shin (2018) proposed a response probability estimation method using strata. In the proposed method by Min and Shin (2018), population is divided into L strata with boundaries obtained using percentiles of given auxiliary variable.

To construct the strata, we determine the total number of strata, L satisfying about r/L ≥ 10 where r is the total number of response samples. Then with L, we calculate the percentiles of the auxiliary variables and use the percentiles as boundaries. For instance, with L = 10, we use 10, 20, . . . , 90 percentiles for boundaries. For univariate case, the percentiles become directly boundaries. However, for multivariate case, we use the Cartesian product with the percentiles of each auxiliary variable to meet the total number of strata L.

Now, let N_h and r_h be the number of population and the number of final response samples in the h^th stratum, respectively. Then the sample weight of y_i in the h^th stratum using the post-stratification method is defined by

where s_h is the index set of h^th stratum sample. Since this method does not use any model to estimate the response probability, it is model free and naturally robust of model misspecification.

2.3. Proposed response probability estimation method

The widely used ${\widehat{w}}_i^{F(P)}$ is obtained using the individual data values. However, since it is based on a model, the results may vary depending on the model used. On the other hand, ${\widehat{w}}_i^{F(D)}$ has the advantage of being easy to calculate because it uses the number of samples in a given stratum rather than the individual data values. Again, the result may depend on the method of constructing the strata and the optimal boundaries of the strata are determined according to the given data.

A new method combining two methods explained in section 2.1 and 2.2 is proposed. First we calculate ${\widehat{w}}_i^{F(P)}$ explained in section 2.1. Then we make sum of ${\widehat{w}}_j^{F(P)}$ in the h^th stratum be N_h. That is, the final sample weight is defined by

where N_h and s_h are the same as those defined in (2.3).

In the case of non-ignorable non-response, since the response probability is a function of the variable of interest, it is necessary to estimate the response probability using the variable of interest. However, as mentioned before it is practically hard to obtain the response probability using the variable of interest. In Section 2, the response probability is estimated using the available auxiliary variables and this usually produces bias in estimation. To obtain better estimation results, the bias should be estimated and corrected.

3.1. Bias estimation

For bias estimation, this section considers a super-population model. Let the inclusion probability, π_i be a function of the variable of interest y_i and f_p(y_i|x_i) be a population distribution with a super-population model. Pfeffermann et al. (1998) showed f_s(y_i|θ^*, x_i) = f (y_i|i ∈ s, x_i) = {Pr(i ∈ s|y_i, x_i) f_p(y_i|θ, x_i)}/{Pr(i ∈ s|x_i)} where θ^* is a function of θ. Since Pr(i ∈ s|y_i, x_i) = E_p(π_i|y_i, x_i) and Pr(i ∈ s|x_i) = E_p(π_i|x_i), finally we have

where f_s(y_i|x_i) is a sample distribution, E_p(π_i|y_i, x_i) is the sample inclusion probability given x_i, y_i and E_p(π_i|x_i) is the sample inclusion probability given x_i. Whenever E_p(π_i|y_i, x_i) = E_p(π_i|x_i), then population distribution and sample distribution are the same. Therefore, with known inclusion probability and population distribution we can obtain the sample distribution and finally, we can calculate the bias. Bias estimation using various inclusion probabilities and population distributions have been obtained in previous studies.

In this study we consider a non-ignorable non-response with a non-informative sampling. Actually the final inclusion probability π_i is obtained by combining the inclusion probability at the time of sampling design and the response probability. As a non-informative sampling design, we use simple random sampling and so w_i in (1.1) is constant. Therefore π_i used in (3.1) is defined by π_i = p_i/w_i = p_i/w where w = N/n. Since p_i is a function of the variable of interest y_i, we have f_s(y_i|x_i) ≠ f_p(y_i|x_i) and so we can apply the bias estimation method developed in the informative sampling design to non-ignorable non-response.

We consider a linear inclusion probability model whose bias is easily estimated theoretically. This is because, if the linear inclusion probability model is effective, it can be expected to be effective when using other inclusion probability models. The linear inclusion probability model considered is as follows,

Then simply we have

Let $E_s\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})={\mu }_i^{(s)}$. Then plugging (3.4) into (3.1) we have,

where μ_i = E_p(y_i|x_i). Since $E_p\hspace{0.17em}(y_i^2\mid {\mathbf{x}}_{\mathbf{i}})={\text{Var}}_p\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})+{\mu }_i^2$, we have ${\mu }_i^{(s)}={\mu }_i+b_1/b_0+b_1{\mu }_i\times {\text{Var}}_p\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})$ and finally the bias is obtained as following,

3.2. Parameter estimation for inclusion probability model

Bias is estimated based on the super-population model and the inclusion probability model established.

In (3.2), E_p(π_i|y_i, x_i) can be calculated using known parameters b₀, b₁ and given y_i for the informative sampling. However, since we consider the non-ignorable non-response situation with non-informative sampling scheme, this situation usually can not happen in practice. The parameters b₀, b₁ should be estimated from the data.

In this study we use the linear probability model (3.2) defined by E_p(π_i|y_i, x_i) = b₀ + b₁y_i. Since E_p(π_i|y_i, x_i) = 1/E_s(w_i|y_i, x_i) and E_s(w_i|y_i, x_i) ≈ w_i from Pffeffermann and Sverchkov (2003) we have 1/w_i ≈ b₀ + b₁y_i. Therefore by plugging the estimated final sample weight ${\widehat{w}}_i^F$ obtained in Section 2.1 into w_i, we have the following model.

So using simple regression analysis we have the estimates of b₀, b₁.

3.3. Parameter estimation for super-population model

Three super-population models are considered in this study.

• Normal distribution

  (3.8)$$f_p\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})=N\hspace{0.17em}\left({\beta }_0+{\beta }_1^T{\mathbf{x}}_{\mathbf{i}},{\sigma }^2\right).$$

  That is, $y_i={\beta }_0+{\beta }_1^T{\mathbf{x}}_{\mathbf{i}}+{\varepsilon }_i$ where ${\beta }_1^T$ is a p-dimensional vector of parameters, ${\varepsilon }_i\stackrel{iid}{\sim }N(0,{\sigma }^2)$ and Var_p(y_i|x_i) = σ².

• Gamma distribution with log-linear model

  (3.9)$$f_p\hspace{0.17em}(y_i\mid {\alpha }_i,{\mu }_i)\propto y_i^{\alpha -1}\hspace{0.17em}\text{exp\hspace{0.17em}}\left(-\alpha \frac{y_i}{{\mu }_i}\right),$$

  where μ_i = E_p(y_i|x_i) and $\text{ln}({\mu }_i)={\beta }_0+{\beta }_1^T{\mathbf{x}}_{\mathbf{i}}$. That is, we use Γ(α, β^*) with β^* = μ_i/α and so we have ${\text{Var}}_p\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})={\mu }_i^2/\alpha$.

• Log-normal distribution with log-linear model

  (3.10)$$f_p\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})=\text{LN}\hspace{0.17em}\left({\beta }_0+{\beta }_1{\mathbf{x}}_{\mathbf{i}},{\sigma }^2\right),$$

  where ${\mu }_i=E_p\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})=\text{exp}({\beta }_0+{\beta }_1^T{\mathbf{x}}_{\mathbf{i}}+{\sigma }^2/2)$ and ${\text{Var}}_p\hspace{0.17em}(y_i\mid {\mathbf{x}}_{\mathbf{i}})={\mu }_i^2(\text{exp}({\sigma }^2)-1)$.

As used in Riddles et al. (2016), using the obtained data (x_i, y_i) and the known super-population model, we can estimate μ_i and Var_p(y_i|x_i).

3.4. Bias corrected propensity-score-adjusted (PSA) estimator

Bias corrected PSA estimator can be obtained based on the bias estimates and the PSA estimator defined in (1.2). That is, the bias corrected PSA estimator is defined by

In (3.11) we can use the final sample weights, ${\widehat{w}}_i^{F(P)},{\widehat{w}}_i^{F(D)}$ and ${\widehat{w}}_i^{F(C)}$ defined in Section 2 and ${\widehat{\text{bias}}}_i$ estimated according to the super-population model. Therefore, we have three bias corrected PSA estimators depending on the final sample weight,

To investigate the finite sample properties of the proposed method, we perform simulation studies. In simulation, for $x_{1i}\stackrel{iid}{\sim }\text{Unif}(100,200)$ and $x_{2i}\stackrel{iid}{\sim }\text{Unif}(30,500)$ we generate a variable of interest y_i according to the following error distribution of the super-population model,

• Normal distribution: y_i = β₀ + β₁x_1i + β₂x_1i + ε_i, ${\varepsilon }_i\stackrel{iid}{\sim }N(0,900)$.

  Here we use β₀ = 10, β₁ = 5 for univariate auxiliary variable case and β₀ = 10, β₁ = 5, β₂ = 3 for bivariate auxiliary variable case.

• Gamma distribution: $y_i\stackrel{iid}{\sim }\mathrm{\Gamma }(10,{\mu }_i/10)$.

  Here, we use μ_i = exp(0.01 + 0.03x_1i) for univariate auxiliary variable case and μ_i = exp(0.01 + 0.03x_1i + 0.002x_2i) for bivariate auxiliary variable case.

• Log-normal distribution: $y_i\stackrel{iid}{\sim }\text{LN}({\mu }_i,0.1)$ Here we use μ_i = 0.01 + 0.03x_1i for univariate auxiliary variable case and μ_i = 0.01 + 0.03x_1i + 0.002x_2i for bivariate auxiliary variable case.

From N = 10, 000 population data, n = 500 samples are drawn using simple random sampling. In order to remove the influence of outliers, 10,100 data are generated and then 100 largest valued data are deleted in the samples of the gamma distribution and log-normal distribution.

We consider three response probability models,

• exponential response: p_i = exp (a₀ + a₁y_i)

• linear response: p_i = b₀ + b₁y_i

• logistic response: $\text{ln\hspace{0.17em}}\left(\frac{p_i}{1-p_i}\right)=c_0+c_1{y}_i$

For the linear response probability model, p_i = b₀ + b₁y_i, p_i ∈ [0, 1], let ${\pi }_y^{\text{min}}$ and ${\pi }_y^{\text{max}}$ denote the response probabilities corresponding to the minimum and maximum values of y_i, respectively. Using given response probabilities, for instance $({\pi }_y^{\text{min}},{\pi }_y^{\text{max}})=(0.9,0.3)$, we calculate b₀, b₁ and then calculate p_i. We consider four cases of response probabilities,

• Case I $({\pi }_y^{\text{min}},{\pi }_y^{\text{max}})=(0.9,0.3)$,

• Case II : $({\pi }_y^{\text{min}},{\pi }_y^{\text{max}})=(0.8,0.3)$,

• Case III :$({\pi }_y^{\text{min}},{\pi }_y^{\text{max}})=(0.3,0.8)$,

• Case IV :$({\pi }_y^{\text{min}},{\pi }_y^{\text{max}})=(0.3,0.9)$.

Response data are obtained according to the calculated response probability p_i. Similarly, the exponential response probability model and the logistic response probability model are used to generate response data.

Three PSA estimators, ${\widehat{\overline{Y}}}_P,{\widehat{\overline{Y}}}_D$ and ${\widehat{\overline{Y}}}_C$ are calculated according to the sample weight estimation method defined in (2.2), (2.3) and (2.4). In (2.3) and (2.4) for univariate case, the number of strata L = 20 is used. For bivariate case, we use L₁ = 5 for x_1i and L₂ = 4 for x_2i. Also three bias corrected PSA estimators, ${\widehat{\overline{Y}}}_P^{\text{BC}},{\widehat{\overline{Y}}}_D^{\text{BC}}$ and ${\widehat{\overline{Y}}}_C^{\text{BC}}$ defined in section 3.4 are calculated. The performance of the estimators is compared using the following comparison statistics; bias, absolute relative bias (ARB), and root mean squared error (RMSE) defined by

Here, R = 1, 000 is used and the comparison statistics are calculated by generating a new population for each iteration. This is to reduce the influence of the specific population generated, so the true value of the r^th repeat population mean is denoted as Ȳ_r.

Tables 1 to 3 contain the results of the PSA estimators. Table 1 shows the results when the super-population model is linear and the error distribution is normal with three response probability models. Even though there are some differences of response rates depending on the response probability cases, in Case I and IV, the response rate of about 53 – 65% is obtained, and in Case II and III, the response rate of about 50 – 56% is obtained. The bias of ${\widehat{\overline{Y}}}_P$ in Case I and II shows smaller than that in the other cases and shows the smallest in the exponential response model and the largest in the logistic response model. Also we have similar results of ${\widehat{\overline{Y}}}_P$ in terms of ARB and RMSE regardless of the response model.

Comparing ${\widehat{\overline{Y}}}_P,{\widehat{\overline{Y}}}_D$ and ${\widehat{\overline{Y}}}_C$ in terms of bias, in some cases ${\widehat{\overline{Y}}}_P$ gives better results than the others, however ${\widehat{\overline{Y}}}_D$ and ${\widehat{\overline{Y}}}_C$ have good results in terms of ARB and RMSE with a very large difference. Also comparing ${\widehat{\overline{Y}}}_D$ and ${\widehat{\overline{Y}}}_C,{\widehat{\overline{Y}}}_C$ is superior to ${\widehat{\overline{Y}}}_D$, although it is not a big difference in the results of ARB and RMSE. In particular, ${\widehat{\overline{Y}}}_C$ is superior to ${\widehat{\overline{Y}}}_D$ in the bias results. Therefore, based on the results of Table 1, we can conclude that the proposed response probability estimation method gives the best results.

Table 2 shows the results when the super-population model is a log-linear model and the error distribution is a gamma distribution. Also three response probability models are considered. In Case I and II, the response rate of about 65 – 80% is obtained, and in Case III and IV, the response rate of about 40 – 45% is obtained. This result comes from the asymmetry of the gamma distribution. Investigating the results of Table 2, one can see that large biases occur in all estimators and the bias has a very large effect on the RMSE. As shown in Table 1, ${\widehat{\overline{Y}}}_C$ is the best through all comparison statistics results in Table 2.

Table 3 shows the results where the super-population model is a log-linear model and the error distribution is log-normal with three response probability models. The response rate is very similar to the result of gamma distribution. Also in Table 3, ${\widehat{\overline{Y}}}_C$ is the best through all comparison statistics results.

Table 4 shows the results of the bias corrected PSA estimator using the linear response probability model. To obtain the bias corrected PSA estimator, known response probability model and known super-population model are required. In this simulation, we use that the error distribution of the linear super-population model is normal and the response probability model is linear. Also we consider gamma and log-normal distributions as the error distribution of the log-linear super-population model.

Through the results, ${\widehat{\overline{Y}}}_P^{\text{BC}}$ gives the worst results in terms of ARB and RMSE and even the results are worse than those of ${\widehat{\overline{Y}}}_C$ in the linear response probability model of Table 2. Comparing ${\widehat{\overline{Y}}}_D^{\text{BC}}$ and ${\widehat{\overline{Y}}}_C^{\text{BC}}$, we can see that ${\widehat{\overline{Y}}}_D^{\text{BC}}$ and ${\widehat{\overline{Y}}}_C^{\text{BC}}$ have similar results in terms of ARB and RMSE. However, the bias of ${\widehat{\overline{Y}}}_C^{\text{BC}}$ is smaller than that of ${\widehat{\overline{Y}}}_D^{\text{BC}}$. Also comparing ${\widehat{\overline{Y}}}_C$ and ${\widehat{\overline{Y}}}_C^{\text{BC}}$, we have that ${\widehat{\overline{Y}}}_C^{\text{BC}}$ have better results with a very large difference. Therefore, it is very reasonable to use a bias corrected PSA estimator when the error distribution is known and the response probability model is linear. Of course, the bias corrected PSA estimator can be applied if bias is known.

Table 5 and Table 6 show the results of bivariate auxiliary variable case. The response rates are similar to the univariate variable cases. In Table 5, the result of the PSA estimator, it can be seen that the bias of ${\widehat{\overline{Y}}}_P$ is relatively small in the normal super-population result. However, the biases of ${\widehat{\overline{Y}}}_P$ in the gamma distribution and log-normal distribution are larger than those of other estimators. Also, in terms of ARB and RMSE, ${\widehat{\overline{Y}}}_P$ gives inferior results when compared to other estimators, so it is not good to use ${\widehat{\overline{Y}}}_P$ for a bivariate auxiliary variable case. Comparing ${\widehat{\overline{Y}}}_D$ and ${\widehat{\overline{Y}}}_C$ in the normal super-population result, it can be seen that the bias of ${\widehat{\overline{Y}}}_C$ is smaller. In addition, ${\widehat{\overline{Y}}}_C$ shows good results in terms of ARB and RMSE. This pattern of results also can be seen in the gamma distribution and log-normal distribution. Therefore, we conclude that ${\widehat{\overline{Y}}}_C$ gives the best result when using PSA estimator in the bivariate auxiliary variable case.

Table 6 shows the results of the bias corrected PSA estimator. As in the case of univariate auxiliary variable, ${\widehat{\overline{Y}}}_P^{\text{BC}}$ gives the best result in terms of all comparison statistics in the case of bivariate auxiliary variable. Therefore, it is appropriate to use ${\widehat{\overline{Y}}}_D^{\text{BC}}$ when using the bias corrected PSA estimator.

By comparing the results of Table 5 and Table 6, the effect of bias correction can be confirmed. The effect of bias correction of ${\widehat{\overline{Y}}}_P^{\text{BC}}$ can be seen in the gamma distribution and the log-normal distribution except for the normal distribution. On the other hand, in the results of ${\widehat{\overline{Y}}}_D^{\text{BC}}$ and ${\widehat{\overline{Y}}}_C^{\text{BC}}$, the bias correction effect can be seen in all distributions. Also the bias correction effect of ${\widehat{\overline{Y}}}_C^{\text{BC}}$ is seen for all distributions, especially the gamma and log-normal distributions. Therefore, it can be confirmed through simulation studies that the bias corrected PSA estimator gives good results when a known super-population model is used in the linear response probability model.

Recently, a lot of non-ignorable non-response has occurred in sample survey and several studies have been conducted to properly deal with it. Non-ignorable non-response is not easy to deal with because it causes bias. In particular, in order to properly handle non-ignorable non-responses, accurate estimation of response probability plays an important role. However, since the response probability is a function of the variable of interest, it is not easy to accurately estimate the response probability in practice. Therefore, practically the method of estimating the response probability should use available auxiliary variables. In this study we propose a new response probability estimation method using available auxiliary variables. It is confirmed that the proposed method provides better result than the existing methods and the bias corrected PSA estimator produces significantly better results.

Here we note that the results in this study are obtained using a linear response probability model with known super population models. A model misspecification of the response probability model is a big concern and the linear model may be neither well motivated nor appropriate in practice. Also the super population model plays an important role and a model missspecification of the super population model is also concerned. Therefore, it is necessary to study a method that can be used for an arbitrary unknown response probability model and an arbitrary unknown super-population model.