Bayesian predictive inference of a proportion under a two-fold small area model with heterogeneous correlations
Section 3. Numerical study and comparisons

• Previous
• Next

In this section, we perform empirical studies to assess the performance of the HeC model that we compare with the HoC model. In Section 3.1, we discuss an illustrative example and, in Section 3.2, we present a simulation study.

3.1 An illustrative example

We use data from the Third Grade US population; see Nandram (2015) for a brief discussion of these data. The dataset, collected in 1999, consists of 2,477 students who participated in the Third International Mathematics and Science Study (TIMSS). Foy, Rust, and Schleicher (1996) described the probability proportional to size (PPS) systematic sampling design used in TIMSS data collection and Caslyn, Gonzales and Frase (1999) gave highlights from TIMSS. Areas are formed crossing four regions (Northeast, South, Central and West) and three communities of the US (village or rural area, outskirts of a town or city and close to the center of a town or city). Thus, there are twelve areas. The binary variable is whether a student’s mathematics score is below average. Clusters are schools while units within the clusters are the students.

To assess the quality of the Bayesian predictive inference, as suggested by a referee, Nandram (2015) took a half sample of the original data, which he called a synthetic sample. The original sample was used as the population, and the half sample was used for analysis, thereby providing a method to assess the predictive power of the models in Nandram (2015). In the current paper, as suggested by a referee, we do not use a half sample and we use the original dataset available to us; see Table 3.1 for the entire dataset which we analyze in this paper. The predictive power of the HeC model is assessed mainly through the simulation study.

Unfortunately, as in many complex surveys, the sample fractions are unknown to secondary data analysts. However, typically for many of these complex surveys, the sample fractions are relatively small. For the TIMSS data we assume that the dataset is a 5% sample of the population. For example, if there are four sampled schools for an area, say $i^{\text{th}}$ area $\left(i\mathrm{<mtext>\hspace{0.17em}</mtext>\; <mo>=</mo>\; <mtext>\hspace{0.17em}</mtext>\; 1,}\dots \mathrm{,}l\right),$ the total number of clusters, $M_i$ is assumed to be 80. If there are 17 observed students within a sampled school, say $j^{\text{th}}$ school, the total number of students, $N_{ij}\left(j\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{m}_i\right)$ is assumed to be 340. For the nonsampled schools, $N_{ij}\left(j\mathrm{<mo>=</mo>}{m}_i+\mathrm{1,}\dots \mathrm{,}{M}_i\right)$ is assumed to be the average of the total number of students within the sampled schools for each area. Moreover, there are many schools in which all or many students were either below or above average. In other words, this dataset is far sparse, thereby making direct estimation difficult.

We perform three goodness-of-fit procedures, the deviance information criterion (DIC), the Bayesian posterior predictive $p-$ value (BPP) and the log pseudo marginal likelihood (LPML), which is a measure based on the same cross-validation (leave-one-out) procedure. We can assess the overall fit of the models with these procedures.

In the HeC model, $s_{ij}|p_{ij}\stackrel{\text{ind}}{\sim }\text{Binomial}\left(n_{ij}\mathrm{,}{p}_{ij}\right)\mathrm{,\; <mtext>\hspace{0.17em}</mtext>\; <mi>\; </mi>}{p}_{ij}\stackrel{\text{ind}}{\sim }\text{Beta}\left({\mu }_i\left(1-{\rho }_i\right)/{\rho }_i\mathrm{,}\left(1-{\mu }_i\right)\left(1-{\rho }_i\right)/{\rho }_i\right).$ Thus, by integrating out the $p_{ij},$ we can obtain the following beta-binomial probability mass function, 

$$f\left(s|\mu \mathrm{,}\rho \right)\mathrm{<mo>=</mo>}{\displaystyle \prod _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle \prod _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\left(\begin{array}{l}\begin{array}{l}{n}_{ij}\hfill \\ s_{ij}\hfill \end{array}\hfill \end{array}\right)\frac{B\left(s_{ij}+{\mu }_i\left(1-{\rho }_i\right)/{\rho }_i\mathrm{,}{n}_{ij}-s_{ij}+\left(1-{\mu }_i\right)\left(1-{\rho }_i\right)/{\rho }_i\right)}{B\left({\mu }_i\left(1-{\rho }_i\right)/{\rho }_i\mathrm{,}\left(1-{\mu }_i\right)\left(1-{\rho }_i\right)/{\rho }_i\right)}.$$

It is also true that $E\left(s_{ij}|{\mu }_i\mathrm{,}{\rho }_i\right)\mathrm{<mo>=</mo>}{n}_{ij}{\mu }_i$ and $\text{Var}\left(s_{ij}|{\mu }_i\mathrm{,}{\rho }_i\right)\mathrm{<mo>=</mo>}{n}_{ij}\left\{1+\left(n_{ij}-1\right){\rho }_i\right\}{\mu }_i\left(1-{\mu }_i\right).$

Let ${\mu }_i^{\left(h\right)}$ and ${\rho }_i^{\left(h\right)}$ $\left(i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}h\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}H\right)$ denote the iterates from the blocked Gibbs sampler. Let ${\overline{\mu }}_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{h\mathrm{<mo>=</mo>\; 1}}^H{\mu }_i^{\left(h\right)}/H}\left(i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\right)$ and ${\overline{\rho }}_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{h\mathrm{<mo>=</mo>\; 1}}^H{\rho }_i^{\left(h\right)}/H}.$ Letting $D\left(\overline{\mu },\overline{\rho }\right)=-2\log \left\{p\left(s|\overline{\mu },\overline{\rho }\right)\right\}$ and $\overline{D}\mathrm{<mo>=</mo>}-2{\displaystyle {\sum }_{h\mathrm{<mo>=</mo>\; 1}}^H\text{log}\left\{p\left(s|{\mu }^{\left(h\right)},{\rho }^{\left(h\right)}\right)\right\}/H},$ deviance information criterion is given by

$$\text{DIC}\mathrm{<mo>=</mo>\; 2}\overline{D}-D\left(\overline{\mu },\overline{\rho }\right)\mathrm{.}$$

Models with smaller DIC are more preferred over those with larger DIC. However, since DIC tends to select over-fitted models, Nandram (2015) described the Bayesian predictive $p-$ values as a backup. For the HeC model, the discrepancy function is

$$T\left(s;\mu ,\rho \right)={\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle \sum _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}}\frac{{\left\{s_{ij}-E\left(s_{ij}|{\mu }_i,{\rho }_i\right)\right\}}^2}{\text{Var}\left(s_{ij}|{\mu }_i,{\rho }_i\right)}.$$

Let $s^{\left(\text{rep}\right)}$ denote repeated (rep) samples from the posterior predictive distribution of $s.$ Then the BPP is $p\left\{T\left(s^{\left(\text{rep}\right)};\mu ,\rho \right)\ge T\left(s^{\left(\text{obs}\right)};\mu ,\rho \right)|s\right\},$ which is calculated over its corresponding iterates $\left({\mu }^{\left(h\right)}\text{},{\rho }^{\left(h\right)}\right),h=1,\dots ,H.$ If the value of this probability is close to 0 or 1, it indicates poor fit of the model. In fact, models with BPPs in (0.05, 0.95) are considered reasonable.

In addition to these quantities, we can evaluate the goodness-of-fit of models with another measure, the LPML which is a summary statistic of the conditional predictive ordinate (CPO) values, and it is based on a cross validation. Unlike the DIC, larger values of LPML indicate better fitting models (e.g., Geisser and Eddy 1979).

For the HeC model, the CPO can be estimated by

$${\widehat{\text{CPO}}}_{ij}={\left[\frac{1}{H}{\displaystyle \sum _{h=1}^H\frac{1}{f\left(s_{ij}|p_{ij}^{\left(h\right)}\right)}}\right]}^{-1}\text{},j=1,\dots ,m_i,i=1,\dots ,l,$$

where $p_{ij}^{\left(h\right)}$ is the samples from $p_{ij}|s_{ij}\mathrm{,}{\mu }_i\mathrm{,}{\rho }_i$ and $s_{ij}|p_{ij}\stackrel{\text{iid}}{\sim }\text{Binomial}\left(n_{ij}\mathrm{,}{p}_{ij}\right).$ Note that for each $\left(i\mathrm{,}j\right),$ ${\widehat{\text{CPO}}}_{ij}$ is the harmonic mean of the likelihoods $f\left(s_{ij}|p_{ij}^{\left(h\right)}\right)\mathrm{,\; <mtext>\hspace{0.17em}</mtext>\; <mi>\; </mi>}h\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}H.$ Then, the LPML is

$$\text{LPML}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>\; 1}}^l}{\displaystyle \sum _{j\mathrm{<mo>=</mo>\; 1}}^{m_i}\log \left({\widehat{\text{CPO}}}_{ij}\right)}\mathrm{.}$$

These three model evaluation measures have similar forms under the HoC model. For the HoC (HeC) model, DIC = 774.421 (773.173), BPP = 0.349 (0.408), LPML = -352.064 (-346.171), thereby indicating that the HeC model gives a better fit. At a finer level, we also looked at the individual CPO values from the two models for each school. In Figure 3.1 we compare the CPOs from the HeC and the HoC models, and we found that generally CPO values for the HeC model are higher that those of HoC model. In fact, under the HoC (HeC) model we found that the percent of the CPOs less than 0.025 is 3.70% (2.96%) and percent of the CPOs less than 0.014 is 0.74% (0.00%). These results do not show any indication of serious departure from model assumptions; see Ntzoufras (2009). Therefore, these measures give prima facie evidence that the HeC model fits the TIMSS data somewhat better than the HoC model.

Now, consider inference about $\theta$ and $\gamma .$ First, consider $\theta .$ Under the HoC model, the posterior mean (PM) is 0.519, posterior standard deviation (PSD) is 0.068 and 95% credible interval (Cre) is (0.390, 0.639). Under the HeC model PM = 0.515, PSD = 0.065 and 95% Cre is (0.383, 0.639). Second, consider $\gamma .$ Under the HoC model PM = 0.207, PSD = 0.011, and 95% Cre is (0.190, 0.224). Under the HeC model $\gamma$ PM = 0.208, PSD = 0.011 and 95% Cre is (0.190, 0.225). Thus, it is good that inference about $\theta$ and $\gamma$ are very close for the two competitors (HoC and HeC models).

In Table 3.2, we present posterior inference about the finite population proportions for mathematics scores by areas. There are differences between the posterior means under the HoC and HeC models. Most of them are small but there are a few large differences. For NC, SR and CR, we have 0.560 (0.543), 0.568 (0.584) and 0.465 (0.445) under the HoC (HeC) model, respectively. The posterior standard deviations are also close but there are a few moderately large differences (e.g., for NR we have 0.113 under the HoC model and 0.077 under the HeC model). These differences are reflected in the credible and highest posterior density (HPD) intervals.

Table 3.3 shows summaries of the PM, PSD and 95% HPD for intracluster correlations under the HeC model. We can see that the intracluster correlations vary over the areas. The largest estimate is 0.337 for NC and the smallest one is 0.073 for SR. Both areas have a few large difference between the posterior means under the HoC and HeC models. The 95% HPD interval for the common correlation in the HoC model is (0.160, 0.260) and this interval is contained by all the intervals except for NR, NC, SR and WC. Thus, it is reasonable to study the HeC model.

In Figure 3.2, we compare the posterior densities of the intracluster correlations (twelve correlations) from the HeC model and the HoC model (one correlation). The distributions under the HeC model are more variable and are mostly to the left or right of those under the HoC model with not much overlap for some areas (e.g., NR, NC and SR).

In Figures 3.3, 3.4 and 3.5, we compare the posterior density plots of the finite population proportions for the mathematics score and all areas for the two models. There are noticeable differences between the HoC and HeC models (e.g., areas NR, NC, SR, CR and WC).

3.2 Simulation study

In order to further assess the performance of the HeC model and to compare it to the HoC model, we perform a simulation study. Here we use two factors, each at three levels, to get nine design points.

We have set 100 as the number of clusters (schools) in each area and 15 as the number of individuals (students) within each cluster. In other words, we take $N_{ij}\mathrm{<mo>=</mo>\; 15,}j\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{M}_i\mathrm{,\; <mtext>\hspace{0.17em}</mtext>\; <mi>\; </mi>}{M}_i\mathrm{<mo>=</mo>\; 100,\; <mtext>\hspace{0.17em}</mtext>\; <mi>\; </mi>}i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l$ where $l\mathrm{<mo>=</mo>\; 12.}$ Let $a$ denote a vector of posterior means and $b$ denote the vector of posterior standard deviations corresponding to the ${\mu }_i$ or the ${\rho }_i.$ Specifically, for the ${\rho }_i$ we use $a_1$ and $b_1$ and for the ${\mu }_i$ we use $a_2$ and $b_2.$ When we simulate data from the HeC model, the levels of the ${\rho }_i$ are $\left(1:a_1-0.5b_1;2:a_1;3:a_1+0.5b_1\right)$ and the levels of the ${\mu }_i$ are $\left(1:a_2-0.5b_2;2:a_2;3:a_2+0.5b_2\right).$ For the twelve areas $a_1$ takes values 0.09, 0.19, 0.32, 0.08, 0.22, 0.18, 0.15, 0.22, 0.23, 0.17, 0.18, 0.30; $b_1$ 0.08, 0.09, 0.08, 0.06, 0.07, 0.08, 0.13, 0.09, 0.07, 0.09, 0.07, 0.06; $a_2$ 0.53, 0.37, 0.54, 0.58, 0.54, 0.65, 0.46, 0.44, 0.55, 0.46, 0.52, 0.66; and $b_2$ 0.08, 0.08, 0.08, 0.06, 0.06, 0.06, 0.12, 0.09, 0.07, 0.08, 0.06, 0.05.

We also take a simple random sample of five clusters among the 100 population clusters, and a simple random sample of ten individuals from each sampled cluster (i.e., $m_i\mathrm{<mo>=</mo>\; 5}$ and $n_{ij}\mathrm{<mo>=</mo>\; 10}).$ These numbers are much smaller than those of data used in Section 3.1, which makes inference a little more challenging (Nandram 2015). Note that the dataset has about 7% of the sampled clusters where all students were either below or above average. We call this quantity the percent of sparseness. The setting of this simulation study also leads to even sparser data. For nine design points, all the average percents of sparseness are greater than 7% and most are around 10%. Figure 3.6 shows the histograms of sparseness percents for each design point.

We consider two scenarios. In the first scenario, we generate data from the HeC model and fit both models, and in the second scenario, we generate data from the HoC model and fit both models. When data are simulated from the HeC model, we have nine design points $\mathrm{<mo\; stretchy="false">[</mo>\; <mo\; stretchy="false">(</mo>\; 1,1\; <mo\; stretchy="false">)</mo>\; ,\; <mo\; stretchy="false">(</mo>\; 1,2\; <mo\; stretchy="false">)</mo>\; ,\; <mo\; stretchy="false">(</mo>\; 1,3\; <mo\; stretchy="false">)</mo>\; ,}\dots \mathrm{,\; <mo\; stretchy="false">(</mo>\; 3,1\; <mo\; stretchy="false">)</mo>\; ,\; <mo\; stretchy="false">(</mo>\; 3,2\; <mo\; stretchy="false">)</mo>\; ,\; <mo\; stretchy="false">(</mo>\; 3,3\; <mo\; stretchy="false">)</mo>\; <mo\; stretchy="false">]</mo>},$ the first factor corresponds to the ${\rho }_i.$ When we simulate data from the HoC model, we have three design points $\left(1:a_2-0.5b_2;2:a_2;3:a_2+0.5b_2\right)$ for the three levels for the ${\mu }_i;$ $\rho$ is kept fixed at its posterior mean.

In the first scenario, at each design point we simulate binary data from the HeC model,

$$\begin{array}{ll}{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{,}\hfill \\ y_{ijk}|p_{ij}\hfill & \stackrel{\text{ind}}{\sim }\text{Bernoulli}\left(p_{ij}\right)\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{N}_{ij}\mathrm{.}\hfill \end{array}$$

So we have the true values of $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}/{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>\; 1}}^{M_i}{N}_{ij}}}$ for $i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l.$ We take 1,000 samples at each of the nine design points. For each sample we perform the blocked griddy Gibbs sampler in the same manner as for the data.

Like Nandram (2015), we calculate ${\text{AB}}_{ih}\mathrm{<mo>=</mo>}\left|{\text{PM}}_{ih}-P_{ih}\right|\mathrm{,}{\text{RAB}}_{ih}\mathrm{<mo>=</mo>}{\text{AB}}_{ih}/P_{ih}$ and ${\text{RPMSE}}_{ih}\mathrm{<mo>=</mo>}\sqrt{{\text{PSD}}_{ih}^2+{\text{AB}}_{ih}^2}$ to study the frequentist properties of our procedure $\left(i\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}l\mathrm{,}h\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}\text{1,000}\right).$ We also obtain the 95% credible interval and HPD interval for each of the 1,000 simulated runs, and we study the width $W_{ih}$ and the credible incidence $I_{ih}.$ If the 95% credible (or HPD) interval of $h^{\text{th}}$ run contains the true value $P_i,$ $I_{ih}$ is equal to one, otherwise it is equal to zero. Thus, the estimated probability content of the 95% credible interval for the $i^{\text{th}}$ area is $C_i\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{h\mathrm{<mo>=</mo>\; 1}}^{\text{1,000}}{I}_{ih}/\text{1,000}}.$

Table 3.4 shows comparison of the HoC and HeC models. Under the HeC model the coverages are much higher than those under the HoC model. Note that the coverages of HPD intervals for the HeC model are much closer to the nominal value of 95% and they are conservative. However, the 95% credible and HPD intervals are wider than those from the HoC model. These effects are much larger as $\rho$ becomes larger. All measures AB, RAB and RPMSE under the HeC model are smaller than those under the HoC model. Thus, based on these measures, the HeC model is preferred over the HoC model.

In Table 3.5 we compare summaries of DIC, BPP and LPML. All the DICs under the HeC model are smaller than the corresponding ones under the HoC model and all the LPMLs under the HeC model are larger than those under the HoC model. Under the HoC model, all the BPPs vary in (0.06, 0.09) but under the HeC model they vary in (0.2, 0.4). Again, these measures show that the HeC model is superior to the HoC model.

In a similar manner, for the second scenario we generate binary data from

$$\begin{array}{ll}{p}_{ij}|{\mu }_i\mathrm{,}\rho \hfill & \stackrel{\text{ind}}{\sim }\text{Beta}\left[{\mu }_i\frac{1-\rho }{\rho }\mathrm{,}\left(1-{\mu }_i\right)\frac{1-\rho }{\rho }\right]\mathrm{,}\hfill \\ y_{ijk}|p_{ij}\hfill & \stackrel{\text{ind}}{\sim }\text{Bernoulli}\left(p_{ij}\right)\mathrm{,}k\mathrm{<mo>=</mo>\; 1,}\dots \mathrm{,}{N}_{ij}\mathrm{.}\hfill \end{array}$$

In Table 3.6 we present comparison of the HoC and HeC models. Here AB, RAB and RPMSE are only slightly smaller under the HoC model. The coverages of the credible and HPD intervals under the HeC model are closer to the nominal value of 95%, while those under the HoC model are smaller. Table 3.7 shows summaries of DIC, BPP and LPML. All the DICs under the HeC model are smaller than those under the HeC model, while the BPPs and LPMLs are similar for the two models, with those under the HoC model being slightly better.

Thus, when data actually come from the HeC model, there are some important differences among the two models, with the HeC model being preferred. However, when data actually come from the HoC model, there are minor differences between the two models. Of course, the HeC model (unequal correlations) has more parameters than the HoC model (one correlation).

• Previous
• Article main page
• Next