Bayesian inference for a variance component model using pairwise composite likelihood with survey data
Section 4. Extension to unequal probability sampling designs

• Previous
• Next

An important extension of our setting is to a complex sampling framework, where frequentist parameter estimation through estimation of a population-level pairwise composite likelihood is now in fairly common use. RVH and YRL have shown that an approach based on applying a frequentist pairwise composite likelihood works well for estimating multilevel model variance components in the case of certain unequal probability sampling designs, and avoids the issue of inconsistency when the second stage sample sizes are small. The uncertainty estimation in this approach uses estimating function theory and may not require the adjustments we consider in this paper. However, it would be desirable to formulate a Bayesian counterpart of this method. If a Bayesian formulation were agreed upon, the results of our paper would predict a need for adjustment of the pseudo-log-pairwise-composite-likelihood to align it with an appropriate log full likelihood function.

Suppose that the purpose is still analytic, that the model for $Y_{ij}$ is (1.1), and the objects of inference are the mean $\theta$ and the variance component ${\sigma }_u^2$ or its square root. The survey population has $N$ first stage units with sizes $M_i,$ $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,N,$ and the first-stage sample consists of $n$ of these, selected with an unequal probability sampling design. At the second stage, $m_i$ elementary units are selected by simple random sampling from the $i^{\text{th}}$ first stage unit, if that unit has been sampled at the first-stage. If the sizes $M_i$ and $m_i$ and the sampling design probabilities $p\left(s\right)$ (where $s$ runs through the two-stage subsets of the population satisfying the sample size specifications) do not depend on the $u_i$ or $e_{ij}$ values, the likelihood function can be taken to be of the form of (2.3), with $m$ replaced by $m_i,$ and the extension of our work is straightforward in principle. However, if the sizes or sampling design probabilities do depend on the values of $u_i$ or $e_{ij},$ they will be informative about the parameters of interest. The sample-level likelihood function from the combination of multilevel model and sampling design may be ill-defined or intractable. From a Bayesian perspective we then need to consider what can reasonably substitute for the true likelihood, and how closely that substitute can be approximated by an adjusted pairwise composite likelihood. The answers may depend upon the preferred method of using the sampling design probabilities in inference, and there are several possibilities. Pursuing these possibilities would be a fruitful avenue for future research.

One method, with limited applicability, would be based on the approach of Léon-Novelo and Savitsky (2019). Assuming single stage Bernoulli sampling (so that the sampling probabilities are fully determined by the inclusion probabilities) they model the joint distribution of the outcome variable, $Y,$ and the inclusion probability, $\pi ,$ using the model generating $Y$ from $x$ in the population and a model generating $\pi$ from $x$ and $Y.$ To make computations feasible there are restrictions on the form of this model; see their Theorem 1 and, especially, the special case in their Section 2.1.

We can extend the model in Section 2.1 of Léon-Novelo and Savitsky (2019) to two-stage cluster sampling. A further extension, i.e., replacing the sampling density of $Y$ with a pairwise composite likelihood analogous to the likelihood part of (2.6), can be made. Thus, subject to the limitations in Theorem 1 of Léon-Novelo and Savitsky (2019), there are counterparts to the posterior densities, (2.5) and (2.6), that include the inclusion probabilities.

Another method, not fully Bayesian, but perhaps the most widely applicable extension of our approach, is to consider the population (census) log likelihood function ((2.5) and (2.6) of RVH) to be correct, and formulate a corresponding census log pairwise composite likelihood function as in our Section 2. We would then try to estimate the latter from the sample using sampling weights ((4.2) of RVH), and make adjustments such as appropriate weight normalization, or “scaling” as in Pfeffermann, Skinner, Holmes, Goldstein and Rasbash (1998), and curvature adjustments to the resulting estimated log pairwise composite likelihood function. This would produce a log pseudo-pairwise-likelihood function that could be used as an approximate log likelihood function in Bayesian inference. It would yield a Bayesian counterpart to the frequentist method put forward by RVH and YRL, and would extend the method of this paper to the unequal probability sampling situation.

We have obtained some preliminary details for this second approach. That is, if ${\sigma }_u^2$ is known, analytic expressions for the full likelihood and pairwise composite likelihood are available for $\theta$ at the census level. For the partial likelihood we alter (2.8) by taking ${\sigma }_u$ fixed and add the weights $w_i$ and $w_{jk|i}$ as in (4.2) of RVH. With a locally uniform prior for $\theta ,$

$$p_{\text{PL}}\left(\theta |y\left(n\right)\right)\propto \text{exp}\left\{-\text{0}\text{.5}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n{\displaystyle \sum _{j<k}{w}_i{w}_{jk|i}\left(y_{ij}-\theta y_{ik}-\theta \right){\Sigma }_2^{-1}{\left(y_{ij}-\theta y_{ik}-\theta \right)}^{\text{T}}}}\right\}$$

where

$${\Sigma }_2^{-1}\mathrm{<mo>=</mo>}\left[\begin{array}{ll}{\sigma }^{11}\hfill & {\sigma }^{12}\hfill \\ {\sigma }^{21}\hfill & {\sigma }^{22}\hfill \end{array}\right]$$

with

$${\sigma }^{11}\mathrm{<mo>=</mo>}{\sigma }^{22}\mathrm{<mo>=</mo>}{\sigma }_e^{-2}\left(1-\frac{{\sigma }_u^2}{{\sigma }_e^2+2{\sigma }_u^2}\right)$$

and

$${\sigma }^{12}\mathrm{<mo>=</mo>}{\sigma }^{21}\mathrm{<mo>=</mo>}-\frac{{\sigma }_u^2}{{\sigma }_e^2\left({\sigma }_e^2+2{\sigma }_u^2\right)}\mathrm{.}$$

After some algebra,

$$p_{\text{PL}}\left(\theta |y\left(n\right)\right)\propto \text{exp}\left\{-\text{0}\text{.5}\frac{2{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n{\displaystyle {\sum }_{j<k}{w}_i{w}_{jk|i}}}}{{\sigma }_e^2+2{\sigma }_u^2}{\left[\theta -\frac{{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n{\displaystyle {\sum }_{j<k}{w}_i{w}_{jk|i}}}\mathrm{<mo\; stretchy="false">(</mo>}{y}_{ij}+y_{ik}\mathrm{<mo\; stretchy="false">)</mo>}/2}{{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n{\displaystyle {\sum }_{j<k}{w}_i{w}_{jk|i}}}}\right]}^2\right\}\mathrm{.}$$

Similarly, we alter (2.7) by taking ${\sigma }_u$ fixed and adding the weights. With a locally uniform prior for $\theta ,$

$$p_{\text{FL}}\left(\theta |y\left(n\right)\right)\propto \text{exp}\left\{-\text{0}\text{.5}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n{\displaystyle \sum _{j\mathrm{<mo>=</mo>1}}^m{\displaystyle \sum _{k\mathrm{<mo>=</mo>1}}^m{w}_i{w}_{jk|i}{\sigma }^{jk}\left(y_{ij}-\theta \right)\left(y_{ik}-\theta \right)}}}\right\}\mathrm{.}$$

After some algebra,

$$p_{\text{FL}}\left(\theta |y\left(n\right)\right)\propto \text{exp}\left\{-0.5{\displaystyle \sum _{i\mathrm{<mo>=</mo>}1}^n{w}_i\left\{{\displaystyle \sum _{j\mathrm{<mo>=</mo>}1}^m{w}_{j|i}{a}^{\left(1\right)}}+{\displaystyle \sum _{j\ne k}{w}_{jk|i}{a}^{\left(2\right)}}\right\}{\left(\theta -\widehat{\theta }\right)}^2}\right\}$$

where

$$a^{\left(1\right)}\mathrm{<mo>=</mo>}{\sigma }_e^{-2}\left(1-\frac{{\sigma }_u^2}{{\sigma }_e^2+m{\sigma }_u^2}\right)\mathrm{,}$$

$$a^{\left(2\right)}\mathrm{<mo>=</mo>}-\frac{{\sigma }_u^2}{{\sigma }_e^2\left({\sigma }_e^2+m{\sigma }_u^2\right)}$$

and

$$\widehat{\theta }\mathrm{<mo>=</mo>}\frac{{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n{w}_i\left[{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>}1}^m{a}^{\left(1\right)}{w}_{j|i}{y}_{ij}}+{\displaystyle {\sum }_{j\ne k}{a}^{\left(2\right)}{w}_{jk|i}\left(y_{ij}+y_{ik}\right)/2}\right]}}{{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>}1}^n{w}_i}\left[{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>}1}^m{a}^{\left(1\right)}{w}_{j\mathrm{<mo>|</mo>}i}+{\displaystyle {\sum }_{j\ne k}{a}^{\left(2\right)}{w}_{jk|i}}}\right]}\mathrm{.}$$

Choice of the scaling of the weights will be important. To quantify the overstated precision in the log pairwise composite posterior a numerical evaluation may be required.

An advantage of pursuing extensions of this Bayesian approach further in future research would be that it is focused on inference for the model parameters rather than on finite population quantities, and thus it would not be necessary to bring third- or fourth-order inclusion probabilities into uncertainty estimation for ${\sigma }_u^2$ or ${\sigma }_u.$

• Previous
• Article main page
• Next