COMPARISONS OF SOME WEIGHTING METHODS FOR NONRESPONSE ADJUSTMENT

Sample and population auxiliary information have been demonstrated to be useful and yield approximately equal results in large samples. Several functional forms of weights are suggested in the literature. This paper studies the properties of calibration estimators when the functional form of response probability is assumed to be known. The focus is on the difference between population and sample level auxiliary information, the latter being demonstrated to be more appropriate for estimating the coefficients in the response probability model. Results also suggest a two-step procedure, using sample information for model coefficient estimation in the first step and calibration estimation of the study variable total in the second step.


Introduction
Weighting is widely applied in surveys to adjust for nonresponse and correct other nonsampling errors.The literature contains many different proposals for nonresponse weighting methods.These methods usually treat the set of respondents as a second-phase sample [2], the elements of the response set being tied to a twofold weight compensating for both sampling and nonresponse.These weights, in particular those for nonresponse adjustment, are constructed with the aid of auxiliary information.
Treating the response set as a random subset of the sample set justifies associating each respondent with a probability of being included in the response set.Estimating this probability with aid the of auxiliary information and multiplying it by the sample inclusion probability gives an estimate of the probability of having a unit in the response set.The observations of target variable values are weighted by the reciprocals of these estimated probabilities and summed over the set of respondents, giving an estimated population total.This is known as direct nonresponse weighting adjustment [13].One example of this method is the cell weighting approach described by [11].
Alternatively, the auxiliary information is incorporated into the estimation such that the secondphase weight adjustments are determined implicitly.Such estimators are known as nonresponse weighting adjustments (see [12]), and one example is the calibration method suggested by [18].[5] combine the two approaches.They assume the response probability function to be known, and calibration serves as the means of estimating the parameters of this function.Once the parameters have been determined, the inverse of the estimated response probabilities are used as nonresponse adjustment factors.
The main feature of the calibration approach is to make the best use of available auxiliary information.When the response mechanism is assumed to be known and of the form p(•; g), parameter g is deemed a nuisance parameter [14]; this means that, although the information associated with its estimator ĝ is important, the primary objective is to estimate the target, say, the total Y = ∑ U y k .Using calibration to estimate the unknown parameters confers a different meaning on the estimation problem, in the sense that auxiliary variables are selected to provide good auxiliary information for Lithuanian Statistical Association, Statistics Lithuania Lietuvos statistiku ˛s ąjunga, Lietuvos statistikos departamentas ISSN 2029-7262 online Comparisons of some weighting methods for nonresponse adjustment estimating the parameters with good precision.This will in turn imply good precision for the estimates of response probabilities.Thus, when the response probability function is known, our principle is to view the problem of estimation in two distinct moments: estimation of parameters and estimation of targets respectively.As noted in [4], the probabilities to respond are usually functions of the sample and survey conditions, that is, the response probability for a specific individual may change when the survey conditions also well change (see also [3]).However, the mechanism leading to response/nonresponse for a sampled individual is generally not known [14].Thus, estimation in the presence of nonresponse requires some kind of modeling, explicitly or implicitly (see [5]).An implicit modeling for nonresponse adjustment can be found in [1], while [12] gives an example of explicit modeling.This paper considers nonresponse adjustment methods when the response probability function is assumed to be known up to a set of unknown coefficients.Under this assumption, direct weighting estimators can be used when the response probability model is estimated using, for example, the maximum likelihood estimator.An alternative here is to estimate the response probability model using calibration, as suggested by [5].This calibration estimator requires only the values of the covariates in the response model for the sample units in the response set, while maximum likelihood needs the values of those variables for the whole sample.One issue considered is the level of information used in calibration.An option is to use either sample or population level information when calibrating for response probability coefficient estimates.This paper contributes by demonstrating that the asymptotic variance of the coefficient estimator is smaller when sample level information is used.A simulation study is performed in order to investigate the properties of the estimators for small sample sizes.We also suggest a two-step procedure in which sample level information is used for response probability model estimation in the first step, and population level information is used for estimating population characteristics in the second step.Furthermore, the importance of correlating auxiliary variables with model and study variables is addressed.
The simulation study performed is based on data from a survey on real estate, and the bias and variance properties of the estimators are considered.Several estimators are studied, including the Horvitz-Thompson (HT) estimator using true model coefficients, direct weighting using maximum likelihood (ML) estimates of coefficients, and calibration-estimated coefficients, where calibration uses sample or population information.Two-step estimators using ML-estimated and calibration-estimated coefficients, respectively, are included, as is the linear calibration (LC) estimator [21].
The estimators studied are introduced in the next section.Section 3 compares the variance of the model parameter calibration estimators when based on population and sample level information.The results of a simulation study are reported in Section 4, and a discussion of the findings is saved for the final section.

Estimators under nonresponse
Sample s of size n is drawn from the population U = {1, 2, ..., k, ..., N} of size N using a probability sampling design, p(s), yielding first and second order inclusion probabilities π k = Pr(k ∈ s) > 0 and π kl = Pr(k, l ∈ s) > 0, respectively, for all k, l ∈ U. Let r ⊂ s denote the response set.Units in the sample respond independently with a probability p k = Pr(k ∈ r |k ∈ s) >0, for the known functional form p k = p(z t k g) evaluated at g = g ∞ , an interior point of the parameter space g ∈ G, and z k is a vector of model variables.Both g and z k are column vectors of dimension K. Furthermore, we assume that conditional on the auxiliary variables, the response probability is independent of the survey variable of interest, which is known as MAR assumption (e.g.[23]).Define the indicators: The survey variable of interest is y, and its population total, Y = ∑ U y k , is to be estimated.We can then construct an estimator for Y of the form: The weights, w k , can be defined in various ways but usually have the form w k = d k v k , where d k = 1/π k is the design weight and v k is a factor adjusting for example, for nonresponse.These factors make use of auxiliary information.The auxiliary vector is x k , with dimension P × 1, where P ≥ K and X = ∑ U x k denotes its population total.

a. Direct nonresponse weighting adjustment
One alternative of weights w k in (1) is given by w k = d k h(z t k ĝ), where h(•) = p −1 (•) and ĝ is an estimator of g ∞ .Assume p(z t k g) to be differentiable w.r.t.g and define the weighted log likelihood function of the response distribution The first order conditions for the maximum likelihood estimator (MLE) are given by ∂l(g) ∂g The first order conditions in (3) are nonlinear in g in general, and a numerical optimization method, such as the Newton-Raphson algorithm, is required to obtain the desired ĝML .Observe that ∂l(g) ∂g results in a K-dimensional column vector of partial derivatives, each with respect to one component of g.For matrix derivations, see [19].With a calculated ĝML , the estimator (1) takes the form where the subscript (DN_ML) stands for direct nonresponse weighting by ML.This estimator is asymptotically unbiased for the population total Y under the assumptions established for Theorem 1 by [13].
b.The propensity score calibration estimation [5] propose a calibration direct nonresponse adjusted estimator (1), where the weights w k are the products of the design weight and the reciprocal of the estimated response probability p(z t k ĝCAL ) for the element k in r, i.e., w k = d k h(z t k ĝCAL ), so that the estimator (1) becomes This estimator is similar to (4) in form but makes use of calibration for the estimation of g ∞ instead of ML.The strategy is to estimate g ∞ using the solution to the calibration equation Assuming h(z t k g) to be twice differentiable, [5] suggest an estimator defined by minimizing an objective function derived from (6), assuming the difference e = X − ∑ r d k h(z t k g ∞ )x k to be asymptotically normal distributed.Here, we do not impose normality assumption and derive their estimator slightly differently.
Assume that P ≥ K and define the distance function as Let Σ n be a P × P symmetric nonnegative definite matrix converging in probability to the positive definite matrix Σ, when the sample size grows arbitrarily large.Construct a weighted quadratic distance as follows: Then, the [5] estimator of g ∞ is defined as the minimizer of (8).Note that this estimator is a generalized method of moments (GMM) estimator, where minimizing (8) entails solving the estimating equations ( [7], p. 378) that results in the equation after an initial guess ĝc0 where, h(a) is the first derivative of h(a) and d(g) is assumed to be of full rank.Section 3 provides some details in the derivation of (10).
The [5] propensity calibration estimator is obtained upon the convergence of (10) and is given by: c.The linear calibration estimator The LC estimator is defined as the estimator (1) with the weights, w k , satisfying the calibration constraint where , and z k is a variable vector with the same dimension as x k .z k is assumed known at least up to the set of respondents and is called an instrument vector if it differs from x k .This system yields the vector λ The linear calibration estimator for the total Y is then given by In this setting, no explicit modeling for response or outcome is required.Instead, the method relies on the strength of the available auxiliary information.Although this is not the basic tenet, the v k factor gives the impression of a linear approximation of the reciprocal of the response probability in the sense that a good linear approximation of h(z t k g) brings about a linear calibration estimator with good statistical properties (see [15]).
d.The two-step calibration estimator [21] describe the two-step calibration approach.The first-and second-step weights are constructed according to the principle of combining population and sample levels auxiliary information.In the first step, sample level information is used to construct preliminary weights, w 1k , such that ∑ r w 1k x s k = ∑ s d k x s k , where x s k is a J-dimensional column vector of auxiliary variables with known values for all sampled units.In the second step, weights w 1k replace the design weights in the derivation of the single step calibration estimator (14), and the final weights, w k , satisfy ∑ r w k x k = X.Here, , with x U k being a P-dimensional column vector of auxiliary variables with known values for all respondents; moreover, their population totals are also known.
[16] also suggest a two-step calibration estimation assuming the known functional form of the response mechanism.The estimation process is conceptually different from the one suggested in [21], where the second-step weights are based on the first-step weights.The prediction approach supports the estimation setting suggested by [16].
Here, the concept of two-step estimation is implemented differently to ( [21], p. 88).As in [16], we assume a specified response mechanism, p(z t k g), where initial weights are calculated as w 1k = d k h(z t k ĝ) after calculating ĝ.Depending on whether the auxiliary vector z k is known up to the response set or the sample gives different options for the estimators of the true value of g.For example, if z k is known up to the sample level, then ĝ may be the MLE.If z k is known only up to the response set level, ĝ is estimated using calibration against sample level information, i.e., ∑ In the second step, the population auxiliary data are employed for estimating targets.That is, the second step weights, w k , are given by 3. Asymptotic variance of the estimated response model parameters [12] and [13] provide analytical and empirical justification for the efficiency gain when using estimated response probabilities in place of the true response probabilities, proving what had been noted by [20], namely, the estimated probabilities outperform true probabilities.[12] and [13] demonstrate this feature in a context of direct and regression adjustments where the scores are estimated using an ML procedure.This efficiency gain by using estimated probabilities can be interpreted as resulting from the lack of the location-invariance property of the HT estimator (e.g.[9], p. 10).Using true response probabilities, observations are given weights equal to the reciprocal of the probability of having the unit in the response set.However, the size of the response set is random due to nonresponse, meaning that it is not location invariant.When using ML-estimated response probabilities, estimates satisfy moment conditions at the sample level.This can be expected to reduce variance but will not in general yield an invariance property.Similar to the difference between true and estimated response probabilities, the difference between population and sample level information in the calibration estimator is considered.The precision of model parameters can be expected to affect the precision of target variable estimates.Here precision is auxiliary information dependent.As noted in [4] and [24], the strength of the relationships between the auxiliary variables and the response probabilities or study variables is crucial for the efficient performance of the weighting adjustment methods.Auxiliary information may be available at different levels, such as the population or sample levels [8].Under nonresponse, this auxiliary information is used for correcting nonresponse bias and reducing the variance of the estimator.In particular, as [23] states, sample level information is suited for nonresponse adjustment rather than variance reduction, because nonresponse affects only the location of means and not their variation.
According to the quasi-randomization setup, response set generation is an experiment made conditional on the sample.On the other hand, calibrating weights against population level information means that estimation is made unconditional on the sample.Calibration based on sample level information is therefore expected to yield more efficient estimators of response probability parameters.
Reformulating the calibration equation as illustrates that calibration against population level information brings a source of uncertainty that does not depend on the response probability distribution, i.e., variation due to the first phase sampling represented by the first term of the right-hand side of this equation.Calibrating against sample level information excludes this term, and the single source of randomness involved is the one defined by the conditional response distribution.
For more formal results, assume the asymptotic framework in which both the sample and population sizes are to increase to infinity (see, [10]), and assume further that the minimizer of (8) is consistent.
Using result 9.3.1 in [22], the covariance matrix of d(g) evaluated at the true value g = g ∞ is given by where, x k x t k , with the expectations being taken jointly with respect to the sampling design p(s) and the response distribution p(z t k g).
Consider equation (9) with g replaced by its solution ĝ, and apply the mean value theorem to decompose d(ĝ), obtaining the following equation: Then, we can substitute d(ĝ) in (9) by the r.h.s of (16) and get: where, ḡ lies in the segment between ĝ and g ∞ .We can rewrite (17) as: Under appropriate assumptions, we have that , where a stands for ĝ, ḡ or g ∞ .We replace d in (17) by its corresponding limit and obtain the asymptotic variance of the estimated model parameters as: where Π is the probability limit of which is equivalent to expression (9.80) in [7].Observe that equation (10) results from (17) after replacing d(a) with the computable entity d(ĝ c0 ) = ∑ r d k h(z t k ĝc0 )x k .Now, for calibration at the sample level, rewrite equation (7) as The conditional expectation of d s (g ∞ ) with respect to the response distribution is zero.This implies that the covariance (15) in this case is Then, with arguments similar to those that led to (20) results in asymptotic variance of the response model parameters given by The additional asymptotic variance introduced by calibrating against population level instead of sample level information is expressed by the difference The positive definiteness of the difference (23) is illustrated by the positive definiteness of the difference [6]).This is equivalent to demonstrating that because Π 1 and Π 2 are both positive definite matrices, unless h(z t k g ∞ ) = 1 for all elements in the population, and D is a full rank matrix as a consequence of (11).Observe that proving (24) is in turn equivalent to demonstrating that (Π 1 + Π 2 ) − Π 2 > 0. Thus, inequality (23) follows.

Simulation study
Under assessment are the estimators described in points "a" to "d": the direct nonresponse weighting adjustment (a), the propensity score calibration estimator (b), the linear calibration estimator (c), and the two-step calibration estimator (d).We used data from a real case study with 4228 sampled elements, of which 1783 were nonrespondents.A two-covariate logistic regression was fitted based on this data and used as the true response probability model in the simulations.Next, we created a synthetic population based on the 2445 respondents to the survey; samples were drawn from this population, after which a response set was generated using the estimated response probability model.
Five variables were selected for the study, one categorical and the others numerical.The numerical variables were transformed into logarithmic scales to reduce variability.The categorical variable, denoted γ, was a stratum indicator in the original study having six strata, thus, γ k = (γ 1k , γ 2k , γ 3k , γ 4k , γ 5k , γ 6k ), where γ ik = 1 S i (k) and S i is the i th stratum.Figure 4 presents the relationship among the original quantitative variables transformed into logarithmic form.One of them, left untransformed, was chosen to be study variable y, and estimation concerns estimating the population total Y = 17014, having the three auxiliary variables v 1 , v 2 , and v 3 .Three simulation sets were defined using three criteria.The first criterion addresses the estimator's performance in relation to the quality of auxiliary variables, the second criterion addresses the effect of the sample size, and the last focuses on the effects of model misspecification.The response probability function is defined by the logistic regression model p(R k = 1|k ∈ s) = 1/ 1 + exp(−z t k g) , where z k = (1, v 1k ) t and the parameter values are defined by their ML fit to the original 4228 observations.The samples were selected using simple random sampling without replacement followed by Poisson sampling, in which the probability used for each Bernoulli trial was the one obtained using the response model.Each simulation result was based on 1000 replications.Initial trials with higher numbers of replications produced similar results.All estimators under study used the same samples and same response sets.The expected response rate was approximately 57%.The estimators are evaluated in terms of the relative bias (RB), standard error (SE) and mean squared error (MSE).

Correctly specified response probability model
Tables 1 -3 present the results with the model vector defined as z k = (1, v 1k ) t .In Table 1, the auxiliary vector is defined as x k = (γ 1k , ..., γ 6k , γ 1k v 2k , ..., γ 6k v 2k ) t , a setup treated as the base case.As seen in Figure 4, the auxiliary variable, v 2 , correlates well with both the model variable and the study variable.A similar auxiliary vector was defined for the results in Table 2, with the exception that v 3 replaces v 2 , that is, x k = (γ 1k , ..., γ 6k , γ 1k v 3k , ..., γ 6k v 3k ) t .Here the auxiliary variable has a moderate level of correlation with the study variable, but carries much less information on the variation of the model variable in the response probability function.The correlations of v 3 with v 1 is approximately 0.16, with y approximately 0.59.
Again a similar auxiliary vector as in Table 1 was used for the results in Table 3, but here va is used in place of v 2 , i.e. x k = (γ 1k , ..., γ 6k , γ 1k va k , ..., γ 6k va k ) t .The auxiliary variable va has approximately moderate correlation with the model variable (0.45) but low correlation with the study variable (0.04).The purpose of the simulation setup in tables 1 -3 is partly to enable the study of the differences in the effect of having a good auxiliary variable for the model variable and the study variable respectively.
Table 1: Simulated estimates of RB, SE and MSE for the total of the survey variable in case of a true response model with In tables 1 -3, one can observe that the use of true probabilities ( ŶDNTrue ) leads to estimated targets with larger variability than that of the estimated targets obtained using ML-estimated probabilities ( ŶDN_ML ).Observe that the standard error when using true probabilities is 872, which is four times more than the standard error when using estimated probabilities.Note that the results for these two estimators are the same over all three tables, because they are not defined by the benchmark variables used.
As predicted by the results in Section 3, the variance for the calibration estimator of the model coefficients is smaller when sample level information is used rather than population level information.This is observed in all three tables.Also, as expected, the ML estimator is associated with the smallest variance estimates, except the two-step estimators.The results also indicate that the variance decreases with increased correlation between the model and auxiliary variables; the variance estimates are the highest in Table 2.However, the comparison of tables 1 and 3 indicates that the correlation is not the only determinant of variance.SE and MSE estimates than does calibration using sample level information.As previously observed in all tables, the two-step calibration estimators have the smallest SE and MSE estimates.

Discussion
Simulation results are consistent with the principle that estimated probabilities outperform true probabilities in weighting for nonresponse, as was earlier known for ML-estimated probabilities (see [13]).The results presented here suggest that the gain in using estimated probabilities also holds for alternative model parameter estimators.This somewhat surprising principle is here interpreted as due to the random response set size whereby the HT estimator is not location invariant.The results presented also suggest that the gain in using estimated probabilities holds for alternative model parameters.In fact, even under the considered misspecifications of the response probability model, the results indicate the improved performance of the weighting estimators using estimated response probabilities.
The major concern in the paper is the use of sample or population level auxiliary information in the calibration of the response probability function.The simulation results obtained are consistent with the formal asymptotic argument presented, suggesting the use of sample auxiliary information for estimating the response probability function.Results indicate that the response function parameters are estimated with lower variance when using sample auxiliary information instead of population level information.The importance of having auxiliary information highly correlated with the model variables is observed for both levels of auxiliary information.
Using sample or population level information in the calibration estimators of population totals produces similar relative biases and standard errors.However, the sample-based calibration estimator has a smaller MSE than does the population counterpart; this is observed in all cases when the model is correctly specified.However, ML-estimated probabilities yield an estimator with the smallest SE and MSE estimates.
The auxiliary vector used in Table 2 is moderately correlated with the study variable while having virtually no relationship with the model variable; the standard errors for the single step calibration estimators are greater than when the auxiliary variable is correlated with both the study and model variables (Table 1).A much smaller difference is observed when auxiliary variables are correlated with the model variable while having virtually no correlation with the study variable (Table 3).This suggests a preference for auxiliary variables related to response propensity model variables over auxiliary variables related to the study variable.
Response probability function modelling is susceptible to misspecification.Under the erroneous choice of model variables, the major effects observed here are on the bias of the propensity basedestimators.The estimators (i.e., ŶDN_ML , ŶPS_pop and ŶPS_samp ) are associated with larger biases, although still at a low relative level (tables 7 and 8).The major observation is that the population-based calibrated estimator is more effective in error protection than is either the sample-calibrated or MLbased estimator.Still, good auxiliary information is important for the model variables.Although the evidence presented suggests that using sample auxiliary information is superior to using population auxiliary information in propensity calibration estimators, the population level propensity calibration is suggested to be the best alternative for reducing the MSE of the target estimates when the model is misspecified.
An erroneous functional form of the response probability model does not have a great impact on the estimator performance, according to the results in Table 9.One likely reason for this is that the two models are similar.However, the results suggest that the choice of the functional form is less important than is having the right model variables.This is partly supported by the competitive performance of the linear calibration estimator at larger sample sizes.
We suggest that estimation be performed in two steps; in the first step, the sample auxiliary data are used in the propensity calibration for estimating the response probabilities; in the second step, the products of the design weights and the reciprocals of response probabilities replace the design weights in the linear calibration estimator.The two-step estimation is to be performed using sample auxiliary information for estimating the response model through calibration, followed by the use of population auxiliary information for estimating target entities.This will generally produce more Comparisons of some weighting methods for nonresponse adjustment efficient estimates.
The results presented all favor the suggested two-step calibration estimators.In some cases, these estimators are associated with larger bias, thought their relative sizes are small.In terms of MSE, the two-step estimators outperform other estimators.This is also observed when the response probability model is misspecified.A general suggestion would be to use ML to estimate model parameters in the first step, if model variables are available at the sample level.If the model variables are available only at the response set level, the [5] calibration estimator for model parameters is almost equally good.The results of the two-step estimators are of particular interest since response probability functions used in practice are models susceptible to misspecification.The effects of misspecification are usually unknown and can yield an adjusted estimator with a larger bias than the unadjusted one, depending on correlation structures among the study variable, response probability and auxiliary variables ( [17]).Although small, a second calibration step reduces bias estimates in cases with a wrong auxiliary variable in the response probability model (tables 7 and 8).A question for further studies is whether a second step calibration can protect against the misspecification of the response probability function and/or if indicators of misspecification can be developed.
With large sample sizes and carefully chosen auxiliary information, the linear calibration estimator is fairly competitive with the propensity-based estimators.The linear calibration estimator is known to have good properties when good auxiliary information is available.On the other hand, poorly defined auxiliary variables may lead to negative and/or very large weights in the linear calibration ( [21], remark 6.1).These problems may result in very inefficient estimates.A conclusion based on the results presented in Table 1 is that the properties of the linear calibration estimator can be improved by using efficient initial weights.These weights can be derived from a sample-based propensity calibration estimator.The combined approaches produce more efficient estimates.
Tables 1 -3 provide results for, Ŷ2stepB , an estimator not discussed here.It is a version of Ŷ2stepA in which auxiliary information exists only at the sample level, i.e. sample level information is used in both steps.This will generally provide slightly better RB, but the SE and MSE are higher than those provided by Ŷ2stepA .

Limitations
In this article, we use an estimation setting in which only positive correlations among the variables in the study are considered.[17] have noted that the direction of the correlation between the variables involved in the study has an influence on the properties of the estimated entities.This suggests a further investigation whether the results here are the same when the variables are negatively correlated.

Figure 1 :
Figure 1: Pairwise correlations among the original variables used in the simulation study.Correlations calculated on the set of synthetic population.

Figure 2 :
Figure 2: Pairwise correlations among artificial and original variables.Correlations calculated on the set of synthetic population.

Table 2 :
Simulated estimates of RB, SE and MSE for the total of the survey variable in case of a true response model with z k = (1, v 1k ) t , x k = (γ k , v 3k γ k ) t and n=300