Bayesian predictive inference of a proportion under a two-fold small area model with heterogeneous correlations
Section 2. Bayesian two-fold small area models and computations

• Previous
• Next

We consider a finite population of $l$ areas and $M_i$ clusters within the $i^{\text{th}}$ area, and we assume there are $N_{ij}$ individuals in $j^{\text{th}}$ cluster within $i^{\text{th}}$ area. The binary responses are $y_{ijk}$ for $i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}j\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{M}_i\mathrm{,}$ $k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{N}_{ij}.$ We assume that a simple random sample of $m_i$ clusters is taken from the $i^{\text{th}}$ small area and a simple random sample of $n_{ij}$ individuals is taken from the $m_i$ sampled clusters from the $i^{\text{th}}$ area. Here, we assume the survey weights are the same within all clusters in each area. Let $n_i={\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{m_i}{n}_{ij}}\mathrm{,}{s}_{ij}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^{n_{ij}}{y}_{ijk}}$ and $s_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{m_i}{s}_{ij}}.$

Our target is the finite population proportion of the $i^{\text{th}}$ area which is given by

$$P_i\mathrm{<mo>=</mo>}\frac{{\displaystyle \sum _{j\mathrm{<mo>=</mo>\; 1}}^{M_i}}{\displaystyle \sum _{k\mathrm{<mo>=</mo>\; 1}}^{N_{ij}}{y}_{ijk}}}{N_i}\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}$$

where $N_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{M_i}{N}_{ij}}.$ Let $T_{ij}^{\left(1\right)}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>}{n}_{ij}+1}^{N_{ij}}{y}_{ijk}}$ denote the nonsampled totals of the sampled clusters $\left(j\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{m}_i\right),$ and $T_{ij}^{\left(2\right)}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^{N_{ij}}{y}_{ijk}},$ the totals of the nonsampled clusters $\left(j\mathrm{<mo>=</mo>}{m}_i+\mathrm{1,}\dots \mathrm{,}{M}_i\right).$ Letting $n_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{m_i}{n}_{ij}},$ ${\widehat{p}}_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^{n_{ij}}{y}_{ijk}/n_i},$ we can express our target, $P_i$ as

$$P_i\mathrm{<mo>=</mo>}\frac{n_i{\widehat{p}}_i+{\displaystyle \sum _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}{T}_{ij}^{\left(1\right)}+{\displaystyle \sum _{j\mathrm{<mo>=</mo>}{m}_i+1}^{M_i}}{T}_{ij}^{\left(2\right)}}{N_i}\mathrm{,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{.}(2.1)$$

To make inference about the $P_i,$ we fit hierarchical Bayesian models to the data. Using the beta-binomial representation, these models accommodate the two-fold design structure. We describe two models, one with homogeneous correlation and the other with heterogeneous correlations, our main contribution beyond Nandram (2015). In Section 2.1 we review the hierarchical Bayesian model with homogeneous correlation, Nandram (2015) and we show how to make it comparable to our hierarchical Bayesian model with heterogeneous correlations which we describe in Section 2.2. In Section 2.3 we describe the blocked Gibbs sampler to fit our model with heterogeneous correlations.

2.1 A review of two-fold model with homogeneous correlation

Nandram (2015) described the two-fold small area model with homogeneous correlation. Here we give a brief review of its main assumptions which are

$$\begin{array}{ll}{y}_{ijk}|p_{ij}\hfill & \stackrel{\text{ind}}{\sim }\text{Bernoulli}\left(p_{ij}\right)\mathrm{,}(2.2)\hfill \\ {\mu }_i|\theta \mathrm{,}\gamma \hfill & \stackrel{\text{iid}}{\sim }\text{Beta}\left[\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right]\mathrm{,}(2.3)\hfill \\ \rho \mathrm{,}\theta \mathrm{,}\gamma \hfill & \stackrel{\text{iid}}{\sim }\text{Uniform}\left(\mathrm{0,1}\right)\mathrm{,}(2.4)\hfill \end{array}$$

where $\rho$ and $\gamma$ represent the intracluster and intercluster correlation, respectively. It is assumed that $\mathrm{0\; <mo><</mo>}\theta \mathrm{,}\rho \mathrm{,}\gamma \mathrm{<mo><</mo>\; 1}$ strictly. Note that within the same area the intracluster correlation $\rho ,$ the correlation between two units in the same cluster, is $\text{cor}\left(y_{ijk},y_{ij{k}^{\prime }}|{\mu }_i,\gamma \mathrm{,}\rho \right)\mathrm{<mo>=</mo>}\rho \mathrm{,\; <mi>\; </mi>}k\ne k^{\prime }\mathrm{.}$ Similarly, within the same area the intercluster correlation $\gamma ,$ the correlation between two units in two different clusters, is $\text{cor}\left(y_{ijk}\mathrm{,}{y}_{i{j}^{\prime }{k}^{\prime }}|\theta \mathrm{,}\gamma \mathrm{,}\rho \right)\mathrm{<mo>=</mo>}\gamma \mathrm{,}j\ne j^{\prime }\mathrm{,}k\ne k^{\prime }\mathrm{.}$ Here, it is $\rho$ that makes a difference between the one-fold and two-fold models, and when $\rho$ goes to zero, the two-fold model becomes the one-fold model, Nandram (2015).

To fit the model specified by (2.2) $–$ (2.4), Nandram (2015) used random sampling and Gaussian quadrature to perform one-dimensional numerical integrations. He also used Gibbs sampling for comparison and found minor differences. However, our generalization to heterogeneous correlations (increased number of parameters) leads to additional weakly identified parameters and model fitting becomes more difficult. So we incorporate unimodality constraints on the prior distributions of the area parameters, thereby making it possible to analyze sparse data. To make fair comparisons between the two models, one with homogeneous correlations and the other with heterogeneous correlations, we also impose the unimodality constraints in the model specified by (2.2) $–$ (2.4). Our results under this slightly modified homogeneous model are similar to those in Nandram (2015).

The methods introduced in this paper allow unimodality to be imposed on some distributions to assist in the estimation of weakly identified parameters. The unimodality restrictions are flexible enough to avoid over-restricting the models. For a full nonparametric Bayesian procedure, see Damien, Laud and Smith (1997). Thus, throughout all our computations, we apply the unimodality restriction to hyperparameters of ${\mu }_i\left(i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\right),$

$$\frac{\gamma }{1-\gamma }\mathrm{<mo><</mo>}\theta \mathrm{<mo><</mo>}\frac{1-2\gamma }{1-\gamma }\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mtext>\hspace{0.17em}</mtext>\; 0\; <mo><</mo>}\gamma \mathrm{<mo><</mo>}\frac{1}{3}\mathrm{.}(2.5)$$

We also use similar unimodality restrictions in Section 2.2 for the model with heterogeneous correlations. Henceforth, we call the model specified by (2.2) $–$ (2.5) the HoC model.

To fit the model, Nandram (2015) use the multiplication rule by obtaining $p_{ij}$ after drawing random samples of $(\mu ,\rho \mathrm{,}\theta ,\text{and}\gamma )$ from their joint posterior density, where $\mu ={\left({\mu }_1\mathrm{,}\dots \mathrm{,}{\mu }_{\mathcal{l}}\right)}^{\prime }.$ The conditional posterior density of the $p_{ij}$ is given by

$$p_{ij}|s_{ij}\mathrm{,}{\mu }_i\mathrm{,}\rho \stackrel{\text{ind}}{\sim }\text{Beta}\left\{s_{ij}+{\mu }_i\frac{1-\rho }{\rho }\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-\rho }{\rho }\right\}\mathrm{,}$$

and letting $s_{ij}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{k\mathrm{<mo>=</mo>\; 1}}^{n_{ij}}{y}_{ijk}}$ and collapsing over the $p_{ij},$ we get

$$\begin{array}{ll}\pi \left(\mu \mathrm{,}\rho \mathrm{,}\theta \mathrm{,}\gamma |y\right)\hfill & \propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>1}}^l}{\displaystyle \prod _{j\mathrm{<mo>=</mo>1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-\rho }{\rho }\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-\rho }{\rho }\right)}{B\left({\mu }_i\frac{1-\rho }{\rho },\left(1-{\mu }_i\right)\frac{1-\rho }{\rho }\right)}\hfill \\ \hfill & \times \frac{{\mu }_i^{\theta \frac{1-\gamma }{\gamma }-1}{\left(1-{\mu }_i\right)}^{\left(1-\theta \right)\frac{1-\gamma }{\gamma }-1}}{B\left(\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right)}\mathrm{,}\mathrm{0<mo><</mo>}{\mu }_i\mathrm{,}\rho \mathrm{<mo><</mo>1,}i\mathrm{<mo>=</mo>1,}\dots \mathrm{,}l\mathrm{,}\frac{\gamma }{1-\gamma }\mathrm{<mo><</mo>}\theta \mathrm{<mo><</mo>}\frac{1-2\gamma }{1-\gamma }\mathrm{,<mtext>\hspace{0.17em}</mtext><mi>\; </mi><mi>\; </mi>0<mo><</mo>}\gamma \mathrm{<mo><</mo>}\frac{1}{3}\mathrm{.}\hfill \end{array}$$

Because $T_{ij}^{\left(1\right)}|p_{ij}\stackrel{\text{ind}}{\sim }\text{Binomial}\left(N_{ij}-n_{ij}\mathrm{,}{p}_{ij}\right)$ and $T_{ij}^{\left(2\right)}|p_{ij}\stackrel{\text{ind}}{\sim }\text{Binomial}\left(N_{ij}\mathrm{,}{p}_{ij}\right)$ and, given $p_{ij},$ $T_{ij}^{\left(1\right)}$ and $T_{ij}^{\left(2\right)}$ are independent, once samples of the $p_{ij}$ are obtained, it is easy to make Bayesian predictive inference. See Nandram (2015) for details.

2.2 A two-fold model with heterogeneous correlations

We extend the HoC model to accommodate the heterogeneous correlations. Our assumptions are

$$\begin{array}{ll}{y}_{ijk}|p_{ij}\hfill & \stackrel{\text{ind}}{\sim }\text{Bernoulli}\left(p_{ij}\right)\mathrm{,}(2.6)\hfill \\ p_{ij}|{\mu }_i\mathrm{,}{\rho }_i\hfill & \stackrel{\text{ind}}{\sim }\text{Beta}\left[{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right]\mathrm{,}(2.7)\hfill \\ {\mu }_i|\theta \mathrm{,}\gamma \hfill & \stackrel{\text{iid}}{\sim }\text{Beta}\left[\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right]\mathrm{,}(2.8)\hfill \\ {\rho }_i|\varphi \mathrm{,}\delta \hfill & \stackrel{\text{iid}}{\sim }\text{Beta}\left[\varphi \frac{1-\delta }{\delta }\mathrm{,}\left(1-\varphi \right)\frac{1-\delta }{\delta }\right]\mathrm{,}(2.9)\hfill \\ \theta \mathrm{,}\gamma \mathrm{,}\varphi \mathrm{,}\delta \hfill & \stackrel{\text{iid}}{\sim }\text{Uniform}\left(\mathrm{0,1}\right)\mathrm{.}(2.10)\hfill \end{array}$$

Note that the intracluster correlation coefficient $\rho$ introduced in the HoC model is replaced by ${\rho }_i\mathrm{<mi>\; </mi>}\left(i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\right)$ to provide the hierarchical Bayesian model with heterogeneous correlations.

Similar to the HoC model, a priori we also impose two sets of unimodality constraints,

$$\frac{\gamma }{1-\gamma }\mathrm{<mo><</mo>}\mathrm{<mtext>\hspace{0.17em}</mtext>\theta <mtext>\hspace{0.17em}</mtext>}\mathrm{<mo><</mo>}\frac{1-2\gamma }{1-\gamma }\mathrm{,<mi>\; </mi><mi>\; </mi><mtext>\hspace{0.17em}</mtext>0<mtext>\hspace{0.17em}</mtext><mo><</mo>}\gamma \mathrm{<mo><</mo>}\frac{1}{3}\text{and}\mathrm{<mi>\; </mi>}\frac{\delta }{1-\delta }\mathrm{<mo><</mo>}\varphi \mathrm{<mo><</mo>}\frac{1-2\delta }{1-\delta }\mathrm{,<mi>\hspace{0.17em}</mi><mi>\; </mi><mi>\; </mi>0<mo><</mo>}\delta \mathrm{<mo><</mo>}\frac{1}{3}\mathrm{.}(2.11)$$

Appendix B contains simple proofs of the above inequalities as unimodality criterion and how to incorporate these constraints into our computation. Henceforth, we call the hierarchical Bayesian model specified by (2.6) $–$ (2.11) the HeC model.

Again, similar to Nandram (2015), under the HeC model, we show in Appendix A that

$$\begin{array}{ll}\text{cor}\left(y_{ijk}\mathrm{,}{y}_{ij{k}^{\prime }}|{\mu }_i\mathrm{,}\gamma \mathrm{,}{\rho }_i\right)\hfill & \mathrm{<mo>=</mo>}{\rho }_i\mathrm{,}k\ne k^{\prime }\mathrm{,}(2.12)\hfill \\ \text{cor}\left(y_{ijk}\mathrm{,}{y}_{i{j}^{\prime }{k}^{\prime }}|\theta \mathrm{,}\gamma \mathrm{,}{\rho }_i\right)\hfill & \mathrm{<mo>=</mo>}\gamma \mathrm{,}j\ne j^{\prime }\mathrm{,}k\ne k^{\prime }\mathrm{.}(2.13)\hfill \end{array}$$

That is, within the $i^{\text{th}}$ area, the intracluster correlation coefficient is ${\rho }_i$ and the intercluster correlation coefficient is $\gamma .$

Using Bayes’ theorem in the HeC model, the joint posterior density $\pi \left(p,\mu ,\rho ,\theta ,\gamma ,\varphi ,\delta |y\right)$ is easy to write down. (This is the density without the normalization constant.) Henceforth, we would call this joint posterior density the HeC posterior.

In order to make inference about the finite population proportion, $P_i,$ we draw samples from $\pi \left(p,\mu ,\rho ,\theta ,\gamma ,\varphi ,\delta |y\right)$ using the multiplication rule and blocked Gibbs sampler. This procedure is described in Section 2.3.

2.3 Computations in the HeC posterior

First, note that we collapse HeC posterior over the $p_{ij}$ and then use the Gibbs sampler to fit the joint marginal posterior density. After obtaining the samples, we can draw samples of the $p_{ij}$ from the conditional posterior densities of the $p_{ij}$ by applying the multiplication rule.

As in the HoC model, the conditional posterior density of $p_{ij}$ is

$$p_{ij}|{\mu }_i\mathrm{,}{\rho }_i\mathrm{,}\theta \mathrm{,}\gamma \mathrm{,}\varphi \mathrm{,}\delta \mathrm{,}y\stackrel{\text{ind}}{\sim }\text{Beta}\left\{s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right\}\mathrm{,}\mathrm{0\; <mo><</mo>}{p}_{ij}\mathrm{<mo><</mo>\; 1.}$$

Thus, it is easy to draw samples of the $p_{ij}$ once samples are obtained from the joint posterior density of $\left(\mu ,\rho ,\theta ,\gamma ,\varphi ,\delta \right).$ After integrating out the $p_{ij}$ from the HeC posterior, the marginal joint posterior density is given by

$$\begin{array}{ll}\pi \left(\mu \mathrm{,}\rho \mathrm{,}\theta \mathrm{,}\gamma \mathrm{,}\varphi \mathrm{,}\delta |y\right)\hfill & \propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}\hfill \\ \hfill & \times \frac{{\mu }_i^{\theta \frac{1-\gamma }{\gamma }-1}{\left(1-{\mu }_i\right)}^{\left(1-\theta \right)\frac{1-\gamma }{\gamma }-1}}{B\left(\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right)}\times \frac{{\rho }_i^{\varphi \frac{1-\delta }{\delta }-1}{\left(1-{\rho }_i\right)}^{\left(1-\varphi \right)\frac{1-\delta }{\delta }-1}}{B\left(\varphi \frac{1-\delta }{\delta }\mathrm{,}\left(1-\varphi \right)\frac{1-\delta }{\delta }\right)},\mathrm{0\; <mo><</mo>}{\mu }_i\mathrm{,}{\rho }_i\mathrm{<mo><</mo>\; 1,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}\hfill \\ \hfill & \frac{\gamma }{1-\gamma }\mathrm{<mo><</mo>}\theta \mathrm{<mo><</mo>}\frac{1-2\gamma }{1-\gamma }\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mtext>\hspace{0.17em}</mtext>\; 0\; <mo><</mo>}\gamma \mathrm{<mo><</mo>}\frac{1}{3}\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>}\frac{\delta }{1-\delta }\mathrm{<mo><</mo>}\varphi \mathrm{<mo><</mo>}\frac{1-2\delta }{1-\delta }\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mtext>\hspace{0.17em}</mtext>\; 0\; <mo><</mo>}\delta \mathrm{<mo><</mo>}\frac{1}{3}\mathrm{.}\hfill \end{array}$$

The conditional posterior densities are

$$\begin{array}{ll}\pi \left({\mu }_i|{\rho }_i\mathrm{,}\theta \mathrm{,}\gamma \mathrm{,}\varphi \mathrm{,}\delta \mathrm{,}y\right)\hfill & \propto {\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}\times {\mu }_i^{\theta \frac{1-\gamma }{\gamma }-1}{\left(1-{\mu }_i\right)}^{\left(1-\theta \right)\frac{1-\gamma }{\gamma }-1}\mathrm{,}\hfill \\ \pi \left({\rho }_i|{\mu }_i,\theta \mathrm{,}\gamma \mathrm{,}\varphi \mathrm{,}\delta \mathrm{,}y\right)\hfill & \propto {\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}\times {\rho }_i^{\varphi \frac{1-\delta }{\delta }-1}{\left(1-{\rho }_i\right)}^{\left(1-\delta \right)\frac{1-\delta }{\delta }-1}\mathrm{,}\hfill \end{array}$$

and letting $G_1\mathrm{<mo>=</mo>}{\left\{{\displaystyle {\prod }_{i\mathrm{<mo>=</mo>\; 1}}^l{\mu }_i}\right\}}^{1/l}$ and $G_2\mathrm{<mo>=</mo>}{\left\{{\displaystyle {\prod }_{i\mathrm{<mo>=</mo>\; 1}}^l}\left(1-{\mu }_i\right)\right\}}^{1/l},$

$$\pi \left(\theta |\mu \mathrm{,}\rho \mathrm{,}\gamma \mathrm{,}\varphi \mathrm{,}\delta \mathrm{,}y\right)\propto {\left\{\frac{G_1^{\theta \frac{1-\gamma }{\gamma }-1}{G}_2^{\left(1-\theta \right)\frac{1-\gamma }{\gamma }-1}}{B\left(\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right)}\right\}}^l\mathrm{,}$$

and

$$\pi \left(\gamma |\mu \mathrm{,}\rho \mathrm{,}\theta ,\varphi \mathrm{,}\delta \mathrm{,}y\right)\propto {\left\{\frac{G_1^{\theta \frac{1-\gamma }{\gamma }-1}{G}_2^{\left(1-\theta \right)\frac{1-\gamma }{\gamma }-1}}{B\left(\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right)}\right\}}^l.$$

Similarly, letting $H_1\mathrm{<mo>=</mo>}{\left\{{\displaystyle {\prod }_{i\mathrm{<mo>=</mo>\; 1}}^l{\rho }_i}\right\}}^{1/l}$ and $H_2\mathrm{<mo>=</mo>}{\left\{{\displaystyle {\prod }_{i\mathrm{<mo>=</mo>\; 1}}^l}\left(1-{\rho }_i\right)\right\}}^{1/l},$

$$\pi \left(\varphi |\mu \mathrm{,}\rho \mathrm{,}\theta ,\gamma \mathrm{,}\delta \mathrm{,}y\right)\propto {\left\{\frac{H_1^{\varphi \frac{1-\delta }{\delta }-1}{H}_2^{\left(1-\varphi \right)\frac{1-\delta }{\delta }-1}}{B\left(\varphi \frac{1-\delta }{\delta }\mathrm{,}\left(1-\varphi \right)\frac{1-\delta }{\delta }\right)}\right\}}^l\mathrm{,}$$

and

$$\pi \left(\delta |\mu \mathrm{,}\rho \mathrm{,}\theta ,\gamma \mathrm{,}\varphi \mathrm{,}y\right)\propto {\left\{\frac{H_1^{\varphi \frac{1-\delta }{\delta }-1}{H}_2^{\left(1-\varphi \right)\frac{1-\delta }{\delta }-1}}{B\left(\varphi \frac{1-\delta }{\delta }\mathrm{,}\left(1-\varphi \right)\frac{1-\delta }{\delta }\right)}\right\}}^l.$$

The problem with this procedure is that $\theta$ and $\gamma$ are correlated because intuitively they both depend on only $\left\{{\mu }_i\right\}$ through two numbers, $G_1$ and $G_2,$ not the data, $y.$ This gives poor mixing in the Gibbs sampler. For instance, $E\left({\mu }_i|\theta ,\gamma \right)=\theta ,$ $\text{Std}\left({\mu }_i|\theta ,\gamma \right)=\theta \sqrt{\gamma \left(1-\theta \right)/\theta }$ and ${\mu }_i\approx \theta \left\{1+z_i\sqrt{\gamma \left(1-\theta \right)/\theta }\right\},$ where $E\left(z_i\right)=0$ and $\text{Var}\left(z_i\right)=1,$ Nandram (2015). That is, $\left\{{\mu }_i\right\}$ is correlated with $\theta$ and $\gamma .$ Similar problems occur in $\left(\rho ,\varphi ,\delta \right).$ Therefore, in order to solve these weak identifiability problems, we use the blocked Gibbs sampler to draw random samples of $\left(\mu ,\rho ,\theta ,\gamma ,\varphi ,\delta \right).$

The blocked Gibbs sampler is obtained by drawing from the conditional posterior density $\left(\mu ,\theta ,\gamma |\rho ,\varphi ,\delta ,y\right)$ and $\left(\rho ,\varphi ,\delta |\mu ,\theta ,\gamma ,y\right)$ each in turn until convergence as we describe below. The two joint conditional posterior densities are

$$\begin{array}{ll}{\pi }_1\left(\mu \mathrm{,}\theta \mathrm{,}\gamma |\rho ,\varphi ,\delta ,y\right)\hfill & \propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}\hfill \\ \hfill & \times \frac{{\mu }_i^{\theta \frac{1-\gamma }{\gamma }-1}{\left(1-{\mu }_i\right)}^{\left(1-\theta \right)\frac{1-\gamma }{\gamma }-1}}{B\left(\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right)},\mathrm{0\; <mo><</mo>}{\mu }_i\mathrm{<mo><</mo>\; 1,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}\frac{\gamma }{1-\gamma }\mathrm{<mo><</mo>}\theta \mathrm{<mo><</mo>}\frac{1-2\gamma }{1-\gamma }\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mtext>\hspace{0.17em}</mtext>\; 0\; <mo><</mo>}\gamma \mathrm{<mo><</mo>}\frac{1}{3}\hfill \end{array}$$

and

$$\begin{array}{ll}{\pi }_2\left(\rho ,\varphi ,\delta |\mu \mathrm{,}\theta \mathrm{,}\gamma ,y\right)\hfill & \propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}\hfill \\ \hfill & \times \frac{{\rho }_i^{\varphi \frac{1-\delta }{\delta }-1}{\left(1-{\rho }_i\right)}^{\left(1-\varphi \right)\frac{1-\delta }{\delta }-1}}{B\left(\varphi \frac{1-\delta }{\delta }\mathrm{,}\left(1-\varphi \right)\frac{1-\delta }{\delta }\right)},\mathrm{0\; <mo><</mo>}{\rho }_i\mathrm{<mo><</mo>\; 1,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}\frac{\delta }{1-\delta }\mathrm{<mo><</mo>}\varphi \mathrm{<mo><</mo>}\frac{1-2\delta }{1-\delta }\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mtext>\hspace{0.17em}</mtext>\; 0\; <mo><</mo>}\delta \mathrm{<mo><</mo>}\frac{1}{3}.\hfill \end{array}$$

To run the blocked Gibbs sampler, we apply the multiplication rule in ${\pi }_1\left(\mu ,\theta ,\gamma |\rho ,\varphi ,\delta ,y\right)$ and ${\pi }_2\left(\rho ,\varphi ,\delta |\mu ,\theta ,\gamma ,y\right);$ see, for example, Molina et al. (2014) and Toto and Nandram (2010).

First, we consider ${\pi }_1\left(\mu ,\theta ,\gamma |\rho ,\varphi ,\delta ,y\right).$ We integrate out $\mu$ and obtain the joint conditional posterior density of $\left(\theta ,\gamma \right)$ given $\rho ,\varphi ,\delta$ and $y,$

$$\begin{array}{ll}p\left(\theta ,\gamma |\rho ,\varphi ,\delta ,y\right)\hfill & \propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}\{{\displaystyle {\int }_0^1\left[{\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}\right]}\hfill \\ \hfill & \times \frac{{\mu }_i^{\theta \frac{1-\gamma }{\gamma }-1}{\left(1-{\mu }_i\right)}^{\left(1-\theta \right)\frac{1-\gamma }{\gamma }-1}}{B\left(\theta \frac{1-\gamma }{\gamma }\mathrm{,}\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right)}d{\mu }_i\},\mathrm{0\; <mo><</mo>}{\mu }_i\mathrm{<mo><</mo>\; 1,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}\hfill \\ \hfill & \frac{\gamma }{1-\gamma }\mathrm{<mo><</mo>}\theta \mathrm{<mo><</mo>}\frac{1-2\gamma }{1-\gamma }\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mtext>\hspace{0.17em}</mtext>\; 0\; <mo><</mo>}\gamma \mathrm{<mo><</mo>}\frac{1}{3}.\hfill \end{array}$$

Here, the middle Riemann sum method is used to integrate out all ${\mu }_i,i=1,\dots ,l.$ We partition the interval (0, 1) into $G$ subintervals $\left(a_0,a_1\right],\left(a_1,a_2\right],\dots ,\left[a_{G-1},a_G\right],$ where $a_0=0,a_i=i/G,i=1,\dots G.$ Then we can compute the joint conditional posterior distribution of $\left(\theta ,\gamma \right),$ as follows.

$$p\left(\theta ,\gamma |\rho ,\varphi ,\delta ,y\right)\propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}\left[\underset{G\to \infty }{\lim }{\displaystyle \sum _{v=1}^G{g}_i\left(\frac{a_{v-1}+a_v}{2}\right)}\left\{F_1\left(a_{v-1}\right)-F_1\left(a_v\right)\right\}\right],$$

$$g_i\left({\mu }_i\right)={\displaystyle \prod _{j=1}^{m_i}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}}$$

and $F_1(\cdot )$ is the cdf corresponding to $f_1(\cdot )$ which is a density function of $\text{Beta}\left(\theta \frac{1-\gamma }{\gamma },\left(1-\theta \right)\frac{1-\gamma }{\gamma }\right).$ Next, we also integrate out $\theta$ by using Gaussian quadrature via Legendre orthogonal polynomials,

$$p\left(\gamma |\rho ,\varphi ,\delta ,y\right)\approx {\displaystyle \sum _{g=1}^G{\omega }_g}\left\{{\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle {\int }_0^1{\pi }_1\left({\mu }_i,x_g,\gamma |{\rho }_i,\varphi ,\delta ,y\right)d{\mu }_i}\right\},$$

where $\left\{{\omega }_g\right\}$ are the weights and $\left\{x_g\right\}$ are roots of the Legendre polynomial with the interval $\left[\frac{\gamma }{1-\gamma },\frac{1-2\gamma }{1-\gamma }\right].$ We have taken $G=20$ in our computations (larger values of $G$ make little difference).

Now, we can use a univariate grid method (e.g., Molina, Nandram and Rao 2014 and Toto and Nandram 2010) in order to draw samples of the posterior density of $\gamma$ conditional on $\rho ,\varphi ,\delta$ and $y;$ see Ritter and Tanner (1992) for a description of the griddy Gibbs sampler. Then, conditional on $\gamma ,$ we get the posterior density of $\theta ,$ as follows,

$$p\left(\theta |\gamma ,\rho ,\varphi ,\delta ,y\right)\approx {\displaystyle \sum _{g=1}^G{\omega }_g}\left\{{\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle {\int }_0^1{\pi }_1\left({\mu }_i,\theta |\gamma ,{\rho }_i,\varphi ,\delta ,y\right)d{\mu }_i}\right\}.$$

Samples are obtained from the conditional posterior density of $\theta$ by using the univariate grid sampler again. Subsequently, conditional on $\left(\theta \mathrm{,}\gamma \right),\mu$ is drawn from $p\left(\mu |\theta ,\gamma ,\rho ,\varphi ,\delta ,y\right)$ using the univariate grid sampler.

For the grid method, we divide the unit interval into sub-intervals of 0.01 width, and the joint posterior density is approximated by a discrete distribution with probabilities proportional to the heights of the continuous distribution at the mid-points of these sub-intervals. Note that a uniform jittering is done within each selected interval to allow different deviates with probability one (Nandram 2015). Even when we used finer sub-intervals (e.g., using 0.005 width), the inference results turned out to be almost same. Thus, we use the sub-intervals of 0.01 width; see Molina et al. (2014). When most of the distribution is near one of the boundaries (e.g., 0 or 1), we make intervals with much smaller widths to capture small or large values of the parameter.

Second, we consider ${\pi }_2\left(\rho ,\varphi ,\delta |\mu ,\theta ,\gamma ,y\right).$ We integrate out $\rho$ and obtain the joint conditional posterior density of $\left(\varphi ,\delta \right)$ given $\mu ,\theta ,\gamma$ and $y,$

$$\begin{array}{ll}p\left(\varphi ,\delta |\mu ,\theta ,\gamma ,y\right)\hfill & \propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}\{{\displaystyle {\int }_0^1\left[{\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}\right]}\hfill \\ \hfill & \times \frac{{\rho }_i^{\varphi \frac{1-\delta }{\delta }-1}{\left(1-{\rho }_i\right)}^{\left(1-\varphi \right)\frac{1-\delta }{\delta }-1}}{B\left(\varphi \frac{1-\delta }{\delta }\mathrm{,}\left(1-\varphi \right)\frac{1-\delta }{\delta }\right)}\},\mathrm{0\; <mo><</mo>}{\rho }_i\mathrm{<mo><</mo>\; 1,}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}\frac{\delta }{1-\delta }\mathrm{<mo><</mo>}\varphi \mathrm{<mo><</mo>}\frac{1-2\delta }{1-\delta }\mathrm{,\; <mi>\; </mi>\; <mi>\; </mi>\; <mtext>\hspace{0.17em}</mtext>\; 0\; <mo><</mo>}\delta \mathrm{<mo><</mo>}\frac{1}{3}.\hfill \end{array}$$

Again, we apply the middle Riemann sum method to integrate out all ${\rho }_i,i=1,\dots ,l$ and compute the joint conditional posterior distribution of $\left(\varphi ,\delta \right),$

$$p\left(\varphi ,\delta |\mu ,\theta ,\gamma ,y\right)\propto {\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}\left[\underset{G\to \infty }{\lim }{\displaystyle \sum _{v=1}^G{h}_i\left(\frac{a_{v-1}+a_v}{2}\right)}\left\{F_2\left(a_{v-1}\right)-F_2\left(a_v\right)\right\}\right],$$

where

$$h_i\left({\rho }_i\right)={\displaystyle \prod _{j=1}^{m_i}\frac{B\left(s_{ij}+{\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}{B\left({\mu }_i\frac{1-{\rho }_i}{{\rho }_i}\mathrm{,}\left(1-{\mu }_i\right)\frac{1-{\rho }_i}{{\rho }_i}\right)}}$$

and $F_2(\cdot )$ is the cdf corresponding to $f_2(\cdot )$ which is a density function of $\text{Beta}\left(\varphi \frac{1-\delta }{\delta },\left(1-\varphi \right)\frac{1-\delta }{\delta }\right).$ Using Gaussian quadrature via Legendre orthogonal polynomials, we can integrate out $\varphi$ and obtain the conditional posterior density of $\delta ,$

$$p\left(\delta |\mu ,\theta ,\gamma ,y\right)\approx {\displaystyle \sum _{g=1}^G{{\omega }^{}}_g}\left\{{\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle {\int }_0^1{\pi }_2\left({\rho }_i,{x^{}}_g,\delta |{\mu }_i,\theta ,\gamma ,y\right)d{\rho }_i}\right\},$$

where $\left\{{{\omega }^{}}_g\right\}$ are the weights and $\left\{{x^{}}_g\right\}$ are roots of the Legendre polynomial with the interval $\left[\frac{\delta }{1-\delta },\frac{1-2\delta }{1-\delta }\right].$

Then, we use the univariate grid method in order to draw samples of the posterior density of $\delta$ conditional on $\mu ,\theta \mathrm{,}\gamma$ and $y.$ Therefore, the conditional posterior density of $\varphi$ can be represented as

$$p\left(\varphi |\delta ,\mu ,\theta ,\gamma ,y\right)\approx {\displaystyle \sum _{g=1}^G{{\omega }^{}}_g}\left\{{\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle {\int }_0^1{\pi }_2\left({\rho }_i,\varphi |\delta ,{\mu }_i,\theta ,\gamma ,y\right)d{\rho }_i}\right\},$$

and we can get samples of $\theta$ by using the univariate grid sampler again. Finally, conditional on $\left(\varphi ,\delta \right),\rho$ can be drawn from $p\left(\rho |\mu ,\theta ,\gamma ,\varphi ,\delta ,y\right),$ where we also use the univariate grid method.

This algorithm samples ${\pi }_1\left(\mu ,\theta ,\gamma |\rho ,\varphi ,\delta ,y\right)$ by first drawing an iterate from ${\pi }_1\left(\gamma |\rho ,\varphi ,\delta ,y\right),$ an iterate from ${\pi }_1\left(\theta |\gamma ,\rho ,\varphi ,\delta ,y\right)$ and then an iterate from ${\pi }_1\left(\mu |\theta ,\gamma ,\rho ,\varphi ,\delta ,y\right).$ Then, it samples ${\pi }_2\left(\rho ,\varphi ,\delta |\mu ,\theta ,\gamma ,y\right)$ by first drawing an iterate from ${\pi }_2\left(\delta |\mu ,\theta ,\gamma ,y\right),$ an iterate from ${\pi }_2\left(\varphi |\delta ,\mu ,\theta ,\gamma ,y\right)$ and then an iterate from ${\pi }_2\left(\rho |\varphi ,\delta ,\mu ,\theta ,\gamma ,y\right).$ The entire procedure continues until convergence. It is like using a Gibbs sampler with two conditional posterior densities which is, in fact, the blocked Gibbs sampler. The construction of the blocked Gibbs sampler is very efficient and it is one of our key contributions in this paper. In fact, we might call the blocked Gibbs sampler the blocked griddy Gibbs sampler (Ritter and Tanner 1992).

We have monitored the convergence of the blocked Gibbs sampler using trace plots, autocorrelation plots and Geweke test of stationarity. The trace plot, iterates versus time, gives information about how long a burn-in period is required to remove the effect of initial values. The autocorrelation plots display dependence in the chain, and thus, in the plots high correlations between long lags indicate a poor mixing chain. The Geweke test compares the means from the early and latter part of the Markov chain by using a $z-$ score statistic, where the null hypothesis is that the chain is stationary; the $p-$ values are all larger than 0.10. We have used the trace plots, autocorrelation plots, and Geweke test for each parameter to study convergence of each run of the blocked Gibbs sampler. For our data, we draw 2,000 samples and burn in 1,000 in order to obtain a sample of 1,000 iterates for inference. This burn-in period, which is based on the trace plots and Geweke test, is long enough to get random samples. The correlations are all nonsignificant, and interestingly, we do not have to thin the iterates. Also, Geweke test demonstrates stationarity of our sampler. Thus, we have a highly efficient blocked Gibbs sampler. The procedure takes a few minutes on R. We have applied the same procedure in our simulation study.

• Previous
• Article main page
• Next