# Block Gibbs samplers for logistic mixed models: convergence properties and a
comparison with full Gibbs samplers

Yalin Rao and Vivekananda Roy
Department of Statistics, Iowa State University, USA

###### Abstract

The logistic linear mixed model (LLMM) is one of the most widely used
statistical models. Generally, Markov chain Monte Carlo algorithms are used to
explore the posterior densities associated with the Bayesian LLMMs. Polson et
al.’s (2013) Pólya-Gamma data augmentation (DA) technique can be used to
construct full Gibbs (FG) samplers for the LLMMs. Here, we develop efficient
block Gibbs (BG) samplers for Bayesian LLMMs using the Pólya-Gamma DA method.
We compare the FG and BG samplers in the context of a real data example, as
the correlation between the fixed effects and the random effects changes as
well as when the dimensions of the design matrices vary. These numerical
examples demonstrate superior performance of the BG samplers over the FG
samplers. We also derive conditions guaranteeing geometric ergodicity of the
BG Markov chain when the popular improper uniform prior is assigned on the
regression coefficients, and proper or improper priors are placed on the
variance parameters of the random effects. This theoretical result has
important practical implications as it justifies the use of asymptotically
valid Monte Carlo standard errors for Markov chain based estimates of the
posterior quantities.

Key words: Data augmentation; drift condition; geometric ergodicity; GLMM;
Markov chain CLT; MCMC; standard errors

##  1 Introduction

The logistic linear mixed model (LLMM) is an extensively used generalized
linear mixed model for binary response data. Let $(Y_{1},Y_{2},...,Y_{n})$
denote the vector of Bernoulli responses. Let $X$ and $Z$ be the $n\times p$
and $n\times q$ known design matrices corresponding to the fixed effects and
the random effects, respectively. Suppose $x_{i}^{\top}$ and $z_{i}^{\top}$
indicate the $i^{th}$ row of $X$ and $Z$, respectively, for $i=1,...,n$. Let
$\beta\in\mathbb{R}^{p}$ be the regression coefficients vector and
$u\in\mathbb{R}^{q}$ be the random effects vector. In general, a generalized
linear mixed model (GLMM) can be built with a link function that connects the
probability that the response variable $Y$ equals 1 (that is, the expectation
of $Y$) with $X$ and $Z$. For the LLMM,
$P(Y_{i}=1)=F(x_{i}^{\top}\beta+z_{i}^{\top}u)$, where $F$ indicates the
cumulative distribution function for the standard logistic random variable,
that is, $F(t)=e^{t}/(1+e^{t}),\,t\in\mathbb{R}$. Also, assume that we can
divide $u$ into $r$ independent random effects. Let
$u=(u_{1}^{\top},\dots,u_{r}^{\top})^{\top}$, where $u_{j}$ is a $q_{j}\times
1$ vector with $q_{j}>0$, $j=1,\dots,r,$ and $\sum_{j=1}^{r}q_{j}=q$. Assume,
$u_{j}\overset{ind}{\sim}N(0,(1/\tau_{j}){\mathrm{I}}_{q_{j}})$, where
$\tau_{j}>0,$ for $j=1,\dots,r$. Let $\tau=(\tau_{1},...,\tau_{r})$. Thus the
data model for the LLMM is

$\displaystyle Y_{i}\mid\beta,u,\tau$
$\displaystyle\overset{ind}{\sim}Ber(F(x_{i}^{\top}\beta+z_{i}^{\top}u))\quad\text{for}\quad
i=1,...,n,$ $\displaystyle u_{j}\mid\tau_{j}$
$\displaystyle\overset{ind}{\sim}N(0,(1/\tau_{j}){\mathrm{I}}_{q_{j}}),\quad
j=1,...,r.$ (1)

Let $y=(y_{1},y_{2},...,y_{n})$ denote the vector of observed Bernoulli
responses. Then, the likelihood function for $(\beta,\tau)$ is

$\displaystyle L(\beta,\tau\mid
y)=\int_{\mathbb{R}^{q}}\prod_{i=1}^{n}\big{[}F(x_{i}^{\top}\beta+z_{i}^{\top}u)\big{]}^{y_{i}}\big{[}1-F(x_{i}^{\top}\beta+z_{i}^{\top}u)\big{]}^{1-y_{i}}\phi_{q}(u;0,D(\tau)^{-1})du,$
(2)

where $D(\tau)^{-1}=\oplus_{j=1}^{r}\frac{1}{\tau_{j}}{\mathrm{I}}_{q_{j}}$,
and $\oplus$ indicates the direct sum. Here, $\phi_{q}(u;0,D(\tau)^{-1})$
denotes the probability density function of the $q$-dimensional normal
distribution with mean vector $0$, covariance matrix $D(\tau)^{-1}$, and
evaluated at $u$.

When working in a Bayesian framework, one must specify priors for $\beta$ and
$\tau$. Here, we consider the prior for $\beta$ as given by

$\displaystyle\pi(\beta)\propto\exp\big{[}-\frac{1}{2}(\beta-\mu_{0})^{\top}Q(\beta-\mu_{0})\big{]},$
(3)

where $\mu_{0}\in\mathbb{R}^{p}$ and $Q$ is a $p\times p$ positive definite
matrix (proper normal prior) or a zero matrix (improper uniform prior). Thus,
if $Q=0$, then $\pi(\beta)\propto 1$. The prior for $\tau_{j}$ is

$\displaystyle\pi(\tau_{j})\propto\tau_{j}^{a_{j}-1}e^{-\tau_{j}b_{j}},\;j=1,...,r,$
(4)

which may be proper or improper depending on the values of $a_{j}$ and
$b_{j}$. Finally, we assume that $\beta$ and $\tau$ are apriori independent,
and all the $\tau_{j}$s are also apriori independent. Hence, the joint
posterior density for $(\beta,\tau)$ is

$\displaystyle\pi(\beta,\tau\mid y)=\frac{1}{c(y)}L(\beta,\tau\mid
y)\pi(\beta)\pi(\tau),$ (5)

where $c(y)=\int_{\mathbb{R}_{+}^{r}}\int_{\mathbb{R}^{p}}L(\beta,\tau\mid
y)\pi(\beta)\pi(\tau)d\beta d\tau$ is the marginal pmf of $y$. If $c(y)$ is
finite, then the posterior density is proper. Since we consider both proper
and improper priors on $(\beta,\tau)$, if improper priors are used, then
$c(y)$ is not necessarily finite. Conditions for the posterior propriety of
nonlinear mixed models with general link functions are given in Chen et al.
(2002) and Wang and Roy (2018b). Theorem 1, given in Section 4 of this paper,
provides easily verifiable sufficient conditions for propriety of (5) when
$\pi(\beta)\propto 1$.

Since the likelihood function $L(\beta,\tau\mid y)$ is not available in a
closed form, the posterior density for $(\beta,\tau)$ is not tractable for any
choice of priors on these parameters. Markov chain Monte Carlo (MCMC)
algorithms can be used to explore the posterior density $\pi(\beta,\tau|y)$.
Even in the absence of the random effects, MCMC algorithms are generally used
for exploring the posterior densities corresponding to the basic logistic
model or other generalized linear models (GLMs). Using the data augmentation
(DA) technique (van Dyk and Meng, 2001), in a highly cited paper, Albert and
Chib (1993) constructed a Gibbs sampler for GLMs with the probit link. Since
then there have been several attempts to construct such a DA Gibbs sampler for
the logistic model (see e.g. Holmes and Held (2006) and Frühwirth-Schnatter
and Frühwirth (2010)). Recently, using the Pólya-Gamma (PG) latent variables,
Polson et al. (2013) (denoted as PS&W hereafter) have proposed an efficient DA
Gibbs sampler for the Bayesian logistic models. A random variable $\omega$
follows a PG distribution with parameters $a>0,\,b\geq 0$, if
$\omega\overset{d}{=}(1/(2\pi^{2}))\sum_{i=1}^{\infty}g_{i}/[(i-1/2)^{2}+b^{2}/(4\pi^{2})]$,
where $g_{i}\overset{iid}{\sim}$Gamma$(a,1)$. We denote this as
$\omega\sim$PG$(a,b)$. PS&W’s DA technique can be extended to construct a
Gibbs sampler for the LLMMs. Indeed, with PG latent variables
$\omega=(\omega_{1},\omega_{2},...,\omega_{n})$, one can construct a joint
posterior density $\pi(\beta,u,\omega,\tau\mid y)$ (details are given in
Section 2) such that

$\displaystyle\int_{\mathbb{R}^{q}}\int_{\mathbb{R}_{+}^{n}}\pi(\beta,u,\omega,\tau\mid
y)d\omega\,du=\pi(\beta,\tau\mid y),$ (6)

where $\mathbb{R}_{+}=(0,\infty)$, and $\pi(\beta,\tau\mid y)$ is given in
(5). Using the conditional distributions of the joint density
$\pi(\beta,u,\omega,\tau\mid y)$, a full Gibbs sampler can be formed (The
details for this Gibbs sampler are given in Section 2.1.). It is known that
blocking parameters can improve the performance of the Gibbs sampler in terms
of reducing its operator norm (Liu et al., 1994). In general, when one or more
variables are correlated, sampling them jointly can improve the efficiency of
MCMC algorithms (Roberts and Sahu, 1997; Chib and Ramamurthy, 2010; Turek et
al., 2017). On the other hand, blocking may result in complex conditional
distributions that are not easy to sample from. For the LLMMs, it turns out
that an efficient two-block Gibbs sampler can be constructed by using the two
blocks, $\eta\equiv(\beta^{\top},u^{\top})^{\top}$ and $(\omega,\tau)$. We
derive this block Gibbs sampler in Section 2.2. Using numerical examples, we
show that blocking can lead to great gains in efficiency in Monte Carlo
estimation for the LLMMs.

The block Gibbs Markov chain is Harris ergodic. Thus, the sample (time)
averages are consistent estimators of the means with respect to the posterior
density (5). On the other hand, in practice, it is important to ascertain the
errors associated with the Monte Carlo estimates. A valid standard error for
the Monte Carlo estimate can be formed if a central limit theorem (CLT) is
available for the time average estimator (Jones and Hobert, 2001).
Establishing geometric ergodicity (GE) of the underlying Markov chain is the
most standard method for guaranteeing CLTs hold for MCMC estimators, and is
also used for consistently estimating the asymptotic variance in the CLT (Vats
et al. (2018), Vats et al. (2019)). GE of Gibbs samplers for probit and
logistic GLMs under different priors has been established in the literature
(Roy and Hobert, 2007; Choi and Hobert, 2013; Chakraborty and Khare, 2017;
Wang and Roy, 2018c). Also, GE of Gibbs samplers for probit mixed models and
normal linear mixed models under improper priors on the regression
coefficients and variance components is considered in Wang and Roy (2018b) and
Román and Hobert (2012), respectively. From Roy (2012b) it follows that Wang
and Roy’s (2018b) GE result also holds for parameter expansion for DA
algorithms for the probit mixed models. Wang and Roy (2018a) prove uniform
ergodicity of a Gibbs sampler for the LLMMs under a proper normal prior on
$\beta$ and a truncated proper prior on $\tau$. Wang and Roy (2018a) establish
a minorization condition to prove uniform ergodicity of the Gibbs sampler (see
e.g. Roberts and Rosenthal, 2004, Theorem 8). Wang and Roy’s (2018a) analysis
of the Gibbs Markov chain puts a constraint on the support of the posterior
density (5), in that the parameter $\tau$ is bounded away from zero. Here,
using a drift condition, we establish geometric convergence rates for the
block Gibbs sampler, in the case when the popular improper uniform prior is
assigned on $\beta$ and the commonly used proper gamma priors or the improper
power priors are assigned on $\tau$. In contrast to Wang and Roy (2018a), our
result does not put any restriction on the support of the variance components,
that is, the common support of $\tau_{j}$’s is the entire positive real line.

The rest of the article is organized as follows. In Section 2, we provide
details on PG data augmentation and construct the full and the block Gibbs
samplers. Section 3 contains numerical examples. These examples are used to
compare the performance of the block and the full Gibbs samplers. In Section
4, we consider geometric convergence of the block Gibbs sampler under improper
priors. Some concluding remarks are provided in Section 5. Finally, several
theoretical results along with proofs of the results appear in the appendices.

##  2 Gibbs samplers

In this section, using the PG variables, we discuss DA for the LLMMs, and
construct Gibbs samplers for (5). Following (2) and (5), the joint posterior
density for $(\beta,\tau)$ is

$\displaystyle\pi(\beta,\tau\mid
y)=\frac{\pi(\beta)\pi(\tau)}{c(y)}\int_{\mathbb{R}^{q}}\prod_{i=1}^{n}\frac{\exp\\{y_{i}(x_{i}^{\top}\beta+z_{i}^{\top}u)\\}}{1+\exp(x_{i}^{\top}\beta+z_{i}^{\top}u)}\phi_{q}(u;0,D(\tau)^{-1})du.$

By Theorem 1 in Polson et al. (2013)

$\displaystyle\pi(\beta,\tau\mid
y)=\int_{\mathbb{R}^{q}}\int_{\mathbb{R}_{+}^{n}}\bigg{[}\prod_{i=1}^{n}$
$\displaystyle\frac{\exp\\{k_{i}(x_{i}^{\top}\beta+z_{i}^{\top}u)-\omega_{i}(x_{i}^{\top}\beta+z_{i}^{\top}u)^{2}/2\\}}{2}p(\omega_{i})\bigg{]}d\omega\;$
$\displaystyle\times\phi_{q}(u;0,D(\tau)^{-1})du\times\frac{\pi(\beta)\pi(\tau)}{c(y)},$
(7)

where $\omega=(\omega_{1},\omega_{2},...\omega_{n})$,
$k_{i}=y_{i}-1/2,\,i=1,...,n$ and $p(\omega_{i})$ is the pdf of
$\omega_{i}\sim$PG$(1,0)$ given by,

$p(\omega_{i})=\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi\omega_{i}^{3}}}\exp\bigg{[}-\frac{(2\ell+1)^{2}}{8\omega_{i}}\bigg{]},\,\omega_{i}>0.$
(8)

We now define the joint posterior density of $\beta,u,\omega,\tau$ given $y$
mentioned in (6) as

$\displaystyle\pi(\beta,u,\omega,\tau\mid y)$
$\displaystyle\propto\bigg{[}\prod_{i=1}^{n}\exp\\{k_{i}(x_{i}^{\top}\beta+z_{i}^{\top}u)-\omega_{i}(x_{i}^{\top}\beta+z_{i}^{\top}u)^{2}/2\\}p(\omega_{i})\bigg{]}\pi(\beta)\pi(\tau)$
$\displaystyle\quad\times\phi_{q}(u;0,D(\tau)^{-1})$
$\displaystyle=\bigg{[}\prod_{i=1}^{n}\exp\\{k_{i}(x_{i}^{\top}\beta+z_{i}^{\top}u)-\omega_{i}(x_{i}^{\top}\beta+z_{i}^{\top}u)^{2}/2\\}p(\omega_{i})\bigg{]}$
$\displaystyle\quad\times\phi_{q}(u;0,D(\tau)^{-1})\times\exp\Big{[}-\frac{1}{2}(\beta-\mu_{0})^{\top}Q(\beta-\mu_{0})\Big{]}$
$\displaystyle\quad\times\prod_{j=1}^{r}\tau_{j}^{a_{j}-1}\exp(-b_{j}\tau_{j}),$
(9)

where (2) follows from the priors on $\beta$ and $\tau$ given in (3) and (4),
respectively.

###  2.1 A full Gibbs sampler

Let $\Omega(\omega)$ be the $n\times n$ diagonal matrix with the $i^{th}$
diagonal element $\omega_{i}$. Let $M=[XZ]$, and $m_{i}^{\top}$ indicates the
$i^{th}$ row of $M$ for $i=1,\dots,n$. Let
$\kappa=(k_{1},k_{2},...,k_{n})^{\top}$. We begin with deriving the
conditional densities required for the full Gibbs (FG) sampler. Based on (2),
the conditional density of $\beta$ given $u,\omega,\tau,y$ is

$\displaystyle\pi(\beta\mid u,\omega,\tau,y)$
$\displaystyle\propto\prod_{i=1}^{n}\exp\Big{[}k_{i}x_{i}^{\top}\beta-\omega_{i}(x_{i}^{\top}\beta)^{2}/2-\omega_{i}x_{i}^{\top}\beta
z_{i}^{\top}u\Big{]}$
$\displaystyle\quad\times\exp\Big{[}-\frac{1}{2}(\beta-\mu_{0})^{\top}Q(\beta-\mu_{0})\Big{]}$
$\displaystyle\propto\exp\Big{[}-\frac{1}{2}\beta^{\top}(X^{\top}\Omega(\omega)X+Q)\beta+\beta^{\top}(X^{\top}\kappa+Q\mu_{0}-X^{\top}\Omega(\omega)Zu)\Big{]}.$

Hence,

$\displaystyle\beta\mid u,\omega,\tau,y\sim
N((X^{\top}\Omega(\omega)X+Q)^{-1}(X^{\top}\kappa+Q\mu_{0}-X^{\top}\Omega(\omega)Zu),(X^{\top}\Omega(\omega)X+Q)^{-1}).$
(10)

Also from (2), the conditional density of $u$ given $\beta,\omega,\tau,y$ is

$\displaystyle\pi(u\mid\beta,\omega,\tau,y)$
$\displaystyle\propto\prod_{i=1}^{n}\exp\Big{\\{}k_{i}z_{i}^{\top}u-\frac{\omega_{i}}{2}\Big{[}(z_{i}^{\top}u)^{2}+2z_{i}^{\top}ux_{i}^{\top}\beta\Big{]}\Big{\\}}\exp\Big{[}-\frac{1}{2}u^{\top}D(\tau)u\Big{]}$
$\displaystyle=\exp\Big{[}-\frac{1}{2}u^{\top}(Z^{\top}\Omega(\omega)Z+D(\tau))u+u^{\top}(Z^{\top}\kappa-Z^{\top}\Omega(\omega)X\beta)\Big{]}.$

Thus, it follows that

$\displaystyle u\mid\beta,\omega,\tau,y\sim
N((Z^{\top}\Omega(\omega)Z+D(\tau))^{-1}(Z^{\top}\kappa-Z^{\top}\Omega(\omega)X\beta),(Z^{\top}\Omega(\omega)Z+D(\tau))^{-1}).$
(11)

Also from (2), the conditional density of $\omega$ and $\tau$ given $\eta$ and
$y$ is as follows

$\displaystyle\pi(\omega,\tau\mid\eta,y)$
$\displaystyle\propto\prod_{i=1}^{n}\exp(-\omega_{i}(m_{i}^{\top}\eta)^{2}/2)p(\omega_{i})\left|D(\tau)\right|^{\frac{1}{2}}\exp(-u^{\top}D(\tau)u/2)$
$\displaystyle\quad\times\prod_{j=1}^{r}\tau_{j}^{a_{j}-1}\exp(-b_{j}\tau_{j}),$
(12)

where $\left|D(\tau)\right|$ is the determinant of $D(\tau)$. From the above,
we see that $\omega_{i}$’s, $i=1,...,n$ are conditionally independent given
$(\eta,y)$. Also, $(\omega,\tau)$ are conditionally mutually independent given
$(\eta,y)$. The conditional density for $\omega_{i}$ is

$\displaystyle\pi(\omega_{i}\mid\eta,y)$
$\displaystyle\propto\exp(-\omega_{i}(m_{i}^{\top}\eta)^{2}/2)p(\omega_{i})$
$\displaystyle=\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi\omega_{i}^{3}}}\exp\Big{(}-\frac{(2\ell+1)^{2}}{8\omega_{i}}-\frac{\omega_{i}(m_{i}^{\top}\eta)^{2}}{2}\Big{)},$
(13)

where the equality follows from (8). From Wang and Roy (2018c), the pdf for
PG$(a,b),a>0,\,b\in\mathbb{R}$ is

$\displaystyle p(x\mid
a,b)=\Big{[}\cosh\big{(}\frac{b}{2}\big{)}\Big{]}^{a}\frac{2^{a-1}}{\Gamma{(a)}}\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{\Gamma{(\ell+a)}}{\Gamma{(\ell+1)}}\frac{(2\ell+a)}{\sqrt{2\pi
x^{3}}}\exp\Big{(}-\frac{(2\ell+a)^{2}}{8x}-\frac{xb^{2}}{2}\Big{)},$

for $x>0,$ where the hyperbolic cosine function $\cosh(t)=(e^{t}+e^{-t})/2$.
Hence, from (2.1) we have

$\displaystyle\omega_{i}\mid\eta,y\overset{ind}{\sim}PG\big{(}1,\mathinner{\\!\left\lvert
m_{i}^{\top}\eta\right\rvert}\big{)},\,i=1,....n.$ (14)

From (2.1), the conditional density for $\tau_{j}$ is given by

$\displaystyle\pi(\tau_{j}\mid\eta,y)\propto\tau_{j}^{q_{j}/2+a_{j}-1}\exp\big{[}-\tau_{j}(b_{j}+u_{j}^{\top}u_{j}/2)\big{]},\,j=1,...,r.$
(15)

Thus, we have

$\displaystyle\tau_{j}\mid\eta,y\overset{ind}{\sim}\text{Gamma}(a_{j}+q_{j}/2,b_{j}+u_{j}^{\top}u_{j}/2),\,j=1,...,r.$

We will allow $a_{j}>-q_{j}/2$, and $b_{j}\geq 0$. The case where $b_{j}=0$ is
more complicated and will be discussed at length in Section 4.

Let $(\beta^{(m)},u^{m},\omega^{(m)},\tau^{(m)})$ denote the $m^{th}$ element
for $(\beta,u,\omega,\tau)$ in the FG chain. Thus, a single iteration of the
full Gibbs sampler
$\\{\beta^{(m)},u^{m},\omega^{(m)},\tau^{(m)}\\}_{m=0}^{\infty}$ consists of
the following four steps:

$\bm{Algorithm\,1}$ The (m+1)st iteration of the full Gibbs sampler
---
1\. Draw $\tau_{j}^{(m+1)}\overset{ind}{\sim}$
Gamma$(a_{j}+q_{j}/2,b_{j}+u_{j}^{\top}u_{j}/2),\,j=1,...,r$ with $u=u^{(m)}$.
2\. Draw
$\omega_{i}^{(m+1)}\overset{ind}{\sim}PG\big{(}1,\mathinner{\\!\left\lvert
m_{i}^{\top}\eta^{(m)}\right\rvert}\big{)},\,i=1,....n$.
3\. Draw $u^{(m+1)}\sim\eqref{eq:pu}$ with
$\tau=\tau^{(m+1)},\,\omega=\omega^{(m+1)}$, and $\beta=\beta^{(m)}$.
4\. Draw $\beta^{(m+1)}\sim\eqref{eq:pbeta}$ with $\omega=\omega^{(m+1)}$.

###  2.2 A two-block Gibbs sampler

In this section, we construct a block Gibbs (BG) sampler for (5). As mentioned
before, $M=[XZ]$ with the $i^{th}$ row $m_{i}^{\top}$ for $i=1,\dots,n$. Note
that $x_{i}^{\top}\beta+z_{i}^{\top}u=m_{i}^{\top}\eta$. From (2), the
conditional density of $\eta$ given $\omega,\tau,y$ is given by

$\displaystyle\pi(\eta\mid\omega,\tau,y)$
$\displaystyle\propto\prod_{i=1}^{n}\exp\big{[}k_{i}m_{i}^{\top}\eta-\omega_{i}(m_{i}^{\top}\eta)^{2}/2\big{]}\exp\Big{[}-\frac{1}{2}u^{\top}D(\tau)u\Big{]}$
$\displaystyle\quad\times\exp\Big{[}-\frac{1}{2}(\beta-\mu_{0})^{\top}Q(\beta-\mu_{0})\Big{]}$
$\displaystyle\propto\exp\Big{[}-\frac{1}{2}(\eta-\Sigma(M^{\top}\kappa+b))^{\top}\Sigma^{-1}(\eta-\Sigma(M^{\top}\kappa+b))\Big{]},$
(16)

where $\Sigma^{-1}=M^{\top}\Omega(\omega)M+A(\tau)$,

$b_{(p+q)\times 1}=\begin{pmatrix}Q\mu_{0}\\\ 0_{q\times
1}\end{pmatrix}\;\mbox{and}\;A(\tau)_{(p+q)(p+q)}=\begin{pmatrix}Q&0\\\
0&D(\tau)\end{pmatrix}.$

Hence,

$\displaystyle\eta\mid\omega,\tau,y\sim
N((M^{\top}\Omega(\omega)M+A(\tau))^{-1}(M^{\top}\kappa+b),(M^{\top}\Omega(\omega)M+A(\tau))^{-1}).$
(17)

In the FG sampler in Section 2.1, $\tau,\,\omega,\,u$ and $\beta$ are drawn
sequentially, whereas, in this section, we show that the conditional
distribution of $\eta$ given $\omega,\,\tau,\,y$ is normal. From (2.1), we can
see, conditional on $(\eta,y)$, $\omega$ and $\tau$ are independent. Thus,
$\tau$ and $\omega$ can be drawn jointly as a block and we have a two-block
Gibbs sampler.

Let $\eta^{(m)},\omega^{(m)},$ and $\tau^{(m)}$ denote the values of
$\eta,\omega,$ and $\tau$, respectively, in the $m^{th}$ iteration of the BG
sampler. A single iteration of the block Gibbs sampler
$\\{\eta^{(m)},\omega^{(m)},\tau^{(m)}\\}_{m=0}^{\infty}$ consists of the
following two steps:

$\bm{Algorithm\,2}$ The (m+1)st iteration of the two-block Gibbs sampler
---
1\. Draw $\tau_{j}^{(m+1)}\overset{ind}{\sim}$
Gamma$(a_{j}+q_{j}/2,b_{j}+u_{j}^{\top}u_{j}/2),\,j=1,...,r$ with $u=u^{(m)}$,
and independently draw
$\omega_{i}^{(m+1)}\overset{ind}{\sim}PG\big{(}1,\mathinner{\\!\left\lvert
m_{i}^{\top}\eta^{(m)}\right\rvert}\big{)},\,i=1,....n$.
2\. Draw $\eta^{(m+1)}\sim\eqref{eq:distrieta}$ with $\tau=\tau^{(m+1)}$ and
$\omega=\omega^{(m+1)}$.

The conditional distribution of $\eta$ in the BG sampler and the conditional
distributions of $\beta$ and $u$ in the FG sampler are all normal
distributions of the form $N(S^{-1}t,S^{-1})$ for some matrix $S$ and a vector
$t$. Note that, for the conditional distribution of $\eta$, $S$ is a
$(p+q)\times(p+q)$ matrix, whereas for $\beta$ and $u$ this is a $p\times p$
and $q\times q$ matrix, respectively. Thus, a naive method of drawing from
$N(S^{-1}t,S^{-1})$ is inefficient especially if $p$ and/or $q$ is large as it
involves calculating inverse of the matrix $S$. Here, we use a known method of
drawing from $N(S^{-1}t,S^{-1})$ that does not require computing $S^{-1}$. The
method is as follows:

Algorithm for drawing from $N(S^{-1}t,S^{-1})$
---
1\. Let $S=LL^{\top}$ be the Cholesky decomposition of $S$.
2\. Solve $Lw=t$.
3\. Draw $z\sim N(0,{\mathrm{I}}_{k})$ where $k$ is the dimension of $S$.
4\. Solve $L^{\top}x=w+z$. Then $x\sim N(S^{-1}t,S^{-1})$.

##  3 A real data example

We consider the student performance data set from Cortez and Silva (2008).
This data set includes $n=649$ observations and $33$ variables including
several categorical variables. As in Cortez and Silva (2008), the binary
response is defined as $1$ if the final grade is greater than or equal to
$10$, otherwise, it is defined as $0$. Recall that $p$ denotes the number of
columns for the design matrix $X$. Also, note that, categorical variables are
incorporated into the LLMM as sets of dichotomous variables through what is
known as dummy coding. We consider different subsets of variables while
fitting the LLMM to compare the BG and FG samplers for different dimensions.
In particular, we consider $p=3,7,23$, including an intercept term. We also
keep one random effect ‘school’ with $2$ levels in the LLMM. As mentioned in
Chib and Ramamurthy (2010), a key principle to block sampling in MCMC is that
‘parameters in different blocks are not strongly correlated whereas those
within a block are’ (see also Roberts and Sahu, 1997; Turek et al., 2017). As
we observe later in this Section, when $p$ varies, the average absolute
(posterior) correlations between $\beta$ and $u$ change greatly. The specific
values of $p$ we choose here are irrelevant for the general conclusions of
this section. Indeed, for some other $p$ values, close to the ones we choose
here, the general pattern of the different empirical measures remains the
same.

We analyze the data set by fitting the LLMM with a proper normal prior (3) on
$\beta$ with $\mu_{0}=0$ and $Q=0.001{\mathrm{I}}_{3}$. Also, we use a proper
Gamma prior (4) on $\tau_{1}$ where the (prior) mean and variance of
$\tau_{1}$ are $1.2$ and $100$, respectively ($a_{1}=0.0144$ and
$b_{1}=0.012$). We ran the BG sampler for $m=120,000$ iterations, starting at
an initial value $\eta^{(0)}=(\beta^{(0)},u^{(0)})$ with a burn-in of
$B=20,000$ iterations. Here, $\beta^{(0)}$ is the estimate of $\beta$ obtained
by fitting a logistic linear model without any random effect. For $p=3,7$, the
initial value $u^{(0)}$ is a sample drawn from
$N(0,(1/\tau_{1}^{(0)}){\mathrm{I}}_{2})$, where $1/\tau_{1}^{(0)}$ is the
estimate of the random effect variance component obtained from the R package
lme4. For $p=23$, $1/\tau_{1}^{(0)}$ is the estimate of the random effect
variance component obtained from $p=7$ as lme4 did not run successfully in the
case of $p=23$. The FG sampler was also run for $m=120,000$ iterations with a
burn-in of $B=20,000$ iterations.

The BG and FG samplers are compared using the lag $k$ autocorrelation function
(ACF) values $k=1,...,5$, the effective sample size (ESS) and the multivariate
ESS (mESS) (See Roy (2020) for a simple introduction to these convergence
diagnostic measures.). The ESS and mESS are calculated using the R package
mcmcse. We also compute the mean squared jumps (MSJ) defined as
$\sum_{i=B+1}^{m}\mathinner{\\!\left\lVert\beta^{(i+1)}-\beta^{(i)}\right\rVert}^{2}/(m-B)$
for the $\beta$ variable, and similarly for the other variables. Here,
$\mathinner{\\!\left\lVert\cdot\right\rVert}$ denotes the Euclidean norm.
Tables 1–3 provide the values of the ACF for $\beta_{0},\beta_{1},\beta_{2}$
and $\tau_{1}$ for the BG and FG samplers, as $p$ varies. Better performance
of the BG sampler compared to the FG sampler is observed from its mostly
smaller ACF values. Table 4 provides the ESS values of the intercept
parameter, first two regression coefficients and $\tau_{1}$. It also gives the
mESS values for $u$ and $(\beta,\tau_{1})$, as $p$ varies. Again, better
efficiency of the BG sampler compared to the FG sampler is observed from its
generally larger ESS and mESS values. We calculate the average of the absolute
(posterior) correlations between the coordinates of the $\beta$ vector and
those of the $u$ vector computed based on the BG samples mentioned before. For
$p=3,7,$ and 23, these values are 0.2601, 0.1143, and 0.0206, respectively.
From the mESS values of the parameters of interest, namely $(\beta,\tau_{1})$,
given in Table 4, we see that the increase in the efficiency of the BG sampler
compared to the FG sampler is higher for smaller $p$ values. From Table 5, it
can be seen that the BG sampler leads to higher MSJ values than the FG sampler
with the exception of $\tau_{1}$ in some cases. Thus, Table 5 also
corroborates better mixing of the BG sampler than the FG sampler. So, in
practice, the BG sampler can provide significant gains compared to the FG
sampler.

Finally, Table 4 also provides the time normalized efficiency (ESS and mESS
values per second) for the two samplers, as $p$ varies. From this Table, we
see that the BG sampler has always resulted in larger ESS and mESS values per
second than the FG sampler. Indeed, the BG sampler results in higher time
normalized ESS values even in the cases when the FG sampler performs better in
terms of the ESS (See e.g. ESS ($\beta_{1}$) for $p=7$.). Recall that in every
iteration, the BG sampler makes a draw from a $(p+q)$ dimensional normal
distribution, whereas the FG sampler draws from a $q$ dimensional normal
distribution and then a $p$ dimensional normal distribution. Other draws are
the same for both the BG and FG samplers. Using the Cholesky decomposition
method mentioned in Section 2.2, we observe that for all values of $p$
considered here, the BG sampler takes less time than the FG sampler to
complete a certain number of iterations. On the other hand, when $(p+q)$ takes
much larger values, the BG sampler takes more running time than the FG
sampler.

Table 1: ACF for the BG and FG samplers for the student performance data with $p=3$ Parameter | Sampler | lag 1 | lag 2 | lag 3 | lag 4 | lag 5
---|---|---|---|---|---|---
$\beta_{0}$ | BG | 0.434 | 0.385 | 0.359 | 0.333 | 0.319
| FG | 0.985 | 0.974 | 0.964 | 0.956 | 0.948
$\beta_{1}$ | BG | 0.597 | 0.380 | 0.258 | 0.192 | 0.153
| FG | 0.613 | 0.402 | 0.282 | 0.212 | 0.173
$\beta_{2}$ | BG | 0.836 | 0.734 | 0.667 | 0.621 | 0.586
| FG | 0.838 | 0.739 | 0.671 | 0.623 | 0.587
$\tau_{1}$ | BG | 0.374 | 0.212 | 0.134 | 0.091 | 0.071
| FG | 0.405 | 0.285 | 0.224 | 0.186 | 0.155
Table 2: ACF for the BG and FG samplers for the student performance data with $p=7$ Parameter | Sampler | lag 1 | lag 2 | lag 3 | lag 4 | lag 5
---|---|---|---|---|---|---
$\beta_{0}$ | BG | 0.463 | 0.409 | 0.365 | 0.340 | 0.320
| FG | 0.924 | 0.867 | 0.821 | 0.781 | 0.747
$\beta_{1}$ | BG | 0.436 | 0.191 | 0.087 | 0.036 | 0.018
| FG | 0.437 | 0.193 | 0.088 | 0.038 | 0.017
$\beta_{2}$ | BG | 0.419 | 0.185 | 0.085 | 0.047 | 0.031
| FG | 0.420 | 0.187 | 0.090 | 0.053 | 0.032
$\tau_{1}$ | BG | 0.379 | 0.216 | 0.142 | 0.098 | 0.063
| FG | 0.372 | 0.244 | 0.190 | 0.150 | 0.122
Table 3: ACF for the BG and FG samplers for the student performance data with $p=23$ Parameter | Sampler | lag 1 | lag 2 | lag 3 | lag 4 | lag 5
---|---|---|---|---|---|---
$\beta_{0}$ | BG | 0.870 | 0.811 | 0.761 | 0.714 | 0.669
| FG | 0.933 | 0.874 | 0.820 | 0.770 | 0.723
$\beta_{1}$ | BG | 0.637 | 0.437 | 0.323 | 0.256 | 0.213
| FG | 0.655 | 0.463 | 0.348 | 0.278 | 0.235
$\beta_{2}$ | BG | 0.885 | 0.807 | 0.753 | 0.713 | 0.681
| FG | 0.880 | 0.801 | 0.744 | 0.702 | 0.669
$\tau_{1}$ | BG | 0.384 | 0.228 | 0.147 | 0.098 | 0.070
| FG | 0.395 | 0.263 | 0.198 | 0.164 | 0.141
Table 4: Multivariate ESS or ESS for the BG and FG samplers for the student performance data with $p=3,7,$ and $23$. The numbers inside the parentheses are the corresponding values per second. The column name MC stands for Markov chain. $p$ | MC | mESS ($\beta\;\tau_{1}$) | mESS ($\beta$) | ESS ($\beta_{0}$) | ESS ($\beta_{1}$) | ESS ($\beta_{2}$) | mESS($u$) | ESS ($\tau_{1}$)
---|---|---|---|---|---|---|---|---
$3$ | BG | 19,012 | 15,979 | 4,229 | 9,268 | 2,586 | 34,623 | 31,688
(321) | (270) | (71) | (156) | (44) | (585) | (535)
FG | 1,539 | 978 | 40 | 8,089 | 2,560 | 76 | 2,037
(22) | (14) | (1) | (115) | (36) | (1) | (29)
$7$ | BG | 27,474 | 26,702 | 4,494 | 36,015 | 31,294 | 34,844 | 32,334
(426) | (414) | (70) | (559) | (486) | (541) | (502)
FG | 13,533 | 12,793 | 1,894 | 38,572 | 29,826 | 1,931 | 11,401
(179) | (169) | (25) | (509) | (394) | (25) | (150)
$23$ | BG | 23,068 | 22,966 | 3,031 | 7,272 | 2,130 | 31,529 | 28,561
(295) | (294) | (39) | (93) | (27) | (404) | (366)
FG | 18,016 | 17,770 | 3,040 | 5,912 | 1,554 | 1,022 | 7,911
(187) | (184) | (32) | (61) | (16) | (11) | (82)
Table 5: Mean squared jumps for the BG and FG samplers for the student performance data with $p=3,7,$ and $23$ p | BG | FG
---|---|---
| $\beta$ | $u$ | $\tau$ | $\beta$ | $u$ | $\tau$
3 | 12.94 | 24.41 | 1319.13 | 0.73 | 0.10 | 1071.40
7 | 15.52 | 24.97 | 1192.74 | 3.03 | 0.10 | 1309.20
23 | 90.19 | 29.68 | 1240.48 | 74.84 | 0.11 | 1249.00

##  4 Geometric ergodicity of the block Gibbs sampler

We begin this section with a discussion on the conditional density
$\pi(\tau\mid\eta,y)$. Since we allow the prior rate parameter $b_{j}$ for
$\tau_{j}$ to be zero, define
$A=\\{j\in\\{1,2,...,r\\}\mathrel{\mathop{\mathchar 58\relax}}b_{j}=0\\}$.
Recall from (15) that $\tau_{j}\mid\eta,y\overset{ind}{\sim}$
Gamma$(a_{j}+q_{j}/2,b_{j}+u_{j}^{\top}u_{j}/2),\,j=1,...,r$ when
$a_{j}+q_{j}/2>0$ and $b_{j}+u_{j}^{\top}u_{j}/2>0$. The density
$\pi(\tau\mid\eta,y)=\prod_{j=1}^{r}\pi(\tau_{j}\mid\eta,y)$ is not defined
when $A$ is not empty and $\|u_{j}\|=0$ for $j\in A$. Let
$N=\\{\eta\in\mathbb{R}^{p+q},\prod_{j\in A}\|u_{j}\|=0\\}$. The fact that
$\pi(\tau\mid\eta,y)$ is not defined on $N$ is irrelevant for simulating the
BG sampler as $N$ is a measure zero set with respect to the Lebesgue measure
on $\mathbb{R}^{p+q}$. But, for a theoretical analysis of the BG chain,
$\pi(\tau\mid\eta,y)$ needs to be defined for all $\eta\in\mathbb{R}^{p+q}$.
For all $\eta\in\mathbb{R}^{p+q}$, we define

$\pi(\tau\mid\eta,y)=\begin{cases}\prod_{j=1}^{r}f_{G}(\tau_{j},a_{j}+\frac{q_{j}}{2},b_{j}+\frac{1}{2}u_{j}^{\top}u_{j})&\text{if}\quad\eta\not\in
N\\\ \prod_{j=1}^{r}f_{G}(\tau_{j},1,1)&\text{if}\quad\eta\in N\end{cases}.$
(18)

Here, $f_{G}(x,a,b)$ denotes the pdf of a gamma random variable with the shape
parameter $a$, the rate parameter $b$, and evaluated at $x$. Thus,
$f_{G}(x,a,b)=(b^{a}/\Gamma(a))x^{a-1}\exp(-bx)$. The Markov transition
density (Mtd) of the BG chain
$\\{\eta^{(m)},\omega^{(m)},\tau^{(m)}\\}_{m=0}^{\infty}$ is

$\displaystyle
k(\eta,\omega,\tau\mid\eta^{\prime},\omega^{\prime},\tau^{\prime})$
$\displaystyle=\pi(\eta\mid\omega,\tau,y)\pi(\omega,\tau\mid\eta^{\prime},y)$
$\displaystyle=\pi(\eta\mid\omega,\tau,y)\bigg{[}\prod_{i=1}^{n}\pi(\omega_{i}\mid\eta^{\prime},y)\bigg{]}\pi(\tau\mid\eta^{\prime},y),$
(19)

where the conditional densities
$\pi(\eta\mid\omega,\tau,y),\pi(\omega_{i}\mid\eta^{\prime},y),$ and
$\pi(\tau\mid\eta^{\prime},y)$ on the right side of (4) are given in (2.2),
(2.1), and (18), respectively. It is easy to see that the joint density (2) is
the invariant density of $k$, and $k$ is $\varphi$-irreducible. Thus, if (2)
is a proper density, that is, if $c(y)<\infty$ in (5), then the BG chain
$\\{\eta^{(m)},\omega^{(m)},\tau^{(m)}\\}_{m=0}^{\infty}$ is Harris ergodic
(Meyn and Tweedie, 1993, Chap 10), and hence, it can be used to consistently
estimate means with respect to (2). Let
$S=\mathbb{R}^{p+q}\times\mathbb{R}_{+}^{n}\times\mathbb{R}_{+}^{r}$. In fact,
if $g\mathrel{\mathop{\mathchar 58\relax}}S\rightarrow\mathbb{R}$ is
integrable with respect to (2), that is, if

$\displaystyle{\mathrm{E}}_{\pi}|g(\eta,\omega,\tau)|\mathrel{\mathop{\mathchar
58\relax}}=\int_{S}\mathinner{\\!\left\lvert
g(\eta,\omega,\tau)\right\rvert}\pi(\beta,u,\omega,\tau\mid
y)d\eta\,d\omega\,d\tau<\infty,$

then $\overline{g}_{m}\mathrel{\mathop{\mathchar
58\relax}}=\sum_{i=0}^{m-1}g(\eta^{(i)},\omega^{(i)},\tau^{(i)})/m\rightarrow{\mathrm{E}}_{\pi}\,g$
almost surely as $m\rightarrow\infty$. On the other hand, even when
${\mathrm{E}}_{\pi}\,g^{2}<\infty$, Harris ergodicity of $k$ does not
guarantee a CLT holds for $\overline{g}_{m}$, which is used to obtain valid
standard errors of $\overline{g}_{m}$. We say a CLT for $\overline{g}_{m}$
exists if
$\sqrt{m}(\overline{g}_{m}-{\mathrm{E}}_{\pi}\,g)\stackrel{{\scriptstyle
d}}{{\rightarrow}}N(0,\sigma_{g}^{2})\;as\;m\rightarrow\infty$ for some
$\sigma^{2}_{g}\in(0,\infty)$. Certain convergence rates of the BG chain, as
we explain next, ensure a CLT holds for $\overline{g}_{m}$.

Let $\mathcal{B}(S)$ denotes the Borel $\sigma$-algebra of $S$. Let
$K^{(m)}\mathrel{\mathop{\mathchar
58\relax}}S\times\mathcal{B}(S)\rightarrow[0,1]$ denotes the $m$-step Markov
transition function (Mtf) corresponding to the Mtd (4), that is,

$\displaystyle
K^{(m)}((\eta^{\prime},\omega^{\prime},\tau^{\prime}),B)=P((\eta^{(m+j)},\omega^{(m+j)},\tau^{(m+j)})\in
B\,|\,(\eta^{(j)},\omega^{(j)},\tau^{(j)})=(\eta^{\prime},\omega^{\prime},\tau^{\prime})),$

for any $j\in\\{1,2,\dots\\}$ and for any measurable set $B\in\mathcal{B}(S)$.
The BG chain is geometrically ergodic if there exist a function
$H\mathrel{\mathop{\mathchar 58\relax}}S\rightarrow[0,\infty)$ and a constant
$\rho\in(0,1)$ such that for all $m=0,1,2,...$,

$\displaystyle\|K^{m}((\eta^{\prime},\omega^{\prime},\tau^{\prime}),\cdot)-\Pi(\cdot)\|_{{\mathrm{TV}}}\mathrel{\mathop{\mathchar
58\relax}}=\sup_{B\in\mathcal{B}(S)}\|K^{m}((\eta^{\prime},\omega^{\prime},\tau^{\prime}),B)-\Pi(B)\|\leq
H(\eta^{\prime},\omega^{\prime},\tau^{\prime})\rho^{m},$ (20)

where $\Pi(\cdot)$ denotes the probability measure corresponding to the joint
posterior density (2), and $\|\cdot\|_{{\mathrm{TV}}}$ denotes the total
variation norm. Harris ergodicity of $k$ implies the TV norm in (20)
$\downarrow 0$ as $m\rightarrow\infty$, but does not ascertain any rate at
which this convergence takes place. On the other hand, (20) guarantees a CLT
for $\overline{g}_{m}$ if ${\mathrm{E}}_{\pi}\,|g|^{2+\delta}<\infty$ for some
$\delta>0$ (Roberts and Rosenthal, 2004). If (20) holds, it also implies that
consistent batch means and a spectral variance estimator
$\hat{\sigma}^{2}_{g}$ of $\sigma^{2}_{g}$ are available (Vats et al. (2018),
Vats et al. (2019)), and thus a valid standard error (SE)
$\hat{\sigma}_{g}/\sqrt{m}$ for $\overline{g}_{m}$ can be calculated. An
advantage of being able to calculate a valid SE is that it can be used to
decide ‘when to stop’ running the BG chain (Roy, 2020).

An important property that we are going to use in this article is that the
marginal sequences $\\{\eta^{(m)}\\}_{m=0}^{\infty}$,
$\\{(\omega^{(m)},\tau^{(m)})\\}_{m=0}^{\infty}$ of the BG chain
$\\{\eta^{(m)},\omega^{(m)},\tau^{(m)}\\}_{m=0}^{\infty}$ are themselves
Markov chains, and either all three chains are geometrically ergodic or none
of them (Liu et al., 1994; Roberts and Rosenthal, 2001). Thus, we are free to
analyze any of these chains to study their geometric convergence properties.
Indeed, here we analyze the $\\{\eta^{(m)}\\}_{m=0}^{\infty}$ marginal chain.

We denote the Markov chain $\\{\eta^{(m)}\\}_{m=0}^{\infty}$ on
$\mathbb{R}^{p+q}$ by $\Psi$ while the Markov chain
$\\{\eta^{(m)}\\}_{m=0}^{\infty}$ on $\mathbb{R}^{p+q}\backslash N$ is denoted
by $\tilde{\Psi}$. From (4) it follows that the Mtd of the $\Psi$ chain is

$\displaystyle\tilde{k}(\eta\mid\eta^{\prime})=\int_{\mathbb{R}_{+}^{r}}\int_{\mathbb{R}_{+}^{n}}\pi(\eta\mid\omega,\tau,y)\pi(\omega,\tau\mid\eta^{\prime},y)d\omega\,d\,\tau,$
(21)

We can verify that $\tilde{k}(\eta\mid\eta^{\prime})\pi(\eta^{\prime}\mid
y)=\tilde{k}(\eta^{\prime}\mid\eta)\pi(\eta\mid y)$ for all
$\eta,\eta^{\prime}\in\mathbb{R}^{p+q}$ where $\pi(\eta\mid
y)=\int_{\mathbb{R}_{+}^{r}}\int_{\mathbb{R}_{+}^{n}}\pi(\eta,\omega,\tau\mid
y)d\omega\,d\,\tau$ is the $\eta$ marginal density of (2). Hence, (21) is
reversible with respect to $\pi(\eta\mid y)$, and thus, $\pi(\eta\mid y)$ is
the invariant density for the Markov chain $\\{\eta^{(m)}\\}_{m=0}^{\infty}$.
Also, since $\\{\eta^{(m)}\\}_{m=0}^{\infty}$ is reversible, GE of the chain
implies that CLTs hold for all square integrable functions with respect to
$\pi(\eta\mid y)$ (Roberts and Rosenthal, 1997). We first establish GE of the
$\tilde{\Psi}$ chain. As explained in the proof of Theorem $1$, GE of
$\tilde{\Psi}$ implies that of $\Psi$.

######  Theorem 1.

If $\pi(\beta)\propto 1$, that is, if $Q=0$ in (3), the Markov chain
underlying the block Gibbs sampler is geometrically ergodic if the following
conditions hold:

  1. 1.

$a_{j}<b_{j}=0$ or $b_{j}>0$ for $j=1,...,r$;

  2. 2.

$a_{j}+q_{j}/2>0$ for $j=1,...,r$;

  3. 3.

$M$ has full rank;

  4. 4.

There exists a positive vector $e>0$ such that $e^{\prime}M^{*}=0$ where
$M^{*}$ is an $n\times(p+q)$ matrix with $i$th row $c_{i}m_{i}^{\top}$, where
$c_{i}=1-2y_{i}$, $i=1,\dots,n$.

The proof of Theorem 1 is given in the Appendix C. The condition $4$ can be
easily checked by an optimization method presented in Roy and Hobert (2007).

######  Remark 1.

The conditions in Theorem $1$ are the same as the conditions assumed in Wang
and Roy’s (2018b) Theorem $2$ that establishes GE of Gibbs samplers for the
probit linear mixed model.

######  Remark 2.

As mentioned before, Wang and Roy (2018a) analyzed the PG sampler for LLMMs
with proper normal priors on $\beta$ and a truncated gamma prior on $\tau$.
Wang and Roy’s (2018a) proof established uniform ergodicity which is stronger
than geometric ergodicity, but they assumed a stronger prior on $\tau$.
Indeed, their proof involving a minorization condition requires that the
support of $\tau$ is bounded away from zero. Our analysis of the BG Markov
chain does not entail any minorization condition, and does not put any
restriction, other than being positive, on the support of the variance
components.

##  5 Conclusion

In this article, we consider an efficient block Gibbs sampler based on the
Pólya-Gamma DA (Polson et al., 2013) for one of the most widely used
statistical models, namely the LLMMs. Through numerical examples, we observe
that blocking can improve performance of the Pólya-Gamma Gibbs samplers. We
hope that the article will encourage development and use of efficient blocking
strategies for Monte Carlo estimation of other GLMMs, including spatial GLMMs
where MCMC algorithms are known to suffer from slow mixing as noted in
Evangelou and Roy (2019).

Undertaking a Foster-Lyapunov drift analysis, we establish CLTs for the BG
sampler based MCMC estimators under the improper uniform prior on the
regression coefficients and improper or proper priors on the variance
components. These theoretical results are crucial for obtaining standard
errors for MCMC estimates of posterior means. In the process of our proof for
demonstrating CLTs for the BG sampler, we also establish some general results
on the Pólya-Gamma distribution. A potential future problem is to construct
and study block Gibbs samplers for other GLMMs, including the mixed models
with the robit link (Roy, 2012a).

## Acknowledgment

The authors thank the editor and two anonymous reviewers for several helpful
comments and suggestions that led to an improved revision of the paper.

Appendices

##  Appendix A Some useful results

Recall from Section 2.2 that if $Q=0$, that is, if $\pi(\beta)\propto 1$, then
$b=0$ and $A(\tau)=\begin{pmatrix}0&0\\\ 0&D(\tau)\end{pmatrix}=B(\tau),$ say.
In that case, (17) becomes

$\eta\mid\omega,\tau,y\sim
N((M^{\top}\Omega(\omega)M+B(\tau))^{-1}M^{\top}\kappa,(M^{\top}\Omega(\omega)M+B(\tau))^{-1}).$

By using the method of calculating the inverse of a partitioned matrix, the
covariance matrix is

$\displaystyle(M^{\top}\Omega(\omega)M+B(\tau))^{-1}$
$\displaystyle=\begin{pmatrix}X^{\top}\Omega(\omega)X&X^{\top}\Omega(\omega)Z\\\
Z^{\top}\Omega(\omega)X&Z^{\top}\Omega(\omega)Z+D(\tau)\end{pmatrix}^{-1}$
$\displaystyle=\begin{pmatrix}(\tilde{X}^{\top}\tilde{X})^{-1}+\tilde{R}\tilde{S}^{-1}\tilde{R}^{\top}&-\tilde{R}\tilde{S}^{-1}\\\
-\tilde{S}^{-1}\tilde{R}^{\top}&\tilde{S}^{-1}\end{pmatrix},$

where $\tilde{X}=\Omega(\omega)^{\frac{1}{2}}X$,
$\tilde{Z}=\Omega(\omega)^{\frac{1}{2}}Z$,
$P_{\tilde{X}}=\tilde{X}(\tilde{X}^{\top}\tilde{X})^{-1}\tilde{X}^{\top}$,
$\tilde{S}=\tilde{Z}^{\top}(I-P_{\tilde{X}})\tilde{Z}+D(\tau)$, and
$\tilde{R}=(\tilde{X}^{\top}\tilde{X})^{-1}\tilde{X}^{\top}\tilde{Z}$. For the
mean vector, it follows that

$\displaystyle(M^{\top}\Omega(\omega)M+B(\tau))^{-1}M^{\top}\kappa$
$\displaystyle=\begin{pmatrix}(\tilde{X}^{\top}\tilde{X})^{-1}+\tilde{R}\tilde{S}^{-1}\tilde{R}^{\top}&-\tilde{R}\tilde{S}^{-1}\\\
-\tilde{S}^{-1}\tilde{R}^{\top}&\tilde{S}^{-1}\end{pmatrix}\begin{pmatrix}X^{\top}\kappa\\\
Z^{\top}\kappa\end{pmatrix}$
$\displaystyle=\begin{pmatrix}(\tilde{X}^{\top}\tilde{X})^{-1}X^{\top}\kappa+\tilde{R}\tilde{S}^{-1}\tilde{R}^{\top}X^{\top}\kappa-\tilde{R}\tilde{S}^{-1}Z^{\top}\kappa\\\
-\tilde{S}^{-1}\tilde{R}^{\top}X^{\top}\kappa+\tilde{S}^{-1}Z^{\top}\kappa\end{pmatrix}.$
(22)

The first element in the right-hand side of (A) is the mean vector for
$\beta$, while the second element in it is the mean vector for $u$. Thus,

$\displaystyle u\mid\omega,\tau,y\sim
N(-\tilde{S}^{-1}\tilde{R}^{\top}X^{\top}\kappa+\tilde{S}^{-1}Z^{\top}\kappa,\tilde{S}^{-1}).$
(23)

######  Lemma 1.

Let $R_{j}$ be a $q_{j}\times q$ matrix consisting of $0$’s and $1$’s such
that $R_{j}u=u_{j}$. Then

$\displaystyle(R_{j}\tilde{S}^{-1}{R_{j}}^{\top})^{-1}\preceq\bigg{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)+\tau_{j}\bigg{)}{\mathrm{I}}_{q_{j}}.$

Here, for two matrices A and B, $A\preceq B$ means $B-A$ is a positive
semidefinite matrix.

###### Proof.

Let $\lambda_{max}$ denote the largest eigenvalue for
$\tilde{Z}^{\top}({\mathrm{I}}-P_{\tilde{X}})\tilde{Z}$, then

$\displaystyle\tilde{S}=\tilde{Z}^{\top}({\mathrm{I}}-P_{\tilde{X}})\tilde{Z}+D(\tau)\preceq\lambda_{max}{\mathrm{I}}_{q}+D(\tau)\preceq{\mathrm{tr}}(\tilde{Z}^{\top}({\mathrm{I}}-P_{\tilde{X}})\tilde{Z}){\mathrm{I}}_{q}+D(\tau),$

where $D(\tau)=\oplus_{j=1}^{r}\tau_{j}{\mathrm{I}}_{q_{j}}$ as defined
before, and the second inequality follows from the fact that
$\tilde{Z}^{\top}({\mathrm{I}}-P_{\tilde{X}})\tilde{Z}$ is a positive
semidefinite matrix. Now

$\displaystyle{\mathrm{tr}}(\tilde{Z}^{\top}({\mathrm{I}}-P_{\tilde{X}})\tilde{Z})\leq{\mathrm{tr}}(\tilde{Z}^{\top}\tilde{Z})$
$\displaystyle={\mathrm{tr}}(Z^{\top}\Omega(\omega)Z)$
$\displaystyle={\mathrm{tr}}\bigg{(}\sum_{i=1}^{n}\omega_{i}z_{i}z_{i}^{\top}\bigg{)}$
$\displaystyle=\sum_{i=1}^{n}{\mathrm{tr}}(\omega_{i}z_{i}z_{i}^{\top})=\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(z_{i}z_{i}^{\top})\leq\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z),$

where $z_{i}^{\top}$ denotes the $i^{th}$ row of the $Z$ matrix, and the first
inequality is due to the fact that $\tilde{Z}^{\top}P_{\tilde{X}}\tilde{Z}$ is
a positive semidefinite matrix. Thus
$\tilde{S}\preceq\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z){\mathrm{I}}_{q}+D(\tau)$.
Hence,
$\tilde{S}^{-1}\succeq(\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z){\mathrm{I}}_{q}+D(\tau))^{-1}$.
Recall that $R_{j}u=u_{j}$. Extracting the result of the $j^{th}$ random
effect, we obtain:

$\displaystyle R_{j}\tilde{S}^{-1}{R_{j}}^{\top}\succeq
R_{j}\bigg{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z){\mathrm{I}}_{q}+D(\tau)\bigg{)}^{-1}{R_{j}}^{\top}=\bigg{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)+\tau_{j}\bigg{)}^{-1}{\mathrm{I}}_{q_{j}}.$

Thus, we have
$(R_{j}\tilde{S}^{-1}{R_{j}}^{\top})^{-1}\preceq(\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)+\tau_{j}){\mathrm{I}}_{q_{j}}.$
∎

##  Appendix B Some properties of the Pólya-Gamma distributions

######  Lemma 2.

Suppose $\omega\sim PG(a,b)$.

  1. 1.

If $a\geq 1,\,b\geq 0$, then for $0<s\leq 1$, ${\mathrm{E}}(\omega^{-s})\leq
2^{s}b^{s}+L(s)$, where $L(s)$ is a constant depending on $s$.

  2. 2.

If $a<1,\,b\geq 0$, then for $0<s<a$, ${\mathrm{E}}(\omega^{-s})\leq
2^{-s}(\pi^{2}+b^{2})^{s}\frac{\Gamma(a-s)}{\Gamma(a)}$.

###### Proof.

We first prove part $1$ for $a=1$. The probability density function of a
$PG(1,b)$ random variable is

$\displaystyle f(x\mid
1,b)=\cosh(b/2)\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi
x^{3}}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}-\frac{b^{2}}{2}x\Big{]},\,x>0.$

We consider the two cases $b=0$ and $b>0$ separately.
Case 1: $b=0$. Since $0<s\leq 1$, for any $x>0$, we have $x^{-s}\leq
x^{-1}+1$. Then,

$\displaystyle{\mathrm{E}}(\omega^{-s})\leq\int_{0}^{\infty}(x^{-1}+1)f(x\mid
1,0)dx=\int_{0}^{\infty}x^{-1}f(x\mid 1,0)dx+1.$

Now,

$\displaystyle\int_{0}^{\infty}x^{-1}f(x\mid 1,0)dx$
$\displaystyle=\int_{0}^{\infty}x^{-1}\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi
x^{3}}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}\Big{]}dx$
$\displaystyle=\int_{0}^{\infty}\sum_{\ell=0}^{\infty}(-1)^{\ell}x^{-\frac{5}{2}}\frac{(2\ell+1)}{\sqrt{2\pi}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}\Big{]}dx.$
(24)

Let
$h_{1}(x,\ell)=(-1)^{\ell}x^{-\frac{5}{2}}\frac{(2\ell+1)}{\sqrt{2\pi}}\exp\big{[}-\frac{(2\ell+1)^{2}}{8x}\big{]}$,
then

$\displaystyle\sum_{\ell=0}^{\infty}\int_{0}^{\infty}\mathinner{\\!\left\lvert
h_{1}(x,\ell)\right\rvert}dx$
$\displaystyle=\sum_{\ell=0}^{\infty}\frac{(2\ell+1)}{\sqrt{2\pi}}\int_{0}^{\infty}x^{-\frac{5}{2}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}\Big{]}dx$
$\displaystyle=8\sum_{\ell=0}^{\infty}\frac{1}{(2\ell+1)^{2}}<\infty.$

Hence, $\mathinner{\\!\left\lvert h_{1}\right\rvert}$ is integrable with
respect to the product measure of the counting measure and the Lebesgue
measure. By Fubini’s Theorem, from (B) we have

$\displaystyle\int_{0}^{\infty}x^{-1}f(x\mid 1,0)dx$
$\displaystyle=\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}\int_{0}^{\infty}x^{-\frac{5}{2}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}\Big{]}dx$
$\displaystyle=8\sum_{\ell=0}^{\infty}(-1)^{\ell}(2\ell+1)^{-2}=8C,$ (25)

where C is Catalan’s constant. Hence, ${\mathrm{E}}(\omega^{-s})\leq 8C+1$.
Case 2: $b>0$. Note that

$\displaystyle{\mathrm{E}}(\omega^{-s})$
$\displaystyle=\int_{0}^{\infty}x^{-s}f(x\mid 1,b)dx$
$\displaystyle=\int_{0}^{\infty}x^{-s-\frac{3}{2}}\cosh(b/2)\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}-\frac{b^{2}}{2}x\Big{]}dx.$
(26)

According to $10.32.10$ in Olver et al. (2010), we have

$\displaystyle\int_{0}^{\infty}x^{-s-\frac{3}{2}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}-\frac{b^{2}}{2}x\Big{]}dx=2K_{s+\frac{1}{2}}\Big{(}\frac{b(2\ell+1)}{2}\Big{)}\Big{(}\frac{2b}{2\ell+1}\Big{)}^{s+\frac{1}{2}},$
(27)

where $K_{v}(\cdot)$ is the modified Bessel function of the second kind of
order $v$. For $x>0$, according to $10.32.8$ in Olver et al. (2010),

$\displaystyle K_{s+\frac{1}{2}}(x)$
$\displaystyle=\frac{\sqrt{\pi}(\frac{1}{2}x)^{s+\frac{1}{2}}}{\Gamma(s+1)}\int_{1}^{\infty}e^{-xt}(t^{2}-1)^{s}dt$
$\displaystyle=\frac{\sqrt{\pi}(\frac{1}{2}x)^{s+\frac{1}{2}}}{\Gamma(s+1)}e^{-x}\int_{0}^{\infty}e^{-xt}(t^{2}+2t)^{s}dt$
(28)
$\displaystyle\leq\frac{\sqrt{\pi}(\frac{1}{2}x)^{s+\frac{1}{2}}}{\Gamma(s+1)}e^{-x}\int_{0}^{\infty}e^{-xt}(t^{2s}+2^{s}t^{s})dt$
$\displaystyle=\frac{\sqrt{\pi}(\frac{1}{2}x)^{s+\frac{1}{2}}}{\Gamma(s+1)}e^{-x}\Big{(}\frac{\Gamma(2s+1)}{x^{2s+1}}+2^{s}\frac{\Gamma(s+1)}{x^{s+1}}\Big{)}$
$\displaystyle=\sqrt{\pi}e^{-x}\Big{[}\frac{\Gamma(2s+1)}{\Gamma(s+1)}2^{-s-1/2}x^{-s-1/2}+2^{-1/2}x^{-1/2}\Big{]}.$
(29)

Also, from (28) we have

$\displaystyle
K_{s+\frac{1}{2}}(x)\geq\frac{\sqrt{\pi}(\frac{1}{2}x)^{s+\frac{1}{2}}}{\Gamma(s+1)}e^{-x}\int_{0}^{\infty}e^{-xt}2^{s}t^{s}dt=\sqrt{\pi}e^{-x}2^{-1/2}x^{-1/2}.$
(30)

Let
$h_{2}(x,\ell)=x^{-s-3/2}\cosh(b/2)(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}\exp\big{[}-\frac{(2\ell+1)^{2}}{8x}-\frac{b^{2}}{2}x\big{]}$,
then

$\displaystyle\sum_{\ell=0}^{\infty}\int_{0}^{\infty}\mathinner{\\!\left\lvert
h_{2}(x,\ell)\right\rvert}dx$
$\displaystyle=\sum_{\ell=0}^{\infty}\cosh(b/2)\frac{(2\ell+1)}{\sqrt{2\pi}}\int_{0}^{\infty}x^{-s-\frac{3}{2}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}-\frac{b^{2}}{2}x\Big{]}dx$
$\displaystyle=\cosh(b/2)\sum_{\ell=0}^{\infty}\frac{(2\ell+1)}{\sqrt{2\pi}}2K_{s+\frac{1}{2}}\Big{(}\frac{b(2\ell+1)}{2}\Big{)}\big{(}\frac{2b}{2\ell+1}\big{)}^{s+\frac{1}{2}}$
$\displaystyle\leq
2\cosh(b/2)\sum_{\ell=0}^{\infty}\frac{(2\ell+1)}{\sqrt{2\pi}}\sqrt{\pi}e^{-\frac{b(2\ell+1)}{2}}\Big{[}\frac{\Gamma(2s+1)}{\Gamma(s+1)}2^{-s-1/2}$
$\displaystyle\hskip
18.06749pt\times\big{(}\frac{b(2\ell+1)}{2}\big{)}^{-s-1/2}+2^{-1/2}\big{(}\frac{b(2\ell+1)}{2}\big{)}^{-1/2}\Big{]}\Big{(}\frac{2b}{2\ell+1}\Big{)}^{s+\frac{1}{2}}$
$\displaystyle=2^{s}(1+e^{-b})\Bigg{[}\sum_{\ell=0}^{\infty}\frac{e^{-b\ell}}{(2\ell+1)^{2s}}\frac{\Gamma(2s+1)}{\Gamma(s+1)}+\sum_{\ell=0}^{\infty}\frac{e^{-b\ell}}{(2\ell+1)^{s}}b^{s}\Bigg{]}<\infty.$

The second equality follows by using (27). The inequality is based on (29).
The convergence of the two series in the last step can be obtained by
utilizing the ratio test. Hence, $\mathinner{\\!\left\lvert
h_{2}\right\rvert}$ is integrable with respect to the product measure of the
counting measure and the Lebesgue measure. By Fubini’s Theorem and (27), (B)
becomes

$\displaystyle{\mathrm{E}}(\omega^{-s})=\cosh(b/2)\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}2K_{s+\frac{1}{2}}\Big{(}\frac{b[2\ell+1]}{2}\Big{)}\Big{(}\frac{2b}{2\ell+1}\Big{)}^{s+\frac{1}{2}}.$
(31)

When $\ell$ is even, applying (29) to (31), and when $\ell$ is odd, applying
(30) to (31), we obtain

$\displaystyle{\mathrm{E}}(\omega^{-s})$ $\displaystyle\leq
2\cosh(b/2)\Bigg{\\{}\sum_{even\;\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}\sqrt{\pi}e^{-\frac{b(2\ell+1)}{2}}\Big{[}\frac{\Gamma(2s+1)}{\Gamma(s+1)}2^{-s-1/2}\Big{\\{}\frac{b(2\ell+1)}{2}\Big{\\}}^{-s-1/2}$
$\displaystyle\quad+(b(2\ell+1))^{-1/2}\Big{]}-\sum_{odd\;\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}\sqrt{\pi}e^{-\frac{b(2\ell+1)}{2}}2^{-1/2}\Big{(}\frac{b(2\ell+1)}{2}\Big{)}^{-1/2}\Bigg{\\}}$
$\displaystyle\quad\quad\times\bigg{(}\frac{2b}{2\ell+1}\bigg{)}^{s+\frac{1}{2}}$
$\displaystyle=e^{b/2}(1+e^{-b})\bigg{\\{}\sum_{even\;\ell}(2\ell+1)^{-2s}2^{s}e^{-b(\ell+1/2)}\frac{\Gamma(2s+1)}{\Gamma(s+1)}$
$\displaystyle\quad+\sum_{even\;\ell}(2\ell+1)^{-s}(2b)^{s}e^{-b(\ell+1/2)}-\sum_{odd\;\ell}(2\ell+1)^{-s}(2b)^{s}e^{-b(\ell+1/2)}\bigg{\\}}$
$\displaystyle=(1+e^{-b})b^{s}\sum_{\ell=0}^{\infty}(-e^{-b})^{\ell}(\ell+1/2)^{-s}+(1+e^{-b})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}$
$\displaystyle\quad\quad\times\sum_{even\;\ell}e^{-b\ell}(\ell+1/2)^{-2s}$
$\displaystyle=(1+e^{-b})b^{s}\Phi(-e^{-b},s,1/2)+(1+e^{-b})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}\sum_{k=0}^{\infty}e^{-2bk}(2k+1/2)^{-2s}$
$\displaystyle=(1+e^{-b})\frac{b^{s}}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-t/2}}{1+e^{-b-t}}dt+(1+e^{-b})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}$
$\displaystyle\quad\quad\times\sum_{k=0}^{\infty}e^{-2bk}(2k+1/2)^{-2s},$ (32)

where $\Phi(\cdot)$ is the Lerch transcendent function.
For fixed $s>0$, let

$\displaystyle f(b)$
$\displaystyle=(1+e^{-b})\frac{b^{s}}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-t/2}}{1+e^{-b-t}}dt-2^{s}b^{s}$
$\displaystyle=\frac{b^{s}}{\Gamma(s)}\Big{[}(1+e^{-b})\int_{0}^{\infty}\frac{t^{s-1}e^{-t/2}}{1+e^{-b-t}}dt-\int_{0}^{\infty}t^{s-1}e^{-t/2}dt\Big{]}$
$\displaystyle=\frac{b^{s}e^{-b}}{\Gamma(s)}\int_{0}^{\infty}\frac{(1-e^{-t})}{1+e^{-b-t}}t^{s-1}e^{-t/2}dt.$

Since $(1-e^{-t})t^{s-1}e^{-t/2}/(1+e^{-b-t})\leq t^{s-1}e^{-t/2}$, which is
integrable, by the Dominated Convergence Theorem (DCT), it follows that $f(b)$
is a continuous function of $b$. Another application of DCT shows that

$\displaystyle\lim_{b\to\infty}\int_{0}^{\infty}\frac{1-e^{-t}}{1+e^{-b-t}}t^{s-1}e^{-t/2}dt=\int_{0}^{\infty}(1-e^{-t})t^{s-1}e^{-t/2}dt\leq\int_{0}^{\infty}t^{s-1}e^{-t/2}dt=2^{s}\Gamma(s).$

Hence, $\lim_{b\to\infty}f(b)=0$. Since $f(b)$ is a continuous function of
$b$, $f(0)=0$ and $\lim_{b\to\infty}f(b)=0$, we can conclude that
$\mathinner{\\!\left\lvert f(b)\right\rvert}$ can be bounded by a positive
constant $f_{0}$, hence,

$\displaystyle(1+e^{-b})\frac{b^{s}}{\Gamma(s)}\int_{0}^{\infty}\frac{t^{s-1}e^{-t/2}}{1+e^{-b-t}}dt\leq
2^{s}b^{s}+f_{0}.$ (33)

As for the second term in (B), we have

$\displaystyle(1+e^{-b})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}\sum_{k=0}^{\infty}e^{-2bk}(2k+1/2)^{-2s}$
$\displaystyle\leq(1+e^{-b})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}[1/(e^{2b}-1)+4^{s}].$
(34)

Here, the inequality is due to the fact $(2k+1/2)^{-2s}\leq 1$ for $k\geq 1$.
Note that for $b\geq\epsilon$, where $\epsilon>0$ is arbitrary, the upper
bound of (B) becomes

$\displaystyle(1+e^{-b})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}[1/(e^{2b}-1)+4^{s}]\leq(1+e^{-\epsilon})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}[1/(e^{2\epsilon}-1)+4^{s}].$

Thus, combining (33) with the above result, from (B) we have for
$b\geq\epsilon$,

$\displaystyle{\mathrm{E}}(\omega^{-s})\leq 2^{s}b^{s}+f_{0}+L(s,\epsilon),$
(35)

where
$L(s,\epsilon)=(1+e^{-\epsilon})2^{-s}\frac{\Gamma(2s+1)}{\Gamma(s+1)}[1/(e^{2\epsilon}-1)+4^{s}]$.

Now, we consider $0<b<\epsilon$. Let $k(b)\equiv{\mathrm{E}}(\omega^{-s})$,
where $\omega\sim PG(1,b)$. Then,

$\displaystyle\lim_{b\to 0}k(b)=\lim_{b\to 0}\cosh(b/2)\lim_{b\to
0}\int_{0}^{\infty}j(b,x)dx=\lim_{b\to 0}\int_{0}^{\infty}j(b,x)dx,$ (36)

where

$\displaystyle
j(b,x)=x^{-s-\frac{3}{2}}\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}-\frac{b^{2}}{2}x\Big{]}.$

Note that

$\displaystyle
j(b,x)\leq(x^{-1}+1)x^{-\frac{3}{2}}\sum_{\ell=0}^{\infty}(-1)^{\ell}\frac{(2\ell+1)}{\sqrt{2\pi}}\exp\Big{[}-\frac{(2\ell+1)^{2}}{8x}\Big{]}=j(x),\,\text{say}.$

From (B), it follows that $\int_{0}^{\infty}j(x)dx\leq 8C+1$. Then by the DCT,
from (36), we have

$\displaystyle\lim_{b\to 0}k(b)=\lim_{b\to 0}\int_{0}^{\infty}j(b,x)dx=k(0).$

So, $k(b)={\mathrm{E}}(\omega^{-s})$ is continuous at $b=0$. Recall that
${\mathrm{E}}(\omega^{-s})\leq 8C+1$ for $b=0$ and $0<s\leq 1$. Thus, there
exists some small $\epsilon>0$ such that ${\mathrm{E}}(\omega^{-s})\leq 8C+2$
for $0<b<\epsilon$. Combining this result with (35), we have
${\mathrm{E}}(\omega^{-s})\leq 2^{s}b^{s}+L(s)$, where
$L(s)=\max\\{f_{0}+L(s,\epsilon),8C+2\\}$. Thus, part $1$ is proved for $a=1$.

Next, we prove the conclusion for $a>1$. From Polson et al. (2013), when
$\omega\sim PG(a,b)$, we have

$\displaystyle\omega\overset{d}{=}\frac{1}{2\pi^{2}}\sum_{\ell=1}^{\infty}\frac{g_{\ell}}{(\ell-1/2)^{2}+b^{2}/(4\pi^{2})},$
(37)

where $g_{\ell}$’s are mutually independent Gamma$(a,1)$ random variables.
Since $a>1$, $g_{\ell}\overset{d}{=}\tilde{g}_{\ell}+g^{*}_{\ell}$, where
$\tilde{g}_{\ell}$ and $g^{*}_{\ell}$ are independent random variables
following Gamma $(a-1,1)$ and Gamma $(1,1)$, respectively. Let
$x_{1}=(1/[2\pi^{2}])\sum_{\ell=1}^{\infty}g^{*}_{\ell}/[(\ell-1/2)^{2}+b^{2}/(4\pi^{2})]$.
Then, $x_{1}\sim PG(1,b)$. Thus, we have ${\mathrm{E}}(x_{1}^{-s})\leq
2^{s}b^{s}+L(s)$. Since for $0<s\leq 1$,
${\mathrm{E}}(\omega^{-s})\leq{\mathrm{E}}(x_{1}^{-s})$, the same conclusion
follows for $\omega\sim PG(a,b)$, where $a>1$. Thus, the proof for part $1$ is
complete.

Next, we prove part $2$. From (37), we have

$\displaystyle{\mathrm{E}}\,\omega^{-s}$
$\displaystyle={\mathrm{E}}\,\Big{[}\frac{1}{2\pi^{2}}\sum_{\ell=1}^{\infty}\frac{g_{\ell}}{(\ell-1/2)^{2}+b^{2}/(4\pi^{2})}\Big{]}^{-s}$
$\displaystyle\leq{\mathrm{E}}\,\Big{[}\frac{1}{2\pi^{2}}\frac{g_{1}}{(1-1/2)^{2}+b^{2}/(4\pi^{2})}\Big{]}^{-s}$
$\displaystyle=\Big{(}\frac{\pi^{2}+b^{2}}{2}\Big{)}^{s}\int_{0}^{\infty}g_{1}^{-s}\frac{1}{\Gamma(a)}g_{1}^{a-1}\exp{(-g_{1})}\,d\,g_{1}$
$\displaystyle=2^{-s}(\pi^{2}+b^{2})^{s}\frac{\Gamma(a-s)}{\Gamma(a)}.$

Thus, the proof for part $2$ is complete. ∎

######  Lemma 3.

If $\omega\sim PG(a,b),\,a>0,\,b\geq 0$, then ${\mathrm{E}}\,\omega\leq a/4$.

###### Proof.

From (37), we have

$\displaystyle{\mathrm{E}}\,\omega$
$\displaystyle={\mathrm{E}}\Big{[}\frac{1}{2\pi^{2}}\sum_{\ell=1}^{\infty}\frac{g_{\ell}}{(\ell-1/2)^{2}+b^{2}/(4\pi^{2})}\Big{]}$
$\displaystyle\leq{\mathrm{E}}\Big{[}\frac{1}{2\pi^{2}}\sum_{\ell=1}^{\infty}\frac{g_{\ell}}{(\ell-1/2)^{2}}\Big{]}$
$\displaystyle=\frac{1}{2\pi^{2}}\sum_{\ell=1}^{\infty}\frac{a}{(\ell-1/2)^{2}}=a{\mathrm{E}}\,\omega_{2},$
(38)

where $\omega_{2}\sim PG(1,0)$. By Polson et al. (2013),
$PG(1,0)=J^{*}(1,0)/4$. From Devroye (2009), the density for $J^{*}(1,0)$ is

$\displaystyle
f^{*}(x)=\pi\sum_{\ell=0}^{\infty}(-1)^{\ell}(\ell+1/2)\exp\Big{[}-\frac{(\ell+1/2)^{2}\pi^{2}x}{2}\Big{]},$

then,

$\displaystyle{\mathrm{E}}\,J^{*}(1,0)$
$\displaystyle=\int_{0}^{\infty}x\,\pi\sum_{\ell=0}^{\infty}(-1)^{\ell}(\ell+1/2)\exp\Big{[}-\frac{(\ell+1/2)^{2}\pi^{2}x}{2}\Big{]}dx$
$\displaystyle=\sum_{\ell=0}^{\infty}\int_{0}^{\infty}x\,\pi(-1)^{\ell}(\ell+1/2)\exp\Big{[}-\frac{(\ell+1/2)^{2}\pi^{2}x}{2}\Big{]}dx$
$\displaystyle=\sum_{\ell=0}^{\infty}\pi(-1)^{\ell}(\ell+1/2)\Gamma(2)\Big{(}\frac{1}{2}\Big{(}\ell+\frac{1}{2}\Big{)}^{2}\pi^{2}\Big{)}^{-2}$
$\displaystyle=\frac{4}{\pi^{3}}\sum_{\ell=0}^{\infty}\frac{(-1)^{\ell}}{(\ell+1/2)^{3}}=\frac{4}{\pi^{3}}\Phi\Big{(}-1,3,\frac{1}{2}\Big{)}=\frac{4}{\pi^{3}}\frac{\pi^{3}}{4}=1,$

where $\Phi(\cdot)$ is the Lerch transcendent function and
$\Phi(-1,3,\frac{1}{2})=\pi^{3}/4$. Note that, we achieve the second equality
above by applying the following logic. Let
$h_{3}(x,\ell)=x\,\pi(-1)^{\ell}(\ell+\frac{1}{2})\exp\big{[}-\frac{(\ell+1/2)^{2}\pi^{2}x}{2}\big{]}$,
then

$\displaystyle\sum_{\ell=0}^{\infty}\int_{0}^{\infty}\mathinner{\\!\left\lvert
h_{3}(x,\ell)\right\rvert}dx$
$\displaystyle=\sum_{\ell=0}^{\infty}\int_{0}^{\infty}x\,\pi\Big{(}\ell+\frac{1}{2}\Big{)}\exp\Big{[}-\frac{(\ell+1/2)^{2}\pi^{2}x}{2}\Big{]}dx$
$\displaystyle=\frac{4}{\pi^{3}}\sum_{\ell=0}^{\infty}\frac{1}{(\ell+1/2)^{3}}<\infty,$

Hence, $h_{3}(x,\ell)$ is integrable with respect to the product measure of
the counting measure and the Lebesgue measure. Thus, by Fubini’s Theorem, the
second equality follows. Consequently,
${\mathrm{E}}\omega_{2}={\mathrm{E}}\,J^{*}(1,0)/4=1/4$. From (B), it follows
${\mathrm{E}}\,\omega\leq a/4$. ∎

######  Remark 3.

Wang and Roy (2018c) proved Lemma 2 in the special case when $a=1$. Although
their result is correct as stated, their proof has an error which can be
repaired following the techniques used in the proof of Lemma 2 here. Lemma 3
for $a=1$ is also proved in Wang and Roy (2018c).

##  Appendix C Proof of Theorem 1

###### Proof.

We first prove the geometric ergodicity of the $\tilde{\Psi}$ chain by
establishing a drift condition. We consider the following drift function

$\displaystyle V(\eta)=\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta\right\rvert}+\sum_{j=1}^{r}(u_{j}^{\top}u_{j})^{-c},$ (39)

where $c\in(0,1/2)$ to be determined later.

Since $M$ has full rank, $V(\eta)\mathrel{\mathop{\mathchar
58\relax}}\mathbb{R}^{p+q}\backslash N\rightarrow[0,\infty)$ is unbounded off
compact sets. We prove that, for any
$\eta,\eta^{\prime}\in\mathbb{R}^{p+q}\backslash N$, there exist constants
$\rho\in[0,1)$ and $L>0$ such that

$\displaystyle{\mathrm{E}}[V(\eta)\mid\eta\prime]={\mathrm{E}}\\{{\mathrm{E}}[V(\eta)\mid\omega,\tau,y]\mid\eta\prime,y\\}\leq\rho
V(\eta\prime)+L.$ (40)

The first term in the drift function is
$\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta\right\rvert}=l^{\top}M\eta$, where
$l=(l_{1},l_{2},...,l_{n})$ is defined as $l_{i}=1$ if $m_{i}^{\top}\eta\geq
0,l_{i}=-1$ if $m_{i}^{\top}\eta<0$. Since $Q=0$, from (17) we have

$\displaystyle{\mathrm{E}}\biggl{[}\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta\right\rvert}\bigg{|}\omega,\tau,y\biggr{]}$
$\displaystyle=l^{\top}M(M^{\top}\Omega(\omega)M+B(\tau))^{-1}M^{\top}\kappa$
$\displaystyle\leq\sqrt{l^{\top}M(M^{\top}\Omega(\omega)M+B(\tau))^{-1}M^{\top}l}\sqrt{\kappa^{\top}M(M^{\top}\Omega(\omega)M+B(\tau))^{-1}M^{\top}\kappa}$
$\displaystyle\leq\sqrt{l^{\top}M(M^{\top}\Omega(\omega)M)^{-1}M^{\top}l}\sqrt{\kappa^{\top}M(M^{\top}\Omega(\omega)M)^{-1}M^{\top}\kappa}$
$\displaystyle=\sqrt{l^{\top}\Omega(\omega)^{-1/2}P_{\Omega(\omega)^{1/2}M}\Omega(\omega)^{-1/2}l}\sqrt{\kappa^{\top}M(M^{\top}\Omega(\omega)M)^{-1}M^{\top}\kappa}$
$\displaystyle\leq\sqrt{\sum_{i=1}^{n}\frac{1}{\omega_{i}}}\sqrt{\kappa^{\top}M(M^{\top}\Omega(\omega)M)^{-1}M^{\top}\kappa},$
(41)

where the first inequality follows from the Cauchy-Schwarz inequality,
$P_{\Omega(\omega)^{1/2}M}\equiv\Omega(\omega)^{1/2}\linebreak
M(M^{\top}\Omega(\omega)M)^{-1}M^{\top}\Omega(\omega)^{1/2}$ is a projection
matrix, and the third inequality follows from the fact that
${\mathrm{I}}\succeq P_{\Omega(\omega)^{1/2}M}$. Recall that
$k_{i}=y_{i}-1/2,\,i=1,...,n$. Define, $v_{i}=-2k_{i}m_{i}$ as the $i^{th}$
row of an $n\times(p+q)$ matrix $V$. Note that
$v_{i}v_{i}^{\top}=m_{i}m_{i}^{\top},\,i=1,...,n$. Since the conditions $3$
and $4$ in Theorem 1 are in force, by Lemma $3$ in Wang and Roy (2018c), for
the second part of (C), we have

$\displaystyle\sqrt{\kappa^{\top}M(M^{\top}\Omega(\omega)M)^{-1}M^{\top}\kappa}=\sqrt{\frac{1}{4}1^{\top}V(V^{\top}\Omega(\omega)V)^{-1}V^{\top}1}\leq\sqrt{\frac{\rho_{1}}{4}\sum_{i=1}^{n}\frac{1}{\omega_{i}}},$
(42)

where $\rho_{1}\in[0,1)$ is a constant. Applying (42) to (C), we have

$\displaystyle{\mathrm{E}}\bigg{[}\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta\right\rvert}\bigg{|}\omega,\tau,y\bigg{]}\leq\frac{\sqrt{\rho_{1}}}{2}\sum_{i=1}^{n}\frac{1}{\omega_{i}}.$
(43)

Next, we consider the inner expectation of the second term in the drift
function (39). Note that, for $c\in(0,1/2)$, we have

$\displaystyle{\mathrm{E}}\Big{[}(u_{j}^{\top}u_{j})^{-c}\Big{|}\omega,\tau,y\Big{]}$
$\displaystyle=\bigg{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)+\tau_{j}\bigg{)}^{c}{\mathrm{E}}\bigg{[}\bigg{\\{}u_{j}^{\top}\bigg{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)+\tau_{j}\bigg{)}{\mathrm{I}}_{q_{j}}u_{j}\bigg{\\}}^{-c}\bigg{|}\omega,\tau,y\bigg{]}$
$\displaystyle\leq\bigg{[}\bigg{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)\bigg{)}^{c}+\tau_{j}^{c}\bigg{]}{\mathrm{E}}\bigg{[}\bigg{\\{}u_{j}^{\top}(R_{j}\tilde{S}^{-1}{R_{j}}^{\top})^{-1}u_{j}\bigg{\\}}^{-c}\bigg{|}\omega,\tau,y\bigg{]}$
$\displaystyle\leq\frac{2^{-c}\Gamma(-c+q_{j}/2)}{\Gamma(q_{j}/2)}\bigg{[}\bigg{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)\bigg{)}^{c}+\tau_{j}^{c}\bigg{]},$
(44)

where the first inequality follows from Lemma 1 and the fact that
$(a+b)^{s}\leq a^{s}+b^{s}$ for $a>0,\,b>0,$ and $0\leq s<1$. For the last
inequality, note that, by (23), we have $u_{j}\mid\omega,\tau,y\sim
N(R_{j}(-\tilde{S}^{-1}\tilde{R}^{\top}X^{\top}\kappa+\tilde{S}^{-1}Z^{\top}\kappa),R_{j}\tilde{S}^{-1}{R_{j}}^{\top})$
where $R_{j}$ is a $q_{j}\times q$ matrix consisting of $0$’s and $1$’s such
that $R_{j}u=u_{j}$. Thus, given $\omega,\tau,y$,
$(R_{j}\tilde{S}^{-1}{R_{j}}^{\top})^{-\frac{1}{2}}u_{j}$ has a multivariate
normal distribution with the identity covariance matrix. Hence, conditional on
$\omega,\tau,y$, the distribution of
$u_{j}^{\top}(R_{j}\tilde{S}^{-1}{R_{j}}^{\top})^{-1}u_{j}$ is
$\chi_{q_{j}}^{2}(w)$, where $w$ is some (unimportant) noncentrality parameter
and $q_{j}$ is the degrees of freedom for the relevant Chi-square
distribution. Therefore, by Lemma 4 in Román and Hobert (2012), we have

$\displaystyle{\mathrm{E}}\Big{[}\Big{\\{}u_{j}^{\top}(R_{j}\tilde{S}^{-1}{R_{j}}^{\top})^{-1}u_{j}\Big{\\}}^{-c}\Big{|}\omega,\tau,y\Big{]}\leq\frac{2^{-c}\Gamma(-c+q_{j}/2)}{\Gamma(q_{j}/2)}.$

Applying the above result, the inequality in (C) is obtained.

Combining (43) and (C), from (39), we have

$\displaystyle
E[V(\eta)\mid\omega,\tau,y]\leq\frac{\sqrt{\rho_{1}}}{2}\sum_{i=1}^{n}\frac{1}{\omega_{i}}+\sum_{j=1}^{r}\frac{2^{-c}\Gamma(-c+q_{j}/2)}{\Gamma(q_{j}/2)}\Big{[}\Big{(}\sum_{i=1}^{n}\omega_{i}{\mathrm{tr}}(Z^{\top}Z)\Big{)}^{c}+\tau_{j}^{c}\Big{]}.$
(45)

Next, we consider the outer expectation in (40). By Lemma 2, we have

$\displaystyle{\mathrm{E}}\Big{[}\frac{\sqrt{\rho_{1}}}{2}\sum_{i=1}^{n}\frac{1}{\omega_{i}}\mid\eta^{{}^{\prime}},y\Big{]}$
$\displaystyle\leq\Big{[}2\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta^{\prime}\right\rvert}+nL(1)\Big{]}\frac{\sqrt{\rho_{1}}}{2}$
$\displaystyle=\sqrt{\rho_{1}}\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta^{\prime}\right\rvert}+\frac{nL(1)\sqrt{\rho_{1}}}{2}.$ (46)

For the outer expectation of the other terms on the right hand side of (45),
we now consider the expectation for $\tau_{j}^{c}$. Recall from section 2.2
that $\tau_{j}\mid\eta^{\prime},y\overset{ind}{\sim}$
Gamma$(a_{j}+q_{j}/2,b_{j}+u_{j}^{\prime\top}u_{j}^{\prime}/2)$, $j=1,...,r$.
Then, it follows that

$\displaystyle{\mathrm{E}}[\tau_{j}^{c}\mid\eta^{\prime},y]=\frac{\Gamma(a_{j}+\frac{q_{j}}{2}+c)}{\Gamma(a_{j}+\frac{q_{j}}{2})}\Big{(}b_{j}+\frac{1}{2}u_{j}^{\prime\top}u_{j}^{\prime}\Big{)}^{-c}.$

Define, $G_{j}(-c)=2^{c}\Gamma(a_{j}+q_{j}/2+c)/\Gamma(a_{j}+q_{j}/2)$,
$j=1,2,...,r$. Hence,

$\displaystyle{\mathrm{E}}[\tau_{j}^{c}\mid\eta^{\prime},y]$
$\displaystyle=2^{-c}G_{j}(-c)\Big{(}b_{j}+\frac{u_{j}^{\prime\top}u_{j}^{\prime}}{2}\Big{)}^{-c}$
$\displaystyle\leq
G_{j}(-c)[(2b_{j})^{-c}{\mathrm{I}}_{(0,\infty)}(b_{j})+(u_{j}^{\prime\top}u_{j}^{\prime})^{-c}{\mathrm{I}}_{\\{0\\}}(b_{j})].$
(47)

Also,

$\displaystyle\sum_{j=1}^{r}\frac{2^{-c}\Gamma(q_{j}/2-c)}{\Gamma(q_{j}/2)}G_{j}(-c)(u_{j}^{\prime\top}u_{j}^{\prime})^{-c}{\mathrm{I}}_{\\{0\\}}(b_{j})\leq\delta_{1}(c)\sum_{j=1}^{r}(u_{j}^{\prime\top}u_{j}^{\prime})^{-c}$
(48)

where $\delta_{1}(c)=2^{-c}\max_{j\in
A}\frac{\Gamma(q_{j}/2-c)}{\Gamma(q_{j}/2)}G_{j}(-c)\geq 0$. Recall that
$A=\\{j\in\\{1,2,...,r\\}\mathrel{\mathop{\mathchar 58\relax}}b_{j}=0\\}$.
From the condition $1$ of Theorem 1, we have $a_{j}<0$ when $b_{j}=0$ (i.e.
when $j\in A$). From the proof of Proposition $2$ in Román and Hobert (2012),
it follows that, there exists $c\in C\equiv(0,1/2)\cap(0,\tilde{a})$, where
$\tilde{a}=-\max_{j\in A}a_{j}$, such that $\delta_{1}(c)<1$.

Using (C), (C), and Jensen’s inequality, from (45), we obtain

$\displaystyle{\mathrm{E}}[V(\eta)\mid\eta\prime]$
$\displaystyle={\mathrm{E}}\\{{\mathrm{E}}[V(\eta\mid\omega,\tau,y)]\mid\eta\prime,y\\}$
(49) $\displaystyle\leq\sqrt{\rho_{1}}\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta^{\prime}\right\rvert}+\frac{nL(1)\sqrt{\rho_{1}}}{2}+\sum_{j=1}^{r}\frac{2^{-c}\Gamma(-c+q_{j}/2)}{\Gamma(q_{j}/2)}$
$\displaystyle\quad\times\big{\\{}\big{(}{\mathrm{tr}}(Z^{\top}Z)\sum_{i=1}^{n}{\mathrm{E}}(\omega_{i}\mid\eta^{\prime},y)\big{)}^{c}+G_{j}(-c)$
$\displaystyle\quad\quad\times\big{[}(2b_{j})^{-c}{\mathrm{I}}_{(0,\infty)}(b_{j})+(u_{j}^{\prime\top}u_{j}^{\prime})^{-c}{\mathrm{I}}_{\\{0\\}}(b_{j})\big{]}\big{\\}}$
$\displaystyle\leq\sqrt{\rho_{1}}\sum_{i=1}^{n}\mathinner{\\!\left\lvert
m_{i}^{\top}\eta^{\prime}\right\rvert}+\delta_{1}(c)\sum_{j=1}^{r}(u_{j}^{\prime\top}u_{j}^{\prime})^{-c}+L\leq\rho
V(\eta\prime)+L,$

where the second inequality is due to (48) and Lemma 3. Here,
$\rho=\max\\{\sqrt{\rho_{1}},\delta_{1}(c)\\}$, and

$\displaystyle L=\frac{nL(1)\sqrt{\rho_{1}}}{2}$
$\displaystyle+\sum_{j=1}^{r}\frac{2^{-c}\Gamma(-c+q_{j}/2)}{\Gamma(q_{j}/2)}\Big{(}\frac{{\mathrm{tr}}(Z^{\top}Z)n}{4}\Big{)}^{c}$
$\displaystyle+\sum_{j=1}^{r}\frac{2^{-c}\Gamma(-c+q_{j}/2)}{\Gamma(q_{j}/2)}G_{j}(-c)(2b_{j})^{-c}{\mathrm{I}}_{(0,\infty)}(b_{j}).$

Recall that $\rho_{1},\delta_{1}(c)\in[0,1)$, thus $\rho<1$. Consequently,
(40) holds. In addition, we can show that $\tilde{\Psi}$ chain is a Feller
Markov chain by the following steps. Let $K(\eta^{\prime},\cdot)$ denote the
Mtf corresponding to (21). To prove that the $\tilde{\Psi}$ chain is a Feller
Markov chain is to show that $K(\eta^{\prime},A)$ is a lower semi-continuous
function on $\mathbb{R}^{p+q}\backslash N$ for each fixed open set $A$ in
$\mathbb{R}^{p+q}\backslash N$. For a sequence $\\{\eta^{\prime}_{m}\\}$,
using (21), Fatou’s Lemma and independence of the conditional distribution of
$\omega$ and $\tau$ given $(\eta^{\prime},y)$, we have

$\displaystyle\liminf_{m\rightarrow\infty}K(\eta^{\prime}_{m},A)$
$\displaystyle=\liminf_{m\rightarrow\infty}\int_{A}\tilde{k}(\eta\mid\eta^{\prime}_{m})d\eta$
$\displaystyle=\liminf_{m\rightarrow\infty}\int_{A}\int_{\mathbb{R}_{+}^{r}}\int_{\mathbb{R}_{+}^{n}}\pi(\eta\mid\omega,\tau,y)\pi(\omega,\tau\mid\eta^{\prime}_{m},y)d\omega\,d\,\tau\,d\eta$
$\displaystyle\geq\int_{A}\int_{\mathbb{R}_{+}^{r}}\int_{\mathbb{R}_{+}^{n}}\pi(\eta\mid\omega,\tau,y)\liminf_{m\rightarrow\infty}\pi(\omega,\tau\mid\eta^{\prime}_{m},y)d\omega\,d\,\tau\,d\eta$
$\displaystyle=\int_{A}\int_{\mathbb{R}_{+}^{r}}\int_{\mathbb{R}_{+}^{n}}\pi(\eta\mid\omega,\tau,y)\liminf_{m\rightarrow\infty}[\pi(\omega\mid\eta^{\prime}_{m},y)\pi(\tau\mid\eta^{\prime}_{m},y)]d\omega\,d\,\tau\,d\eta.$

Considering the conditions $1$ and $2$ in Theorem 1, for any fixed
$(\eta,\omega,y)$, both $\pi(\omega\mid\eta^{\prime}_{m},y)$ and
$\pi(\tau\mid\eta^{\prime}_{m},y)$ are continuous functions on
$\mathbb{R}^{p+q}\backslash N$. Hence, if
$\eta^{\prime}_{m}\rightarrow\eta^{\prime}$, then
$\liminf_{m\rightarrow\infty}K(\eta^{\prime}_{m},A)\geq K(\eta^{\prime},A)$,
and we can conclude that the $\tilde{\Psi}$ chain is a Feller Markov chain.
Thus, by Lemma $15.2.8$ in Meyn and Tweedie (1993), GE of the $\tilde{\Psi}$
chain is proved.

Next, using the similar techniques as in Wang and Roy (2018b) and (Román,
2012, Lemma 12), we can establish that the GE of the original chain $\Psi$
follows from that of $\tilde{\Psi}$. We include a proof here for completeness.
Let
$\mathsf{X}\equiv\mathbb{R}^{p+q},\,\tilde{\mathsf{X}}\equiv\mathbb{R}^{p+q}\backslash
N$. Let $K$ and $\tilde{K}$ denote the Mtfs of $\Psi$ and $\tilde{\Psi}$
chains respectively. Also, since the Lebesgue measure of $N$ is zero,
$\tilde{K}(x,B)=K(x,B)$, for any $x\in\tilde{\mathsf{X}}$ and
$B\in\mathcal{B}_{\tilde{\mathsf{X}}}=\\{\tilde{\mathsf{X}}\cap
A\mathrel{\mathop{\mathchar 58\relax}}A\in\mathcal{B}_{\mathsf{X}}\\}$, where
$\mathcal{B}_{\mathsf{X}}$ denotes the Borel $\sigma$-algebra of
$\mathbb{R}^{p+q}$ and $\mathcal{B}_{\tilde{\mathsf{X}}}$ denotes the Borel
$\sigma$-algebra of $\mathbb{R}^{p+q}\backslash N$, respectively.

Let $\mu$ and $\tilde{\mu}$ be the Lebesgue measures on $\mathsf{X}$ and
$\tilde{\mathsf{X}}$, respectively. As the Mtds are strictly positive for the
two chains, the $\Psi$ chain is $\mu$-irreducible and the $\tilde{\Psi}$ chain
is $\tilde{\mu}$-irreducible. Both chains are aperiodic. Note that, $\mu$ and
$\tilde{\mu}$ are also the maximal irreducibility measures of $\Psi$ and
$\tilde{\Psi}$ chains, respectively. By Theorem $15.0.1$ in Meyn and Tweedie
(1993), since we have established the GE of the $\tilde{\Psi}$ chain, there
exists a $v$-petite set $C\in\mathcal{B}_{\tilde{\mathsf{X}}}$,
$\rho_{C}<1,\,M_{C}<\infty$ and $\tilde{K}^{\infty}(C)>0$ such that
$\tilde{\mu}(C)>0$ and

$\displaystyle\mathinner{\\!\left\lvert\tilde{K}^{m}(x,C)-\tilde{K}^{\infty}(C)\right\rvert}<M_{C}\rho_{C}^{m},$
(50)

for all $x\in C$. Also, it can be shown that

$\displaystyle K^{m}(x,B)=\tilde{K}^{m}(x,B\cap\tilde{\mathsf{X}}),$ (51)

for any $x\in\tilde{\mathsf{X}}$ and $B\in\mathcal{B}_{\mathsf{X}}$. Note
that, $K^{m}$ and $\tilde{K}^{m}$ indicate the corresponding $m$-step Mtfs.
Thus, for $x\in C$, $K^{m}(x,C)=\tilde{K}^{m}(x,C)$. Then, (50) becomes
$\mathinner{\\!\left\lvert
K^{m}(x,C)-\tilde{K}^{\infty}(C)\right\rvert}<M_{C}\rho_{C}^{m}.$ Since
$\mu(N)=0$, we have $\mu(C)=\tilde{\mu}(C)$. Recall that $\tilde{\mu}(C)>0$,
thus $\mu(C)>0$. Note that, $C$ is $v$-petite for the $\tilde{\Psi}$ chain,
then for all $x\in C$ and $B\in\mathcal{B}_{\tilde{\mathsf{X}}}$,

$\displaystyle\sum_{m=0}^{\infty}\tilde{K}^{m}(x,B)a(m)\geq v(B),$ (52)

where $v$ is a nontrivial measure on $\mathcal{B}_{\tilde{\mathsf{X}}}$ and
$a(m)$ is a mass function on $\\{0,1,2,...\\}$. It can be shown that a
nontrivial measure on $\mathcal{B}_{\mathsf{X}}$, which is

$\displaystyle v^{*}(\cdot)=v(\cdot\cap\tilde{\mathsf{X}}),$ (53)

is well defined. Then for any $x\in C$ and any $B\in\mathcal{B}_{\mathsf{X}}$,
using (51), (52) and (53), we have

$\sum_{m=0}^{\infty}K^{m}(x,B)a(m)=\sum_{m=0}^{\infty}\tilde{K}^{m}(x,B\cap\tilde{\mathsf{X}})a(m)\geq
v(B\cap\tilde{\mathsf{X}})=v^{*}(B).$

Hence, $C$ is also a petite set for the $\Psi$ chain. Applying Theorem
$15.0.1$ in Meyn and Tweedie (1993) again, GE of the $\Psi$ chain is proved.
Hence, we show that GE of $\tilde{\Psi}$ implies that of the original chain
$\Psi$.

∎

## References

  * Albert and Chib (1993) Albert, J. H. and Chib, S. (1993).  Bayesian analysis of binary and polychotomous response data.  Journal of the American statistical Association, 88 669–679.
  * Chakraborty and Khare (2017) Chakraborty, S. and Khare, K. (2017).  Convergence properties of Gibbs samplers for Bayesian probit regression with proper priors.  Electronic Journal of Statistics, 11 177–210.
  * Chen et al. (2002) Chen, M.-H., Shao, Q.-M. and Xu, D. (2002).  Necessary and sufficient conditions on the properiety of posterior distributions for generalized linear mixed models.  Sankhyā: The Indian Journal of Statistics, Series A 57–85.
  * Chib and Ramamurthy (2010) Chib, S. and Ramamurthy, S. (2010).  Tailored randomized block MCMC methods with application to DSGE models.  Journal of Econometrics, 155 19–38.
  * Choi and Hobert (2013) Choi, H. M. and Hobert, J. P. (2013).  The Polya-Gamma Gibbs sampler for Bayesian logistic regression is uniformly ergodic.  Electronic Journal of Statistics, 7 2054–2064.
  * Cortez and Silva (2008) Cortez, P. and Silva, A. M. G. (2008).  Using data mining to predict secondary school student performance.  In Proceedings of 5th Annual Future Business Technology Conference (A. Brito and J. Teixeira, eds.). EUROSIS-ETI, 5–12.
  * Devroye (2009) Devroye, L. (2009).  On exact simulation algorithms for some distributions related to Jacobi theta functions.  Statistics & Probability Letters, 79 2251–2259.
  * Evangelou and Roy (2019) Evangelou, E. and Roy, V. (2019).  Estimation and prediction for spatial generalized linear mixed models with parametric links via reparameterized importance sampling.  Spatial statistics, 29 289–315.
  * Frühwirth-Schnatter and Frühwirth (2010) Frühwirth-Schnatter, S. and Frühwirth, R. (2010).  Data augmentation and MCMC for binary and multinomial logit models.  In Statistical Modelling and Regression Structures. Springer, 111–132.
  * Holmes and Held (2006) Holmes, C. C. and Held, L. (2006).  Bayesian auxiliary variable models for binary and multinomial regression.  Bayesian Analysis, 1 145–168.
  * Jones and Hobert (2001) Jones, G. L. and Hobert, J. P. (2001).  Honest exploration of intractable probability distributions via Markov chain Monte Carlo.  Statistical Science, 16 312–334.
  * Liu et al. (1994) Liu, J. S., Wong, W. H. and Kong, A. (1994).  Covariance structure of the Gibbs sampler with applications to the comparisons of estimators and augmentation schemes.  Biometrika, 81 27–40.
  * Meyn and Tweedie (1993) Meyn, S. P. and Tweedie, R. L. (1993).  Markov chains and stochastic stability.  Springer.
  * Olver et al. (2010) Olver, F. W., Lozier, D. W., F., B. R. and Clark, C. W. (2010).  NIST Handbook of Mathematical Functions.  Cambridge University Press.
  * Polson et al. (2013) Polson, N. G., Scott, J. G. and Windle, J. (2013).  Bayesian inference for logistic models using Pólya-Gamma latent variables.  Journal of the American statistical Association, 108 1339–1349.
  * Roberts and Rosenthal (1997) Roberts, G. O. and Rosenthal, J. S. (1997).  Geometric ergodicity and hybrid Markov chains.  Electronic Communications in Probability, 2 13–25.
  * Roberts and Rosenthal (2001) Roberts, G. O. and Rosenthal, J. S. (2001).  Markov chains and de-initializing processes.  Scandinavian Journal of Statistics, 28 489–504.
  * Roberts and Rosenthal (2004) Roberts, G. O. and Rosenthal, J. S. (2004).  General state space Markov chains and MCMC algorithms.  Probability Surveys, 1 20–71.
  * Roberts and Sahu (1997) Roberts, G. O. and Sahu, S. K. (1997).  Updating schemes, correlation structure, blocking and parameterization for the Gibbs sampler.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 59 291–317.
  * Román (2012) Román, J. C. (2012).  Convergence analysis of block Gibbs samplers for Bayesian general linear mixed models.  Ph.D. thesis, University of Florida.
  * Román and Hobert (2012) Román, J. C. and Hobert, J. P. (2012).  Convergence analysis of the Gibbs sampler for Bayesian general linear mixed models with improper priors.  The Annals of Statistics, 40 2823–2849.
  * Roy (2012a) Roy, V. (2012a).  Convergence rates for MCMC algorithms for a robust Bayesian binary regression model.  Electronic Journal of Statistics, 6 2463–2485.
  * Roy (2012b) Roy, V. (2012b).  Spectral analytic comparisons for data augmentation.  Stat. and Prob. Letters, 82 103–108.
  * Roy (2020) Roy, V. (2020).  Convergence diagnostics for Markov chain Monte Carlo.  Annual Review of Statistics and Its Application, 7 387–412.
  * Roy and Hobert (2007) Roy, V. and Hobert, J. P. (2007).  Convergence rates and asymptotic standard errors for Markov chain Monte Carlo algorithms for Bayesian probit regression.  Journal of the Royal Statistical Society: Series B, 69 607–623.
  * Turek et al. (2017) Turek, D., de Valpine, P., Paciorek, C. J. and Anderson-Bergman, C. (2017).  Automated parameter blocking for efficient Markov chain Monte Carlo sampling.  Bayesian Analysis, 12 465–490.
  * van Dyk and Meng (2001) van Dyk, D. A. and Meng, X.-L. (2001).  The art of data augmentation (with discussion).  Journal of Computational and Graphical Statistics, 10 1–50.
  * Vats et al. (2018) Vats, D., Flegal, J. M. and Jones, G. L. (2018).  Strong consistency of multivariate spectral variance estimators in Markov chain Monte Carlo.  Bernoulli, 24 1860–1909.
  * Vats et al. (2019) Vats, D., Flegal, J. M. and Jones, G. L. (2019).  Multivariate output analysis for Markov chain Monte Carlo.  Biometrika, 106 321–337.
  * Wang and Roy (2018a) Wang, X. and Roy, V. (2018a).  Analysis of the Pólya-Gamma block Gibbs sampler for Bayesian logistic linear mixed models.  Statistics & Probability Letters, 137 251–256.
  * Wang and Roy (2018b) Wang, X. and Roy, V. (2018b).  Convergence analysis of the block Gibbs sampler for Bayesian probit linear mixed models with improper priors.  Electronic Journal of Statistics, 12 4412–4439.
  * Wang and Roy (2018c) Wang, X. and Roy, V. (2018c).  Geometric ergodicity of Pólya-Gamma Gibbs sampler for Bayesian logistic regression with a flat prior.  Electronic Journal of Statistics, 12 3295–3311.