Bayesian inference for a variance component model using pairwise composite likelihood with survey data
Section 3. Simulation studies

• Previous
• Next

3.1 Simulation design

Using simulation studies we have evaluated the performance of the proposed method, i.e., pairwise composite likelihood with a curvature adjustment, and compared it with using the full likelihood and the pairwise composite likelihood. We used the model in (1.1) to generate our data, i.e., for $i\mathrm{<mo>=</mo>}\mathrm{1,}\dots \mathrm{,}n$ and $j\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,m$ we simulated values of $Y_{ij}$ from

$$Y_{ij}\mathrm{<mo>=</mo>}\theta +u_i+e_{ij}\mathrm{,}(3.1)$$

where $\theta \mathrm{<mo>=</mo>}1,$ $u_i\stackrel{\text{iid}}{{\displaystyle ~}}N\left(\mathrm{0,}{\sigma }_u^2\right),$ and $e_{ij}\stackrel{\text{iid}}{{\displaystyle ~}}N\left(\mathrm{0,}{\sigma }_e^2\right)\mathrm{.}$ This is equivalent to having applied the superpopulation generation and sampling described in the paragraph surrounding (1.1).

Our first study, not included here, considered inference about $\theta$ with known ${\sigma }_u$ and ${\sigma }_e.$ It showed that using the pairwise composite likelihood for inference about $\theta$ badly overstated the precision, and that the curvature adjustment was successful. Thus, we proceeded to a more thorough study, considering inference for both $\theta$ and ${\sigma }_u.$ To simplify we took ${\sigma }_e\mathrm{<mo>=</mo>}\text{0}\text{.5},$ and considered $n\in \left\{\mathrm{20,}40\right\}$ and $m\in \left\{\mathrm{5,}10\right\}.$ For the half-Cauchy prior defined in (2.4) we took $A\in \left\{\mathrm{5,}\mathrm{10,}15\right\}.$ There were 500 replicate data sets for each setting.

We considered three scenarios: (1) ${\sigma }_u\in \left\{\text{0}\text{.1},\text{0}\text{.5}\right\}$ and the half-Cauchy prior on ${\sigma }_u;$ (2) Signal to Noise Ratio, $\text{SNR}\in \left\{\text{0}\text{.25}\mathrm{,}\text{0}\text{.75}\right\}$ and the half-Cauchy prior on ${\sigma }_u,$ where $\text{SNR}\mathrm{<mo>=</mo>}{\sigma }_u^2/\left({\sigma }_u^2+{\sigma }_e^2\right);$ and (3) ${\sigma }_u\in \left\{\text{0}\text{.1}\mathrm{,}\text{0}\text{.5}\right\}$ and a uniform prior on ${\sigma }_u.$ Throughout, we took a uniform prior on $\theta .$

In Section 3.2 we describe the algorithms for the simulation studies.

3.2 Algorithms

As in Sections 2.1 and 2.2 define $y\left(n\right)\mathrm{<mo>=</mo>}\left\{y_1\mathrm{,}\dots \mathrm{,}{y}_n\right\}$ with $y_i\mathrm{<mo>=</mo>}{\left(y_{i1}\mathrm{,}\dots \mathrm{,}{y}_{im}\right)}^{\text{T}}$ and $\overline{y}\mathrm{<mo>=</mo>}{\displaystyle {\sum }_{i\mathrm{<mo>=</mo>1}}^n{\displaystyle {\sum }_{j\mathrm{<mo>=</mo>1}}^m{y}_{ij}/\left(mn\right)}}.$ Further, ${\eta }^{\left(t\right)}$ denotes the value of $\eta$ at the $t^{\text{th}}$ iteration where $\eta \mathrm{<mo>=</mo>}{\left(\theta \mathrm{,}{\sigma }_u\right)}^{\text{T}}.$ The full likelihood is

$$L_{\text{FL}}\left(\theta \mathrm{,}{\sigma }_u\mathrm{<mo>|</mo>}y\left(n\right)\right)\propto {\left|{\Sigma }_m\right|}^{-n/2}\text{exp}\left[-\frac{1}{2}\text{tr}\left({\Sigma }_m^{-1}{S}_0\right)\right]\mathrm{,}(3.2)$$

as in (2.7).

Using (3.2) together with the prior, $\pi \left(\eta \right),$ yields the posterior density,

$$p_{\text{FL}}\left(\eta |y\left(n\right)\right)\propto L_{\text{FL}}\left(\eta |y\left(n\right)\right)\pi \left(\eta \right)\mathrm{.}$$

Sampling $\theta$ and ${\sigma }_u$ is done in three steps:

Step 1. Sample ${\theta }^{\left(t\right)}$ from $p_{\text{FL}}\left(\theta |y\left(n\right)\mathrm{,}{\sigma }_u^{(t-1)}\right)$ where

$$\theta |\left(y\left(n\right)\mathrm{,}{\sigma }_u\right)~N\left(\overline{y},\frac{{\sigma }_e^2+m{\sigma }_u^2}{mn}\right)\mathrm{.}$$

We set the starting value, ${\sigma }_u^{\left(0\right)},$ to be the maximum likelihood estimate of ${\sigma }_u.$

Step 2. Use the Metropolis-Hastings (MH) algorithm to sample ${\sigma }_u^{\left(t\right)}$ from $p_{\text{FL}}\left({\sigma }_u|y\left(n\right)\mathrm{,}{\theta }^{(t)}\right).$ The latter is easily obtained from $p_{\text{FL}}\left(\eta |y\left(n\right)\right).$ Given $s\mathrm{<mo>></mo>}0,$ the candidate ${\sigma }_u,$ labelled ${\sigma }_u^{*},$ is sampled from the jumping distribution, $N\left({\sigma }_u^{\left(t-1\right)}\mathrm{,}{s}^2\right).$ If ${\sigma }_u^{*}<\mathrm{0,}{\sigma }_u^{\left(t\right)}\mathrm{<mo>=</mo>}{\sigma }_u^{\left(t-1\right)}.$ Otherwise, the procedure is standard with accept/reject ratio $p_{\text{FL}}\left({\eta }_{\text{FL}}^{*}|y\left(n\right)\right)/p_{\text{FL}}\left({\eta }_{\text{FL}}^{(t-1)}|y\left(n\right)\right)$ where $\eta {}_{\text{FL}}^{*}={\left({\theta }^{\left(t\right)}\mathrm{,}{\sigma }_u^{*}\right)}^{\text{T}}$ and $\eta {}_{\text{FL}}^{(t-1)}={\left({\theta }^{(t)}\mathrm{,}{\sigma }_u^{(t-1)}\right)}^{\text{T}}.$

Step 3. Repeat Steps 1 and 2 for $K\mathrm{<mo>=</mo>}\text{1,000}$ times with the first 200 samples used as the burn-in.

The pairwise composite likelihood (PL) is

$$L_{\text{PL}}\left(\theta \mathrm{,}{\sigma }_u|y\left(n\right)\right)\propto {\left|\Sigma {}_2\right|}^{nm\left(m-1\right)/4}\text{exp}\left[-\frac{1}{2}\text{tr}\left({\Sigma }_2^{-1}{S}_{0\text{PL}}\right)\right]\mathrm{,}(3.3)$$

as in (2.8).

Using (3.3) together with the chosen prior, $\pi \left(\eta \right),$ yields the posterior density $p_{\text{PL}}\left(\eta |y\left(n\right)\right).$

Sampling $\theta$ and ${\sigma }_u$ is done in three steps:

Step 1. Sample ${\theta }^{\left(t\right)}$ from $p_{\text{PL}}\left(\theta |y\left(n\right)\mathrm{,}{\sigma }_u^{(t-1)}\right)$ where

$$\theta |\left(y\left(n\right)\mathrm{,}{\sigma }_u\right)~N\left(\overline{y},\frac{{\sigma }_e^2+2{\sigma }_u^2}{nm\left(m-1\right)}\right)\mathrm{.}$$

Step 2. Use the Metropolis-Hastings (MH) algorithm to sample ${\sigma }_u^{\left(t\right)}$ from $p_{\text{PL}}\left({\sigma }_u|y\left(n\right)\mathrm{,}{\theta }^{(t)}\right),$ as described in Step 2 above for the FL (substituting PL for FL in all formulas).

Step 3. Repeat Steps 1 and 2 for $K\mathrm{<mo>=</mo>}\text{1,000}$ times with the first 200 samples used as the burn-in.

The final part is to obtain the (curvature) adjusted pairwise composite likelihood (APL), as described in Section 2.3. This derivation, based on the approach of RCD, exploits ${\widehat{\eta }}_{A\text{PL}},$ the estimated posterior means of $\theta$ and ${\sigma }_u.$

Step 1. Given $\left(s_{\theta }\mathrm{,}{s}_{\sigma }\right)$ sample the candidate ${\eta }^{*}\mathrm{<mo>=</mo>}{\left({\theta }^{*}\mathrm{,}{\sigma }_u^{*}\right)}^{\text{T}}$ from the bivariate normal jumping distribution, $N_2\left(\eta {}^{\left(t-1\right)}\mathrm{,}\Sigma \right)$ where $\Sigma \mathrm{<mo>=</mo>}\text{diag}\left(s_{\theta }^2\mathrm{,}{s}_{\sigma }^2\right).$ If ${\sigma }_u^{*}\mathrm{<mo><</mo>0,}{\eta }^{\left(t\right)}\mathrm{<mo>=</mo>}{\eta }^{\left(t-1\right)}.$ Otherwise, go to Step 2.

Step 2. Define $l_{\text{PL}}\left(y\left(n\right)|\theta \mathrm{,}{\sigma }_u\right)$ as the log pairwise composite likelihood obtained by taking the logarithm of (3.3), and $l_{\text{PL}}\left(y_i|\theta \mathrm{,}{\sigma }_u\right)$ as the log pairwise composite likelihood corresponding to the data from cluster $i,$ i.e., $y_i.$

Step 3. Numerically obtain $\widehat{H}\mathrm{<mo>=</mo>}{\nabla }^2{l}_{\text{PL}}\left(y\left(n\right)|{\widehat{\theta }}_{\text{PL}},{\widehat{\sigma }}_{u\text{PL}}\right)$ and

$$\widehat{J}\mathrm{<mo>=</mo>}{\displaystyle \sum _{i\mathrm{<mo>=</mo>1}}^n\left[\nabla l_{\text{PL}}\left(y_i|{\widehat{\theta }}_{\text{PL}}\mathrm{,}{\widehat{\sigma }}_{u\text{PL}}\right){\left\{\nabla l\left(y_i|{\widehat{\theta }}_{\text{PL}}\mathrm{,}{\widehat{\sigma }}_{u\text{PL}}\right)\right\}}^{\text{T}}\right]}\mathrm{,}$$

where ${\widehat{\theta }}_{\text{PL}}$ and ${\widehat{\sigma }}_{u\text{PL}}$ are the estimated posterior means of $\theta$ and ${\sigma }_u.$

Step 4. Based on the approach of RCD, and using the singular value decomposition, we write $\widehat{H}\mathrm{<mo>=</mo>}{M}^{\text{T}}M$ and $\widehat{H}{\widehat{J}}^{-1}\widehat{H}\mathrm{<mo>=</mo>}{M}_A^{\text{T}}{M}_A$ for some matrices $M$ and $M_A.$ Then define $C\mathrm{<mo>=</mo>}{M}^{-1}{M}_A.$ In our case, $C$ is a $2\times 2$ matrix.

Step 5. From RCD the adjusted log pairwise composite likelihood, $l_{\text{APL}},$ is

$$l_{\text{APL}}\left(y\left(n\right)|\eta \right)\mathrm{<mo>=</mo>}{l}_{\text{PL}}\left(y\left(n\right)|{\eta }^{*}\right)$$

where

$${\eta }^{*}\mathrm{<mo>=</mo>}{\widehat{\eta }}_{\text{PL}}+C\left(\eta -{\widehat{\eta }}_{\text{PL}}\right).$$

Step 6. Define the adjusted pairwise posterior density as

$$p_{\text{APL}}\left(\eta |y\left(n\right)\right)\propto L_{\text{APL}}\left(y\left(n\right)|\eta \right)\pi \left(\theta \mathrm{,}{\sigma }_u\right)$$

where $L_{\text{APL}}\left(y\left(n\right)|\eta \right)\mathrm{<mo>=</mo>}\text{exp}\left(l_{\text{APL}}\left(y\left(n\right)|\eta \right)\right),$ the latter defined in Step 5.

Using the candidate value, ${\eta }^{*},$ from Step 1 define the adjusted candidate value ${\eta }_c^{*}\mathrm{<mo>=</mo>}{\widehat{\eta }}_{\text{PL}}+C\left({\eta }^{*}-{\widehat{\eta }}_{\text{PL}}\right).$ Then the accept/reject ratio is

$$p_{\text{APL}}\left({\text{η}}_c^{*}|y\left(n\right)\right)/p_{\text{APL}}\left({\eta }^{\left(t\right)}|y\left(n\right)\right).$$

The remaining steps are the standard ones for the Metropolis-Hastings algorithm.

3.3 Results from simulations

For each method (FL, PL, APL), each design parameter $\left(m\mathrm{,}n\right)$ and each prior distribution we summarized the simulation results using (a) the credible interval coverage rate in repeated sampling, and (b) the averages of the 0.025, 0.25, 0.50, 0.75 and 0.975 points of the posterior distributions of $\theta$ and ${\sigma }_u.$

There are also graphical summaries, i.e., averaged posterior density estimates for each of the posterior distributions, i.e., $p_{\text{FL}}\left(\eta |y\left(n\right)\right),p_{\text{PL}}\left(\eta |y\left(n\right)\right)$ and $p_{\text{APL}}\left(\eta |y\left(n\right)\right).$ First, consider an interval, say, $\left[a\mathrm{,}b\right],$ that supports most of the mass (e.g., 95%) of the posterior densities. Then divide it into $M\mathrm{<mo>=</mo>}50$ equally spaced subintervals with the cut points $a\mathrm{<mo>=</mo>}{c}_0<c_1<\dots <c_{M-1}<c_M\mathrm{<mo>=</mo>}b.$ For $t\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,T,$ let ${\widehat{f}}_P^{\left(t\right)}\left(\mathrm{.}\right)$ denote the estimate of the posterior density $f_P\left(\mathrm{.}\right),$ derived from the $t^{\text{th}}$ simulation, where $P$ stands for FL, PL or APL, and $T$ is the number of simulations. Next define, for $r\mathrm{<mo>=</mo>}\mathrm{1,}\dots ,M,$

$${\widehat{f}}_P\left(c_r\right)\mathrm{<mo>=</mo>}\frac{1}{T}{\displaystyle \sum _{t\mathrm{<mo>=</mo>}1}^T{\widehat{f}}_P^{\left(t\right)}\left(c_r\right)}\mathrm{.}$$

Then a curve connecting the points $\left\{c_r\mathrm{,}{\widehat{f}}_P\left(c_r\right)\right\}$ for $a\mathrm{<mo>=</mo>}{c}_0<c_1<\dots <c_M\mathrm{<mo>=</mo>}b,$ is taken as the averaged posterior density estimate for $f_P\left(\mathrm{.}\right).$

Table 3.1 presents the coverage rates for $\theta$ and ${\sigma }_u$ for $A\mathrm{<mo>=</mo>}\mathrm{15,}$ $n\in \left\{\mathrm{20,}40\right\},$ $m\in \left\{\mathrm{5,}10\right\},$ and ${\sigma }_u\in \left\{\text{0}\text{.1}\mathrm{,}\text{0}\text{.289}\mathrm{,}\text{0}\text{.5}\mathrm{,}\text{0}\text{.866}\right\}\mathrm{.}$ Figure 3.1 has the average posterior density estimates for $\theta$ and ${\sigma }_u$ for $A\mathrm{<mo>=</mo>}\mathrm{15,}{\sigma }_u\in \left\{\text{0}\text{.1}\mathrm{,}\text{0}\text{.5}\right\}\mathrm{,}n\mathrm{<mo>=</mo>}40,$ and $m\mathrm{<mo>=</mo>}10.$ In both Table 3.1 and Figure 3.1 the summaries are given for the full likelihood (FL), pairwise composite likelihood (CL), and adjusted pairwise composite likelihood (APL).

The following summary includes the results for only the half-Cauchy prior with $A\in \left\{\mathrm{5,}\mathrm{10,}15\right\},$ $m\in \left\{\mathrm{5,}10\right\},$ $n\in \left\{\mathrm{20,}40\right\},$ and ${\sigma }_u\in \left\{\text{0}\text{.1,}\text{0}\text{.289,}\text{0}\text{.5,}\text{0}\text{.866}\right\},$ the second and fourth values corresponding to SNR = 0.25 and SNR = 0.75, respectively. The results are similar for the three choices of $A,$ and for the uniform prior.

Without any adjustment the coverages of PL differ substantially from the nominal 0.95. For example (Table 3.1), for $A\mathrm{<mo>=</mo>}15,$ $n\mathrm{<mo>=</mo>}40,$ $m\mathrm{<mo>=</mo>}10,$ and ${\sigma }_u\mathrm{<mo>=</mo>}\text{0}\text{.5},$ the coverage for $\theta$ is less than 0.30. Considering all values of the design parameters, the largest coverage is 0.70. In most cases, the coverage for $\theta$ is much less than 0.70.

With the curvature adjustment the coverage for $\theta$ is excellent. Of the 48 cases (three choices of $A,$ two choices of $m,$ two choices of $n,$ four choices of ${\sigma }_u),$ thirteen had coverage between 0.93 and 0.95, twenty-two between 0.92 and 0.93, eleven between 0.91 and 0.92, and two below 0.91, with the latter for ${\sigma }_u\mathrm{<mo>=</mo>}\text{0}\text{.1},$ $n\mathrm{<mo>=</mo>}20,$ $m\mathrm{<mo>=</mo>}10,$ and $A=5$ and 15.

With the curvature adjustment the coverage for ${\sigma }_u$ varies considerably, but there is, in almost all cases, a very large improvement in coverage relative to using the uncorrected pairwise composite likelihood.

The plots (Figure 3.1) show that for $\theta$ the posterior distribution corresponding to the adjusted likelihood is very close to the posterior distribution using the full likelihood. For ${\sigma }_u$ there are differences between the posterior distributions corresponding to the full and adjusted likelihoods, most notably a shift to smaller values for the latter.

To investigate the effects of increasing $m$ and $n,$ consider the difference $\delta \mathrm{<mo>=</mo>}{C}_{\text{FL}}-C_{\text{APL}}$ where $C$ denotes coverage and FL and APL refer to the corresponding posterior distributions.

Overall with all $m\mathrm{,}n\mathrm{,}A,$ and ${\sigma }_u,$ for $\theta ,$ $\delta$ decreases as $n$ increases. For the larger values of ${\sigma }_u,$ $\delta$ decreases as $m$ increases, while for the smaller values of ${\sigma }_u,$ $\delta$ tends to increase as $m$ increases. Overall, for ${\sigma }_u,$ $\delta$ decreases as $n$ increases except in the case ${\sigma }_u\mathrm{<mo>=</mo>}\text{0}\text{.1},$ while $\delta$ increases as $m$ increases.

The reason for the deterioration of the adjustment as $m$ increases might be that the number of pairs per cluster is $m\left(m-1\right)/2$ and increases more rapidly, so that the pairwise likelihood quickly becomes more concentrated around its mode; the curvature adjustment may not suffice to compensate for a change in shape of the log pairwise composite likelihood, e.g., an increase in kurtosis.

Table 3.2 presents one-sided non-coverage rates of the 95% Credible Intervals for $\theta$ and ${\sigma }_u$ with $A\mathrm{<mo>=</mo>}15,$ $n\in \left\{\mathrm{20,}40\right\},$ $m\in \left\{\mathrm{5,}10\right\},$ and ${\sigma }_u\in \left\{\text{0}\text{.1,}\text{0}\text{.289,}\text{0}\text{.5,}\text{0}\text{.866}\right\}.$ We observe the following:

• For $\theta ,$ the non-coverage for full likelihood intervals appears symmetric. The adjusted pairwise likelihood has undercoverage for $\theta ,$ and except when ${\sigma }_u$ is 0.1 the non-coverage is symmetric. A dependence of the coverage on $m$ is seen only in the ${\sigma }_u\mathrm{<mo>=</mo>}\text{0}\text{.1}$ case.
• For ${\sigma }_u,$ the full likelihood interval has non-coverage that is close to nominal and not very skewed, except in the case when ${\sigma }_u\mathrm{<mo>=</mo>}\text{0}\text{.1},$ where there is marked over-coverage. For ${\sigma }_u>\text{0}\text{.1}$ and $m\mathrm{<mo>=</mo>}5,$ coverage improves as $n$ moves from 20 to 40, but for ${\sigma }_u>\text{0}\text{.1}$ and $m\mathrm{<mo>=</mo>}10,$ there is little difference in coverage for the two values of $n.$
• For ${\sigma }_u,$ the adjusted pairwise likelihood has asymmetric non-coverage. Except in the case of ${\sigma }_u\mathrm{<mo>=</mo>}\text{0}\text{.1},$ the magnitude of the non-coverage tends to be similar on the left to that of the full likelihood, but much greater on the right, and the coverage improves as $n$ moves from 20 to 40.

Remembering that the adjusted log pairwise likelihood is not explicitly being constructed to approximate the log full likelihood, it does appear in Figure 3.1 that the adjusted log pairwise likelihood falls more quickly in the tails.

We also tried centering the curvature adjustment at the log pairwise posterior mode rather than the log pairwise posterior mean, and found 

• Previous
• Article main page
• Next