Parsimonious and Efficient Likelihood Composition by Gibbs Sampling

Let $\{\mathcal{A}_{1},\dots,\mathcal{A}_{M}\}$ be a set of marginal or conditional sample sub-spaces associated with probability density functions (pdfs) $f_{m}(\boldsymbol{x}\in\mathcal{A}_{m}|\boldsymbol{\theta})$. See ? for interpretation and examples of $\mathcal{A}_{m}$. Given independent $d$-dimensional observations $\boldsymbol{X}^{(1)},\dots,\boldsymbol{X}^{(n)}\sim f(\boldsymbol{x}|\boldsymbol{\theta}_{0})$, $\boldsymbol{\theta}_{0}\in\boldsymbol{\Theta}\subseteq\mathbb{R}^{p}$, $p\geq 1$, we define the composite log-likelihood function:

$\displaystyle\ell_{cl}(\boldsymbol{\theta},\boldsymbol{\omega})=\sum_{m=1}^{M}\omega_{m}\ell_{m}(\boldsymbol{\theta}),$

(2)

where $\boldsymbol{\omega}=(\omega_{1},\dots,\omega_{M})^{T}\in\boldsymbol{\Omega}=\{0,1\}^{M}$, and $\ell_{m}(\boldsymbol{\theta})$ is the partial log-likelihood

$\displaystyle\ell_{m}(\boldsymbol{\theta})=\sum_{i=1}^{n}\log f_{m}(\boldsymbol{X}^{(i)}\in\mathcal{A}_{m}|\boldsymbol{\theta}).$

(3)

Each partial likelihood object (sub-likelihood) $\ell_{m}(\boldsymbol{\theta})$ is allowed to be selected or not, depending on whether $\omega_{m}$ is $1$ or $0$. We aim to approximate the unknown complete log-likelihood function $\ell_{mle}(\boldsymbol{\theta})=\log L_{mle}(\boldsymbol{\theta})$ by selecting a few – yet the most informative – sub-likelihood objects from the $M$ available objects, where $M$ is allowed to be larger than $n$ and $p$. Given the composition rule $\boldsymbol{\omega}\in\boldsymbol{\Omega}$, the maximum composite likelihood estimator (McLE), denoted by $\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})$, is defined by the solution of the following system of estimating equations:

$\displaystyle\mathbf{0}=\sum_{i=1}^{n}\boldsymbol{U}^{(i)}(\boldsymbol{\theta},\boldsymbol{\omega}):=\sum_{i=1}^{n}\sum_{m=1}^{M}\omega_{m}\boldsymbol{U}^{(i)}_{m}(\boldsymbol{\theta}),$

(4)

where $\boldsymbol{U}^{(i)}_{m}(\boldsymbol{\theta})=\nabla_{\boldsymbol{\theta}}\log f_{m}(\boldsymbol{X}^{(i)}\in\mathcal{A}_{m}|\boldsymbol{\theta})$, $m=1,\dots,M$, $i=1,\dots,n$, are unbiased partial scores functions. Under standard regularity conditions on the sub-likelihoods (?), $\sqrt{n}(\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})-\boldsymbol{\theta}_{0}(\boldsymbol{\omega}))\overset{\mathcal{D}}{\rightarrow}N_{p}(\mathbf{0},\mathbf{V}_{0})$ with asymptotic variance given by the $p\times p$ matrix

$\displaystyle\mathbf{V}_{0}(\boldsymbol{\omega})=\mathbf{V}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})=\mathbf{H}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})^{-1}\mathbf{K}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})\mathbf{H}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})^{-1},$

(5)

where

$$\mathbf{H}(\boldsymbol{\theta},\boldsymbol{\omega})=\sum_{m=1}^{M}\omega_{m}\mathbf{H}_{m}(\boldsymbol{\theta}),\ \ \mathbf{K}(\boldsymbol{\theta},\boldsymbol{\omega})=Var\left[\sum_{m=1}^{M}\omega_{m}\boldsymbol{U}_{m}(\boldsymbol{\theta})\right],$$

$\mathbf{H}_{m}(\boldsymbol{\theta})=Var(\boldsymbol{U}_{m}(\boldsymbol{\theta}))$ is the $p\times p$ sensitivity matrix for the $m$th component, and $\boldsymbol{U}_{m}(\boldsymbol{\theta})=\nabla_{\boldsymbol{\theta}}\log f_{m}(\boldsymbol{X}\in\mathcal{A}_{m}|\boldsymbol{\theta})$.

The main approach followed here is to minimize, in some sense, the asymptotic variance, $\mathbf{V}_{0}(\boldsymbol{\omega})$. To this end, theory of unbiased estimating equations suggests that both matrices $\mathbf{H}$ and $\mathbf{K}$ should be considered in order to achieve this goal (e.g., see ?). On one hand, note that $\mathbf{H}$ measures the covariance between the composite likelihood score $\boldsymbol{U}(\boldsymbol{\theta},\boldsymbol{\omega})$$=\sum_{m=1}^{M}\omega_{m}\boldsymbol{U}_{m}(\boldsymbol{\theta})=\sum_{m=1}^{M}\omega_{m}\nabla_{\boldsymbol{\theta}}\log f_{m}(\boldsymbol{X}\in\mathcal{A}_{m}|\boldsymbol{\theta})$ and the MLE score $\boldsymbol{U}_{mle}(\boldsymbol{\theta})=\nabla_{\boldsymbol{\theta}}\log f(\boldsymbol{X};\boldsymbol{\theta})$. To see this, differentiate both sides in $E[\boldsymbol{U}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})]=\mathbf{0}$ and obtain $\mathbf{H}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})=E\left[\boldsymbol{U}_{mle}(\boldsymbol{\theta}_{0})\boldsymbol{U}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})^{T}\right].$ This shows that adding sub-likelihood components is desirable, since it increases the covariance with the full likelihood. On the other hand, including too many correlated sub-likelihoods components inflates the variance through the covariance terms in $\mathbf{K}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})$.

The objective of minimizing the asymptotic variance is still undefined since $\mathbf{V}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})$ in (5) is a $p\times p$ positive semidefinite matrix. Therefore, we consider the following one-dimensional objective function

$\displaystyle g_{0}(\boldsymbol{\omega})=\log\det\{\mathbf{V}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})\}=\log\det\{\mathbf{K}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})\}-2\log\det\{\mathbf{H}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})\}.$

(6)

The minimizer, $\boldsymbol{\omega}_{0}$, of the ideal objective (6) corresponds to the $O_{F}$-optimal solution (fixed-sample optimality) (e.g., see ?, Chapter 2). Clearly, such a program still lacks practical relevance, since (6) depends on the unknown parameter $\boldsymbol{\theta}_{0}$. Therefore, $g_{0}(\boldsymbol{\omega})$ should be replaced by some sample-based estimate, say $\hat{g}_{0}(\boldsymbol{\omega})$.

One option is to use the following consistent estimates of $\mathbf{H}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})$ and $\mathbf{K}(\boldsymbol{\theta}_{0},\boldsymbol{\omega})$ in (6):

$\displaystyle\hat{\mathbf{H}}(\boldsymbol{\omega})=\dfrac{1}{n-1}\sum_{m=1}^{M}\sum_{i=1}^{n}\boldsymbol{\omega}_{m}\boldsymbol{U}^{(i)}_{m}(\hat{\boldsymbol{\boldsymbol{\theta}}})\boldsymbol{U}^{(i)}_{m}(\hat{\boldsymbol{\boldsymbol{\theta}}})^{T},\ \ \hat{\mathbf{K}}(\boldsymbol{\omega})=\dfrac{1}{n}\sum_{i=1}^{n}\boldsymbol{U}^{(i)}(\hat{\boldsymbol{\boldsymbol{\theta}}},\boldsymbol{\omega})\boldsymbol{U}^{(i)}(\hat{\boldsymbol{\boldsymbol{\theta}}},\boldsymbol{\omega})^{T},$

(7)

where the estimator $\hat{\boldsymbol{\boldsymbol{\theta}}}=\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})$ is the McLE. Although this strategy works in simple models when $n$ is relatively large and $M$ is small, the estimator $\hat{\mathbf{K}}(\boldsymbol{\omega})$ is knowingly unstable when $n$ is small compared to $\dim(\boldsymbol{\Theta})$ (?). Another issue with this approach in high-dimensional datasets is that the number of operations required to compute $\hat{\mathbf{K}}(\boldsymbol{\omega})$ (or $\mathbf{K}(\boldsymbol{\theta},\boldsymbol{\omega})$) increases quadratically in $\sum_{m=1}^{M}\omega_{m}$.

To reduce the computational burden and avoid numerical instabilities, we estimate $g_{0}$ by the following one-step jackknife criterion:

$$\hat{g}(\boldsymbol{\omega})=\log\det\left\{\sum_{i=1}^{n}\left(\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})^{(-i)}-\overline{\boldsymbol{\theta}}(\boldsymbol{\omega})\right)\left(\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})^{(-i)}-\overline{\boldsymbol{\theta}}(\boldsymbol{\omega})\right)^{T}\right\},$$

(8)

where the pseudo-value $\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})^{(-i)}$ is a composite likelihood estimator based on a sample without observation $\boldsymbol{X}^{(i)}$, and $\overline{\boldsymbol{\theta}}(\boldsymbol{\omega})=\sum_{i=1}^{n}\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})^{(-i)}/n$. Alternatively, one could use the delete-$k$ jackknife estimate where $k>1$ observations at the time are deleted to compute the pseudo-values. The delete-$k$ version is computationally cheaper than delete-$1$ jackknife and therefore should be preferred when the sample size, $n$, is moderate or large. Other approaches – including bootstrap – should be considered depending on the model set up. For example, block re-sampling techniques such as the block-bootstrap (see ? and subsequent papers) are viable options for spatial data and time series.

When the sub-likelihood scores are in closed form, the pseudo-values can be efficiently approximated by the following one-step Newton-Raphson iteration:

$$\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})^{(-i)}=\tilde{\boldsymbol{\theta}}+\left(\sum_{m=1}^{M}\omega_{m}\sum_{j\neq i}\boldsymbol{U}_{m}^{(j)}(\tilde{\boldsymbol{\theta}})\boldsymbol{U}_{m}^{(j)}(\tilde{\boldsymbol{\theta}})^{T}\right)^{-1}\left(\sum_{m=1}^{M}\omega_{m}\sum_{j\neq i}\boldsymbol{U}_{m}^{(j)}(\tilde{\boldsymbol{\theta}})\right),$$

(9)

where $\tilde{\boldsymbol{\theta}}$ is any root-$n$ consistent estimator of $\boldsymbol{\theta}$. Note that $\tilde{\boldsymbol{\theta}}$ does not need to coincide with the McLE, $\hat{\boldsymbol{\boldsymbol{\theta}}}(\boldsymbol{\omega})$, so a computationally cheap initial estimator – based only on a few small sub-likelihoods subset – may be considered. Remarkably, the number of operations required in the Newton-Raphson iteration (9) grows linearly in the number of sub-likelihood components, so that the one-step jackknife objective has computational complexity comparable to a single evaluation of the composite likelihood function (4). Large sample properties of jackknife estimator of McLE’s asympotic variance can be derived under regularity conditions on $\boldsymbol{U}_{m}$ and $\mathbf{H}_{m}$ analogous to those described in ?. Then, for the one-step jackknife using the root-$n$ consistent starting point $\tilde{\boldsymbol{\theta}}$, we have $n\left(\hat{g}(\boldsymbol{\omega})-g_{0}(\boldsymbol{\omega})\right)\overset{p}{\rightarrow}0,\ \ \text{as}\ \ n\rightarrow\infty.$ uniformly on $\boldsymbol{\Omega}$. Moreover, the estimator $\hat{g}(\boldsymbol{\omega})$ is asymptotically equivalent to the classic jackknife estimator.

When computing the McLE of $\boldsymbol{\theta}$, our objective is to find an optimal binary vector $\boldsymbol{\omega}^{\ast}$ to estimate $\boldsymbol{\omega}_{0}={\text{argmin}}_{\boldsymbol{\omega}\in\Omega}\ g_{0}(\boldsymbol{\omega}),$ where $g_{0}(\cdot)$ is the ideal objective function defined in (6). Typically, the population quantity $g_{0}$ cannot be directly assessed, so we replace $g_{0}$ with the sample-based jacknife estimate $\hat{g}$ described in Section 2 and aim at finding

$$\boldsymbol{\omega}^{\ast}=\underset{\boldsymbol{\omega}\in\Omega}{\text{argmin}}\ \hat{g}(\boldsymbol{\omega}).$$

This task, however, is computationally infeasible through enumerating the space $\boldsymbol{\Omega}$ if $d$ is even moderately large. For example, for composite likelihoods defined based on all pairs of variables $(X_{s},X_{t})$, $1\leq s<t\leq d$ (pairwise likelihood), $\boldsymbol{\Omega}$ contains $2^{{d}\choose{2}}=2^{190}$ elements when $d=20$.

To overcome this enumeration complexity, we carry out a random search method based on Gibbs sampling. We regard the weight vector $\boldsymbol{\omega}$ as a random vector following the joint probability mass function (pmf)

$$\pi_{\tau}(\boldsymbol{\omega})=\dfrac{1}{Z(\tau)}\exp\left\{-\tau\hat{g}(\boldsymbol{\omega})\right\},\ \ \boldsymbol{\omega}\in\boldsymbol{\Omega},$$

(10)

where $Z(\tau)=\sum_{\boldsymbol{\omega}\in\boldsymbol{\Omega}}\exp\{-\tau\hat{g}(\boldsymbol{\omega})\}$ is the normalizing constant. The above distribution depends on the tuning parameter $\tau>0$, which controls the extent to which we emphasize larger probabilities (and reduce smaller probabilities) on $\boldsymbol{\Omega}$. Then $\boldsymbol{\omega}^{*}$ is also the mode of $\pi_{\tau}(\boldsymbol{\omega})$, meaning that $\boldsymbol{\omega}^{*}$ will have the highest probability to appear, and will be more likely to appear earlier rather than later, if a random sequence of $\boldsymbol{\omega}$ is to be generated from $\pi_{\tau}(\boldsymbol{\omega})$.

Therefore, estimating $\boldsymbol{\omega}^{*}$ can be readily done based on the random sequence generated from $\pi_{\tau}(\boldsymbol{\omega})$. But generating a random sample from $\pi_{\tau}(\boldsymbol{\omega})$ directly is difficult because $\pi_{\tau}(\boldsymbol{\omega})$ contains an intractable normalizing constant $Z(\tau)$. Instead, we will generate a Markov chain using the product of all univariate conditional pmf’s with respect to $\pi_{\tau}(\boldsymbol{\omega})$ as the transitional kernel. The stationary distribution of such a Markov chain can be proved to be $\pi_{\tau}(\boldsymbol{\omega})$ (?, Chapter 10). Hence, the part of this Markov chain after reaching equilibrium can be regarded as a random sample from $\pi_{\tau}(\boldsymbol{\omega})$ for most purposes. The MCMC method just described is in fact the so-called Gibbs sampling method. The key factor for the Gibbs sampling to work is all the univariate conditional probability distributions of the target distribution can be relatively easily simulated.

Let us write $\boldsymbol{\omega}=(\omega_{1},\cdots,\omega_{M})$; $\boldsymbol{\omega}_{m_{1}:m_{2}}=(\omega_{m_{1}},\omega_{m_{1}+1},\cdots,\omega_{m_{2}})$, if $m_{1}\leq m_{2}$, and $\boldsymbol{\omega}_{m_{1}:m_{2}}=\emptyset$ otherwise; and $\boldsymbol{\omega}_{-m}=(\omega_{1},\cdots,\omega_{m-1},\omega_{m+1},\cdots,\omega_{M})$. Then it is easy to see that the conditional pmf of $\omega_{m}$ given $\boldsymbol{\omega}_{-m}$, $m=1,\cdots,M$, is

$$\pi_{\tau}(\omega_{m}|\boldsymbol{\omega}_{-m})=\frac{\exp\{-\tau\hat{g}(\boldsymbol{\omega})\}}{\exp\{-\tau\hat{g}(\boldsymbol{\omega}_{m}^{[0]})\}+\exp\{-\tau\hat{g}(\boldsymbol{\omega}_{m}^{[1]})\}},\quad\omega_{m}=0,1,$$

(11)

where $\omega_{m}^{[j]}=(\omega_{1:(m-1)},j,\omega_{(m+1):M})$, $j=0,1$. Note that $\pi_{\tau}(\omega_{m}|\boldsymbol{\omega}_{-m})$ is simply a Bernoulli pmf and so it is easily generated. For each binary vector $\boldsymbol{\omega}$, $\hat{g}(\boldsymbol{\omega})$ is computed using the one-step jackknife estimator described in Section 2. Therefore, the probability mass function $\pi_{\tau}(\boldsymbol{\omega})$ is well defined for any $\tau>0$. On the other hand, we have shown that the Gibbs sampling can be used to generate a Markov chain from $\pi_{\tau}(\boldsymbol{\omega})$, from which we can find a consistent estimator $\hat{\boldsymbol{\omega}}^{*}$ for the mode $\boldsymbol{\omega}^{*}$ if $\boldsymbol{\omega}^{*}$ is unique. Consequently, the universally maximum composite log-likelihood estimator of $\boldsymbol{\theta}$ can be approximated by the McLE, $\hat{\boldsymbol{\boldsymbol{\theta}}}({\hat{\boldsymbol{\omega}}}^{*})$.

Note that the mode $\boldsymbol{\omega}^{*}$ is not necessarily unique. In this case one can still consider consistent estimation for $\boldsymbol{\omega}^{*}$ but in the following meaning. We know $\hat{g}(\boldsymbol{\omega}^{*})$ is the minimum and always unique according to its definition. Thus $\hat{g}(\boldsymbol{\omega}^{*})$ can be consistently and uniquely estimated based on the Markov chain of $\hat{g}(\boldsymbol{\omega})$ induced from the Markov chain of $\boldsymbol{\omega}$ generated from $\pi_{\tau}(\boldsymbol{\omega})$ when the length of the Markov chain goes to infinity. Therefore, any estimator $\hat{\boldsymbol{\omega}}^{*}$ such that $\hat{g}(\hat{\boldsymbol{\omega}}^{*})$ becomes consistent with $\hat{g}(\boldsymbol{\omega}^{*})$ can be regarded as a consistent estimator of $\boldsymbol{\omega}^{*}$. Consequently, the McLE $\hat{\boldsymbol{\boldsymbol{\theta}}}({\hat{\boldsymbol{\omega}}}^{*})$ for each such $\hat{\boldsymbol{\omega}}^{*}$ is still the universally McLE but not necessarily unique. From a practitioner’s view point there is no need to identify all consistent estimators of $\boldsymbol{\omega}^{*}$ and all universally McLE of $\boldsymbol{\theta}$. Finding out one such $\hat{\boldsymbol{\omega}}^{*}$ or a tight superset of it would be a sufficient advance in parsimonious and efficient likelihood composition.

The above discussion motivates the steps of our core Gibbs sampling algorithm for simultaneous composite likelihood estimation and selection. Let $\tau$ be given and fixed.

1. 0.

   For $t=0$, choose an initial binary vector of weights $\boldsymbol{\omega}^{(0)}$ and compute the one-step jackknife estimator $\hat{g}(\boldsymbol{\omega}^{(0)})$. E.g. randomly set 5 elements of $\boldsymbol{\omega}^{(0)}$ to 1 and the rest to 0.

2. 1.

   For each $t=1,\cdots,T$ for a given $T$, obtain $\boldsymbol{\omega}^{(t)}$ by repeating $1.1$ to $1.4$ for each $m=1,\cdots,M$ sequentially.

   1. 1.1.

      Compute, if not available yet, $\hat{g}(\omega_{1:(m-1)}^{(t)},j,\omega_{(m+1):M}^{(t-1)})$, $j=0,1$.

   2. 1.2.

      Compute the conditional pmf of $\omega_{m}$ given $(\omega_{1:(m-1)}^{(t)},\omega_{(m+1):M}^{(t-1)})$:

      $$\pi_{\tau}(\omega_{m}=j|\omega_{1:(m-1)}^{(t)},\omega_{(m+1):M}^{(t-1)})\propto\exp\{-\tau\hat{g}(\omega_{1:(m-1)}^{(t)},j,\omega_{(m+1):M}^{(t-1)})\}$$

      (12)

      where $j=0,1$. Note this is a Bernoulli pmf.

   3. 1.3.

      Generate a random number from the Bernoulli pmf obtained in $1.2$, and denote the result as $\omega_{m}^{(t)}$.

   4. 1.4.

      Set $\boldsymbol{\omega}^{(t)}\leftarrow(\omega_{1:(m-1)}^{(t)},\omega_{m}^{(t)},\omega_{(m+1):M}^{(t-1)})^{T}$. Also compute and record $\hat{g}(\boldsymbol{\omega}^{(t)})$.

3. 2.

   Compute $\hat{\boldsymbol{\omega}}^{*}=\arg\min_{1\leq t\leq T}\hat{g}(\boldsymbol{\omega}^{(t)}),$ and regard it as the estimate of $\boldsymbol{\omega}^{*}$. Alternatively, column-combine $\boldsymbol{\omega}^{(1)},\cdots,\boldsymbol{\omega}^{(T)}$ generated in Step $1$ into an $M\times T$ matrix $\hat{\mathbf{W}}$; then compute row averages of $\hat{\mathbf{W}}$, say $\overline{\omega}_{1},\cdots,\overline{\omega}_{M}$, and set $\tilde{\boldsymbol{\omega}}^{*}=(\tilde{\omega}_{1}^{*},\cdots,\tilde{\omega}_{M}^{*})$ where $\tilde{\omega}_{m}^{*}=1$, if $\overline{\omega}_{m}\geq\xi$, and $\tilde{\omega}_{m}^{*}=0$ if $\overline{\omega}_{m}<\xi$, where $\xi$ is some constant larger than $0.5$.

4. 3.

   Finally, compute $\hat{\boldsymbol{\boldsymbol{\theta}}}(\hat{\boldsymbol{\omega}}^{*})$ (or $\hat{\boldsymbol{\boldsymbol{\theta}}}(\tilde{{\boldsymbol{\omega}}}^{*})$) and $\hat{g}(\hat{\boldsymbol{\omega}}^{*})$ (or $\hat{g}(\tilde{{\boldsymbol{\omega}}}^{*})$).

Firstly, note that Gibbs sampling has been used in various contexts in the literature of model selection. ? used a similar strategy to generate the distribution of the variable indicators in Bayesian linear regression. ? used Gibbs sampler for selecting robust linear regression models. ? used the Gibbs sampler for selecting generalized linear regression models. ? and ? used Gibbs sampling for selection in the context of time series models. However, to our knowledge this is the first work proposing a general-purpose Gibbs sampler for construction of composite likelihoods.

Secondly, the sequence $\boldsymbol{\omega}^{(1)},\dots,\boldsymbol{\omega}^{(T)}$ is a Markov chain, which requires an initial vector $\boldsymbol{\omega}^{(0)}$ and a burn-in period to be in equilibrium. Values of $\boldsymbol{\omega}^{(0)}$ do not affect the eventual attainment of equilibrium so can be arbitrarily chosen. From a computational point of view most components of $\boldsymbol{\omega}^{(0)}$ should be set to 0 to reduce computing load. For example, we can randomly set all but 5 of the components to 0. To assess whether the chain has reached equilibrium, we suggest the control-chart method discussed in ?. For the random variable $\hat{g}({\boldsymbol{\omega}})$, we have the following probability inequality for any given $b>1$:

$$Pr\left(\hat{g}(\boldsymbol{\omega})-\min\hat{g}(\boldsymbol{\omega})\geq b\sqrt{Var[\hat{g}(\boldsymbol{\omega})]+(E[\hat{g}(\boldsymbol{\omega})]-\min\hat{g}(\boldsymbol{\omega}))^{2}}\right)\leq\dfrac{1}{b^{2}}$$

This inequality can be used to find an upper control limit for $\hat{g}({\boldsymbol{\omega}})$. For example, by setting $b=\sqrt{10}$, an at least $90\%$ upper control limit for $\hat{g}({\boldsymbol{\omega}})$ can be estimated as $\hat{g}^{\ast}+\sqrt{10s^{2}+10(\overline{g}-\hat{g}^{\ast})^{2}},$ where $\hat{g}^{\ast}$, $\overline{g}$ and $s^{2}$ are the minimum, sample mean and sample variance based on the first $N$ observations, $\hat{g}(\boldsymbol{\omega}^{(1)}),\dots,\hat{g}(\boldsymbol{\omega}^{(N)})$, $N<T$, where typically we set $N=\lfloor T/2\rfloor$. We then count the number of observations passing the upper control limit in the remaining sample $\hat{g}(\boldsymbol{\omega}^{(N+1)}),\dots,\hat{g}(\boldsymbol{\omega}^{(T)})$. If more than 10% of the points are above the control limit, then at a significance level not more than 90% there is statistical evidence against equilibrium. Upper control limits of different levels for $\hat{g}({\boldsymbol{\omega}})$ can be similarly calculated and interpreted by choosing different values of $b$.

Thirdly, $\hat{\boldsymbol{\omega}}^{*}$ computed in Step 2 is simply a sample mode of $\pi_{\tau}(\boldsymbol{\omega})$ based on its definition in (10). Hence, by the Ergodic Theorem for stationary Markov chain, $\hat{\boldsymbol{\omega}}^{*}$ is a strongly consistent estimator of $\boldsymbol{\omega}$ under minimal regularity conditions. With similar arguments $\overline{\omega}_{m}$ is a strongly consistent estimator of the success probability involved in the marginal distribution of $\omega_{m}$ induced from $\pi_{\tau}(\boldsymbol{\omega})$. Hence, it is not difficult to see that the resultant estimator $\tilde{\boldsymbol{\omega}}^{*}$ should satisfy $\tilde{\omega}_{m}^{*}\geq\omega_{m}^{*}$ without requiring $T$ to be very large. Propositions 1 and 2 in ? provide an exposition of this property. Therefore, the estimator $\tilde{\boldsymbol{\omega}}^{*}$ captures all informative sub-likelihood components with high probability.

Finally, the tuning constant $\tau$ adjusts the mixing behavior of the chain, which has important consequences on the exploration/exploitation trade-off on the search space $\boldsymbol{\Omega}$. If $\tau$ is too small, the algorithm produces solutions approaching the global optimal value $\boldsymbol{\omega}^{\ast}$ slowly; if $\tau$ is large, then the algorithm finds local optima and may not reach $\boldsymbol{\omega}^{\ast}$. The former behavior corresponds to a rapidly mixing chain, while the latter occurs when the chain is mixing too slowly. In the composite likelihood selection setting, the main hurdle is the computational cost, so $\tau$ should be set according to the available computing capacity, after running some graphical or numerical diagnostics (e.g., see ?). We choose to use $\tau=d$ in our empirical study, which does not seem to create adverse effects.

Without additional modifications, Algorithm 1 ignores the likelihood complexity, since solutions with many sub-likelihoods have in principle the same chance to occur as those with fewer components. To discourage selection of overly complex likelihoods, we augment the Gibbs distribution (10) as follows:

$$\pi_{\tau,\lambda}(\boldsymbol{\omega})=Z(\tau,\lambda)^{-1}\exp\{-\tau\hat{g}_{\lambda}(\boldsymbol{\omega})\},$$

(13)

where

$$\hat{g}_{\lambda}(\boldsymbol{\omega})=\hat{g}(\boldsymbol{\omega})+{pen}(\boldsymbol{\omega}|\lambda),\quad\tau>0,\;\lambda>0,$$

(14)

$\hat{g}(\boldsymbol{\omega})$ is the jackknifed variance objective defined in Section 2, $Z(\tau,\lambda)$ is the normalization constant, and $pen(\boldsymbol{\omega})$ is a complexity penalty enforcing sparse solutions when $\dim(\boldsymbol{\Omega})$ is large. Maximization of $\pi_{\tau,\lambda}(\boldsymbol{\omega})$ is interpreted as a maximum a posteriori estimation for $\boldsymbol{\omega}$, where the probability distribution proportional to $\exp\{-pen(\boldsymbol{\omega}|\lambda)\}$ is regarded as a prior pmf over $\boldsymbol{\Omega}$. In this paper, we use the penalty term of form $\text{pen}(\boldsymbol{\omega}|\lambda)=\lambda\sum_{m=1}^{M}\omega_{m}$, since it corresponds to well-established model-selection criteria. For example, choices $\lambda=1$, $\lambda=2^{-1}\log n$ and $\lambda=\log\log n$ correspond to the AIC, BIC and HQC criteria, respectively (e.g., see ?). Other penalties could be considered as well depending on the model structure and available prior information; however, these are not shown to be crucial based on our empirical study so will not be explored in this paper.

To find the optimal solution $\boldsymbol{\omega}^{\ast}$, one could compute a sequence of optimal values $\hat{\boldsymbol{\omega}}^{\ast}_{\lambda_{1}},\dots,\hat{\boldsymbol{\omega}}^{\ast}_{\lambda_{B}}$ and then take $\min_{1\leq b\leq B}\hat{g}(\hat{\boldsymbol{\omega}}^{\ast}_{\lambda_{b}})$. There are, however, various issues in this approach: first, the globally optimal value $\boldsymbol{\omega}^{\ast}$ might not be a member of the set $\{\hat{\boldsymbol{\omega}}^{\ast}_{\lambda_{b}}\}_{b=1}^{B}$, since the mode of $\pi_{\tau,\lambda}(\boldsymbol{\omega})$ is not necessarily the composite likelihood solution which minimizes $\hat{g}(\boldsymbol{\omega})$. Second, even if $\boldsymbol{\omega}^{\ast}$ is in such a set, determining $\lambda$ is typically challenging. To address the above issues, we employ the idea of stability selection, introduced by ? in the context of variable selection for linear models. Given an arbitrary value for $\lambda$, stochastic stability exploits the variability of random samples generated from $\pi_{\tau,\lambda}(\boldsymbol{\omega})$ by the Gibbs procedure, say $\boldsymbol{\omega}^{(1)}_{\lambda},\dots,\boldsymbol{\omega}^{(T)}_{\lambda}$ and choses all the partial likelihoods that occur in a large fraction of generated samples. For a given $0<\xi<1$, we define the set of stable likelihoods by the vector $\hat{\boldsymbol{\omega}}^{\text{stable}}$, with elements

$$\hat{\omega}_{m}^{\text{stable}}=\left\{\begin{array}[]{ll}1,&\text{if }\ \ \dfrac{1}{T}\sum_{t=1}^{T}\omega_{\lambda,m}^{(t)}\geq\xi,\\ 0,&\text{otherwise},\end{array}\right.$$

(15)

so we regard as stable those sub-likelihoods selected more frequently and disregard sub-likelihood items with low selection probabilities. Following ?, we choose the tuning constant $\xi$ using the following bound on the expected number of false selections, $V$:

$$E(V)\leq\dfrac{1}{(2\xi-1)}\dfrac{\eta_{\lambda}}{M},$$

(16)

where $\eta_{\lambda}$ is average number of selected sub-likelihood components. In multiple testing, the quantity $\alpha=E(V)/M$ is sometimes called the per-comparison error rate (PCER). By increasing $\xi$, only few likelihood components are selected, so that we reduce the expected number of falsely selected variables. We choose the threshold $\xi$ by fixing the PCER at some desired value (e.g., $\alpha=0.10$), and then choose $\xi$ corresponding to the desired error rate. The unknown quantity $\eta_{\lambda}$ in our setting can be estimated by the average number of sub-likelihood components over $T$ Gibbs samples.

Finally, note that tuning $\xi$ according to (16) makes redundant the determination of the optimal $\lambda$ value as long as $pen(\boldsymbol{\omega}|\lambda)$ in (14) is not dominant over $\hat{g}({\boldsymbol{\omega}})$. This is further supported by our empirical study where we found the effect of $\lambda$ on $\hat{\boldsymbol{\omega}}^{\text{stable}}$ is negligible.

The preceding discussions lead to Algorithm 2 which is essentially the same as Algorithm 1 with two exceptions: (i) we replace $\hat{g}$ in Algorithm 1 by the augmented objective function $\hat{g}_{\lambda}$ defined in (14); (ii) Steps 2 in Algorithm 1 is replaced by the following stability selection step.

• $2^{\prime}$.

  Column-combine $\boldsymbol{\omega}^{(1)},\cdots,\boldsymbol{\omega}^{(T)}$ generated in Step $1$ into an $M\times T$ matrix $\hat{\mathbf{W}}$. Compute the row averages of $\hat{\mathbf{W}}$, denoted as $(\hat{\omega}_{1},\cdots,\hat{\omega}_{M})$, i.e. $\hat{\omega}_{m}=T^{-1}\sum_{t=1}^{T}\omega_{m}^{(t)}$, $m=1,\cdots,M$. Then set $\hat{\boldsymbol{\omega}}^{\text{stable}}=(\hat{\omega}_{1}^{*},\cdots,\hat{\omega}_{M}^{*})$ where $\hat{\omega}_{m}^{*}=1$ if $\hat{\omega}_{m}\geq\hat{\xi}_{\lambda}$ and $\hat{\omega}_{m}^{*}=0$ if $\hat{\omega}_{m}<\hat{\xi}_{\lambda}$, $m=1,\cdots,M$, where $\hat{\xi}_{\lambda}=\dfrac{1}{2}\left(\dfrac{\hat{\eta}_{\lambda}}{\alpha M^{2}}+1\right),$ and $\alpha$ is a nominal level for per-comparison error rate (e.g., 0.05 or 0.1).

The estimated threshold $\hat{\xi}_{\lambda}$ is obtained from (16) by plugging-in $\hat{\eta}_{\lambda}=\sum_{t=1}^{T}\sum_{m=1}^{M}\omega^{(t)}_{m}/T$, the sample average of the numbers of selected sub-likelihood components in $T$ Gibbs samples.

Let $\boldsymbol{X}\sim N_{d}(\mu\bf 1,\mathbf{\Sigma})$, where the parameter of interest is the common location $\mu$. We study the scenario where many components bring redundant information on $\mu$ by considering covariance matrix $\mathbf{\Sigma}$ with elements $\{\mathbf{\Sigma}\}_{mm}=1$, for all $1\leq m\leq d$, and off-diagonal elements $\{\mathbf{\Sigma}\}_{lm}=\rho>0$ if $l,m\leq d^{\ast}$, for some $d^{\ast}<d$, while $\{\mathbf{\Sigma}\}_{lm}=0$ elsewhere.

We consider one-wise score composite likelihood estimator solving $0=\sum_{m=1}^{d}\omega_{m}U_{m}(\mu)$, where $U_{m}(\mu)=\sum_{i=1}^{n}(X^{(i)}_{m}-\mu)$. It is easy to find the composite likelihood estimator is $\widehat{\mu}_{cl}(\boldsymbol{\omega})={\sum_{m=1}^{d}\omega_{m}\overline{X}_{m}}/{\sum_{m=1}^{d}\omega_{m}}$ where $\overline{X}_{m}=\sum_{i=1}^{n}X_{m}^{(i)}/n$. For this simple model, the jackknife criterion $\hat{g}(\cdot)$ can be easily computed in closed form. The pseudo-values are $\widehat{\mu}^{(-i)}_{cl}(\boldsymbol{\omega})=\sum_{m=1}^{d}\omega_{m}\overline{X}^{(-i)}_{m}/\sum_{m=1}^{d}\omega_{m}$, and the average of pseudo-values is $\widehat{\mu}_{cl}$. It can be shown that

$$\hat{g}(\boldsymbol{\omega})=\log\sum_{i=1}^{n}\left(\sum_{m=1}^{d}\omega_{m}(X_{m}^{(i)}-\overline{X}_{m})\right)^{2}-2\log\left(\sum_{m=1}^{d}\omega_{m}\right),$$

up to a constant not depending on $\boldsymbol{\omega}$. It can also be shown that $O_{F}$-criterion has the following expression

$\displaystyle g_{0}(\boldsymbol{\omega})=\log Var\left(\widehat{\mu}_{cl}(\boldsymbol{\omega})\right)\propto\log\left(\sum_{m=1}^{d}\omega_{m}+2\rho\sum_{l<m\leq d^{\ast}}\omega_{l}\omega_{m}\right)-2\log\left(\sum_{m=1}^{d}\omega_{m}\right),$

(17)

depending on the unknown parameter $\rho$. This suggests that including too many correlated (redundant) components (with $\rho\neq 0$) damages McLE’s variance as $d$ increases. Particularly, setting $\omega_{j}=1$, for all $j=1,\dots,d$, implies $Var(\widehat{\mu}_{cl}(\boldsymbol{\omega}))=O(1)$, while choosing only uncorrelated sub-likelihoods ($\omega_{j}=0$ if $2\leq j\leq d^{\ast}$, and $\omega_{j}=1$ elsewhere), gives $Var(\widehat{\mu}_{cl}(\boldsymbol{\omega}))=O(d^{-1})$.

In Table 1, we show Monte Carlo simulation results from $B=250$ runs of Algorithm 1 (MCMC-CLS1) using the two approaches described in Section 3.2: one consists of choosing the best weights vector (CLS1 min); the other uses thresholding, i.e. selects elements when the weights are selected with a sufficiently large frequency (CLS1 thres.). In the same table we also report results on the Algorithm 2 (MCMC-CLS2) based on the stochastic stability selection approach. Our algorithms are compared with the estimator including all one-wise sub-likelihood components (No selection) and the optimal maximum likelihood estimator (MLE). We compute the MLE of $\mu$ in two ways: either based on using the known $\Sigma$ value or based on using the sample covariance $\hat{\Sigma}=(n-1)^{-1}\sum_{i=1}^{n}(X_{i}-\overline{X})(X_{i}-\overline{X})^{T}$. Note that the latter is not available when $d>n$. Note our algorithms are designed to obtain simultaneous estimation and dimension reduction in the presence of limited information (i.e. $d>n$ or $d\gg n$) where the optimal weights are difficult to estimate from the data. Stochastic selection with $\xi=0.7$ and $T=10d$ was carried out for samples of $n=5,25$ and $100$ observations from a model with $d^{\ast}=0.8\times d$ correlated components, with $d=10,30$. To compare the methods we computed Monte Carlo estimates of the variance ($Var$) and squared bias ($Bias^{2}$) of $\widehat{\mu}_{cl}(\boldsymbol{\omega})$. Table 2 further reports the average number of selected likelihood components (no. comp).

For all considered data dimensions, our selection methods outperform the all one-wise likelihood (AOW or No selection) estimator and show relatively small losses in terms of mean squared error ($Var+Bias^{2}$) compared to the MLEs. The gains in terms of variance reduction are particularly evident for larger data dimensions (e.g., see $d=30$). At the same time, our method tends to select mostly uncorrelated sub-likelihoods, while discarding the redundant components which do not contribute useful information on $\mu$.

Next, we illustrate the selection procedure based on MCMC-CLS2 (Algorithm 2). We draw a random sample of $n=50$ observations from the model $N_{d}(\mu\bf 1,\mathbf{\Sigma}(\rho))$ described above with $d=250$, $d^{\ast}=0.8d$, and $\rho=0.9$. This corresponds to 200 strongly redundant variables and 50 independent variables. We applied Algorithm 2 with $\alpha=0.1$ and $\lambda=1$ (corresponding to the AIC penalty). In Figure 1 (b), we show the relative frequencies of the components, $\overline{\omega}_{m}=\sum_{t=1}^{T}\omega^{(t)}_{m}$, $m=1,\dots,250$, in $T=1000$ MCMC iterations. As expected, the uncorrelated sub-likelihood components (components 201–250) are sampled much more frequently than the redundant ones (components 1–200). Figure 1 (c) shows the objective function $\hat{g}_{\lambda}(\hat{\boldsymbol{\omega}}^{\text{stable}})$ (up to a constant) evaluated at the best solution computed from past MCMC samples ($\xi=0.7$). Figure 1 (d) shows the Hamming distance, $dist(\boldsymbol{\omega},\boldsymbol{\omega}^{\prime})=\sum_{j}I(\omega_{j}\neq\omega_{j}^{\prime})$, between the current selected rule, $\hat{\boldsymbol{\omega}}^{\ast}$ and true optimal value, $\boldsymbol{\omega}^{\ast}=(1,\underbrace{0,\dots,0}_{199},\underbrace{1,\dots,1}_{50})$. Overall, our algorithm quickly approaches the optimal solution and the final selection has 94.0% asymptotic relative efficiency (ARE) compared to MLE. When no stochastic selection is applied and all 250 sub-likelihood components are included, the relative efficiency is only 13%. As far as computing time is concerned, our non-optimized R implementations of the MCMC-CLS1 and MCMC-CLS2 algorithms for the above example takes, respectively 4.25 and 1.87 seconds per MCMC iteration, respectively. The computing time was measured on a laptop computer with Intel ®Core${}^{\text{TM}}$ i7-2620M CPU 2.70GHz.

Let $\boldsymbol{X}\sim N_{d}(\bf 0,\mathbf{\Sigma}(\rho))$, where $\mathbf{\Sigma}(\rho)=\{(1-\rho)\boldsymbol{I}_{d}+\rho{\bf 1}_{d}{\bf 1}_{d}^{T}\}$ and $0\leq\rho\leq 1$ is the unknown parameter of interest. The marginal univariate sub-likelihoods do not contain information on $\rho$, so we consider pairwise sub-likelihoods

$\displaystyle\ell_{lm}(\rho)=-\frac{n}{2}\log(1-\rho^{2})-\frac{(SS_{ll}+SS_{mm})}{2(1-\rho^{2})}+\frac{\rho SS_{lm}}{1-\rho^{2}},\ \ \ 1\leq l<m\leq d,$

(18)

where $SS_{mm}=\sum_{i=1}^{n}(X^{(i)}_{m})^{2}$ and $SS_{lm}=\sum_{i=1}^{n}X^{(i)}_{l}X^{(i)}_{m}$. Given $\boldsymbol{\omega}$, we estimate $\rho$ by solving the composite score equation $0=\sum_{j<k}\omega_{jk}U_{jk}(\rho)$ by Newton-Raphson iterations where each pairwise score

$$U_{jk}(\rho)=(1+\rho^{2})SS_{jk}-\rho(SS_{jj}+SS_{kk})+n\rho(1-\rho^{2})$$

(19)

is a cubic function of $\rho$. It is well known that the composite likelihood estimation can lead to poor results for this model. ? carries out asymptotic variance calculations, showing that efficiency losses compared to MLE occur whenever $d>2$, with more pronounced efficiency losses for $\rho$ near $0.5$. Particularly, if $d=2$ (exactly one pair), $ARE=1$; if $d>2$, ARE=1 if $\rho$ approaches $0$ or $1$. The next simulation results show that in finite samples composite likelihood selection is advantageous even in such a challenging situation.

Since closed-form pairwise score expressions are available for this model, we use the objective function $\hat{g}(\boldsymbol{\omega})$ based on the one-step jackknife with pseudo-values computed as in Equation (9). In Table 3, we present results from $B=250$ Monte Carlo runs of the MCMC-CLS algorithm applied to the model with $d=5,8,10$ dimensions, which correspond to $M={{d}\choose{2}}=10,28,45$ sub-likelihoods, respectively. We compute Monte Carlo estimates for the finite-sample relative efficiency $RE=\widehat{Var}(\hat{\rho}_{apw})/\widehat{Var}(\hat{\rho}_{cl}(\hat{\boldsymbol{\omega}}^{\ast}))$, where $\widehat{Var}$ denotes the Monte Carlo estimate of the finite-sample variance $\hat{\rho}_{cl}(\hat{\boldsymbol{\omega}}^{\ast})$ is the estimator selected by the MCMC-CLS algorithm, while $\hat{\rho}_{apw}$ is the all pairwise (APW) estimator obtained by including all available pairs (thus $RE>1$ indicates that our stochastic selection outperforms no selection). We show values of $\rho$ around $0.5$, since they correspond to the largest asymptotic efficiency losses of pairwise likelihood compared to MLE (see ?, Figure 1). In all considered cases, our stochastic selection method improves the efficiency of the estimator based on all pairwise components; at the same time, our composite likelihoods employ a much smaller number of components. For example, when $\rho=0.6$ the efficiency improvements range from 8% to 39%, using only about half of the available components.