Convergence analysis of block Gibbs samplers for Bayesian
probit linear mixed models

arXiv:1706.01846v3 [math.ST] 10 Aug 2017

Xin Wang∗ and Vivekananda Roy†
Department of Statistics, Iowa State University, Ames, IA
August, 2017

Abstract
In this article, we consider Markov chain Monte Carlo (MCMC) algorithms for
exploring the intractable posterior densities associated with Bayesian probit linear
mixed models under both proper and improper priors on the regression coefficients
and variance components. In particular, we construct two-block Gibbs samplers using the data augmentation (DA) techniques. Furthermore, we prove the geometric
ergodicity of the Gibbs samplers, which is the foundation for building the central limit theorems for MCMC based estimators and subsequent inferences. Under
proper priors, the Gibbs sampler chain is practically always geometrically ergodic.
While under improper priors, the conditions for geometric convergence are similar
to those guaranteeing posterior propriety. We also provide conditions for posterior propriety when the design matrices take commonly observed forms. The Haar
parameter expansion for DA (PX-DA) algorithm is an improvement of the DA algorithm and it has been shown that it is theoretically at least as good as the DA
algorithm. We propose corresponding Haar PX-DA algorithms, which have essentially the same computational cost as the two-block Gibbs samplers. An example is
used to show the efficiency gain of the Haar PX-DA algorithm over the block Gibbs
sampler.

key words: Data augmentation, Drift condition, Geometric ergodicity, Haar PXDA algorithm, Probit linear mixed models, Posterior propriety

1

Introduction

Generalized linear mixed models (GLMMs) are generalized linear models with random terms in the linear predictor. The random effects in the GLMM can accommodate
for overdispersion often present in non-Gaussian data, and dependence among correlated
∗
†

Email: xinwang@iasate.edu
Email: vroy@iastate.edu

1

observations arising from longitudinal or repeated measures studies. GLMM is one of the
most frequently used statistical models. Here, we consider a popular Bayesian GLMM
for binary data, namely, the probit linear mixed model.
Let (Y1 , Y2 , . . . , Yn ) denote the vector of Bernoulli random variables, and let xi and
zi be the p × 1 and q × 1 known covariates and random effect design vectors respectively
associated with ith observation for i = 1, . . . , n. Let β ∈ Rp be the unknown vector of
regression coefficients and u ∈ Rq be the random effects vector. A GLMM can be built
(McCulloch et al., 2011) with a link function that connects the expectation of Yi with
the covariates xi , zi . One of the very popular link functions is the probit link function
Φ−1 , that is,
P (Yi = 1) = Φ(xTi β + ziT u),
where Φ is the cumulative distribution function of the standard normal random variable.
Assume that we have r random effects with u = (uT1 , . . . , uTr )T , where uj is a qj × 1
ind

vector with qj > 0, q1 + · · · + qr = q, and uj ∼ N (0, Iqj 1/τj ), where τj ∈ R+ = (0, ∞) is
the precision parameter associated with uj for j = 1, . . . , r. Let τ = (τ1 , . . . , τr ), thus,
the data model for the probit GLMM is
Yi |β, u

ind

Bern(αi ) for i = 1, . . . , n with

αi

=

uj |τ

ind

Φ(xTi β + ziT u) for i = 1, . . . , n,

 1
N 0, Iqj , j = 1, . . . , r.
τj

∼

∼

(1)

Let y = (y1 , y2 , . . . , yn )T be the observed Bernoulli response variables. Note that, the
likelihood function for (β, τ ) is
L(β, τ |y) =

Z

Rq

n h
Y

iyi h

Φ(xTi β + ziT u)

i=1

i1−yi

1 − Φ(xTi β + ziT u)

φq (u; 0, D −1 )du, (2)

which is not available in closed form, but only as a high dimensional integral. Here,
φq (s; a, B) denotes the probability density function of the q−dimensional normal distribution with mean vector a, covariance matrix B and evaluated at s, where D =
⊕rj=1 τj Iqj .
In Bayesian framework, one needs to specify the prior distributions of β and τ .
Assume β and τ are apriori independent. Let π(β) and π(τ ) be the prior densities of β
and τ respectively. Thus, the joint posterior density of (β, τ ) is
π(β, τ |y) =
where
c(y) =

Z

Rr+

Z

1
L(β, τ |y)π(β)π(τ ),
c(y)

Rp

(3)

L(β, τ |y)π(β)π(τ )dβdτ ,

is the marginal density of y. Since the likelihood function L(β, τ |y) is not available in
closed form, the posterior density is intractable for any choice of the prior distributions
2

of β and τ . In this article, we consider both of the two popular choices of priors for β,
namely a normal prior and the improper flat prior. When the prior distribution of β is
normal (section 2), we do not assume any specific form of the prior distribution of τ .
On the other hand, if π(β) ∝ 1 (section 3), our results hold under popular gamma priors
and power priors on τ . Generally, Markov chain Monte Carlo (MCMC) algorithms are
used for exploring these posterior densities.
Even in the absence of random effects, for the probit regression model, the posterior
distribution of β is difficult to sample from (Roy and Hobert, 2007). Albert and Chib’s
(1993) MCMC algorithm for sampling from the posterior distribution associated with
the probit regression model is the most widely used data augmentation (DA) algorithm.
Albert and Chib’s (1993) DA technique can be also applied to the probit linear mixed
model. Following Albert and Chib (1993), let vi ∈ R be the continuous latent variable
ind
corresponding to binary observation Yi , such that Yi = I(vi > 0), where vi |β, u ∼
N (xTi β + ziT u, 1) for i = 1, . . . , n. Then
P (Yi = 1) = P (vi > 0) = Φ(xTi β + ziT u),

(4)

ind

that is, Yi |β, u ∼ Bern(αi ) as in (1). Let v = (v1 , . . . , vn )T , then v|β, u ∼ N (Xβ +
Zu, In ), where Xn×p = (x1 , . . . , xn )T and Zn×q = (z1 , . . . , zn )T .
Using the latent variables v, we can introduce a joint density π(β, u, v, τ |y) (see
section 2.1 for details) such that
Z

Rq

Z

Rn

π(β, u, v, τ |y)dvdu = π(β, τ |y),

(5)

where π(β, τ |y) is the posterior density defined in (3). If all the full conditionals of
the joint density π(β, u, v, τ |y) are easy to sample from, then a Gibbs sampler can be
run and it can be used to make inference on the posterior density (3). Indeed this
full conditional Gibbs sampler is traditionally used in the analysis of Bayesian probit
linear mixed models (Baragatti, 2011). In this article, instead of using full conditional
distributions, we construct a two-block Gibbs sampler with (β T , uT )T as one block
and (v T , τ T )T as the other block, which is our first contribution. In general, block
Gibbs samplers are more efficient than the Gibbs samplers based on full conditional
distributions (Liu et al., 1994).
The above mentioned block Gibbs sampler has an everywhere strictly positive Markov
transition density, implying the underlying Markov chain is Harris ergodic (Asmussen and Glynn,
2011; Meyn and Tweedie, 1993). Thus, the time average estimators based on the block
Gibbs sampler can be used to consistently estimate the (posterior) means with respect
to the joint density π(β, u, v, τ |y). In practice, it is crucial to know whether the Monte
Carlo errors associated with these estimates are sufficiently small. However, in order to
provide valid standard errors, we need to establish a central limit theorem (CLT) for the
time average estimators. Unlike, in the ordinary Monte Carlo methods based on iid samples, mere existence of the finite second moment does not guarantee a CLT for MCMC
estimators. The only standard method of establishing CLT for MCMC estimators, is
to prove that the underlying Markov chain is geometrically ergodic (Jones and Hobert,
3

2001). Geometric ergodicity is also needed for consistently estimating the asymptotic
variance in the Markov chain CLT (Flegal and Jones, 2010). Roy and Hobert (2007) and
Chakraborty and Khare (2017) proved the geometric ergodicity of Albert and Chib’s
(1993) DA algorithm for the Bayesian probit regression model under improper and proper
priors on the regression coefficients. For linear models, Jones and Hobert (2004) analyzed the Gibbs sampler for one-way random effects models under proper priors and
Tan and Hobert (2009) extended it to improper priors case. Johnson and Jones (2010)
analyed the block Gibbs sampler for Bayeisan linear mixed models under an assumption
X T Z = 0. Román and Hobert (2012) and Román and Hobert (2015) analyzed convergence properties of the Gibbs samplers for Bayesian linear mixed models under improper
and proper priors without the assmuption of X T Z = 0. Doss and Hobert (2010) established geometric convergence for the Gibbs sampler for a linear mixed model where
random effects have a t distribution. However, no such convergence analysis of the Gibbs
samplers for any Bayesian generalized linear mixed model is available in the literature.
Our second contribution, in this paper, is establishing geometric convergence rates for
the block Gibbs sampler for Bayesian probit linear mixed models.
DA algorithms are known to suffer from slow convergence (Meng and van Dyk, 1999;
van Dyk and Meng, 2001). Liu and Wu (1999) proposed the parameter expansion for
data augmentation (PX-DA) algorithm, which can converge faster than the DA algorithm without much extra computational effort (van Dyk and Meng, 2001; Roy, 2014).
Hobert and Marchev (2008) proved that the Haar PX-DA algorithm, that is based on
the Haar measure, is better than any other PX-DA algorithms and the original DA algorithm in both efficiency and operator norm ordering. In the probit regression model,
Roy and Hobert (2007), through an example, showed that PX-DA algorithms can lead
to huge gains in efficiency over Albert and Chib’s (1993) DA algorithm. Our third contribution is to construct Haar PX-DA algorithms improving the block Gibbs samplers
mentioned before. Since the geometric ergodicity of the Haar PX-DA algorithm follows
from geometric ergodicity of the DA algorithm (Hobert and Marchev, 2008), we have
CLTs for the Haar PX-DA algorithm based estimators.
The article is organized as follows. In section 2, we construct the two-block Gibbs
sampler for the Bayesian probit linear mixed model under proper priors, provide its
geometric convergence results and the corresponding Haar PX-DA algorithm. In section
3, we construct a similar two-block Gibbs sampler under improper priors and prove
the geometric ergodicity of the underlying Markov chain. We also establish conditions
for propriety of the posterior distribution under improper priors, when X and Z take
commonly observed forms. A real data example is given in section 4. Section 5 contains
some conclusions and discussions. Finally, the proofs of the geometric ergodicity of the
Gibbs samplers and posterior propriety appear in the appendices.

4

2

Geometric ergodicity of the two block Gibbs sampler
under proper priors

In this section, we first construct a two block Gibbs sampler for the Bayesian probit
linear mixed model in section 2.1. Then, in section 2.2, we describe its geometric rate
of convergence. And in section 2.3, the proposed Haar PX-DA algorithm is presented.

2.1

A two block Gibbs sampler

We begin with deriving the joint density π(β, u, v, τ |y) mentioned in the introduction. Given τ , we assume that u and β are independent with
β|τ ∼ Np (Q−1 µ0 , Q−1 ),
where Q is a p × p known positive definite matrix, µ0 ∈ Rp is known. As mentioned in
Chakraborty and Khare (2017), any vector µ̃ ∈ Rp can be written as µ̃ = Q−1 Qµ̃ =
Q−1 µ0 with µ0 = Qµ̃. Thus, it is not restrictive to assume the above form of the prior
mean. Define
n
h
iyi h
i1−yi
Y
1
π(β, u, v, τ |y) =
φ(vi ; xTi β + ziT u, 1) 1(0,∞) (vi )
1(−∞,0] (vi )
c(y) i=1

"

× φq (u; 0, D −1 )φp (β; Q−1 µ0 , Q−1 )π(τ ).

#

(6)

From (2) and (4) it follows that (5) holds. Since
n 
n h
2 Y
i yi h
i1−yi
1X
vi − xTi β − ziT u
·
1(0,∞) (vi )
1(−∞,0] (vi )
π(β, u, v, τ |y) ∝ exp −
2 i=1
i=1

#

"

·

r
Y

qj
2

τj

j=1


T 
1
1
· exp − τj uTj uj · exp − β − Q−1 µ0 Q β − Q−1 µ0 · π(τ ),
2
2






letting η ≡ (β T , uT )T and Wn×(p+q) ≡ (X, Z), standard calculations show that the
conditional density of η is
1
1
π(η|v, τ , y) ∝ exp − (v − W η)T (v − W η) · exp − η T A(τ )η + η T l ,
2
2








(7)

= W T W + A(τ ) and µ = W T v + l.

(8)

where
l=
Thus,

T
µT0 , 0Tq×1



and A(τ ) =


Q 0p×q
0q×p D

!

.



η|v, τ , y ∼ Np+q Σ−1 µ, Σ−1 ,
where
Σ=

XT X + Q
XT Z
T
T
Z X
Z Z +D

!

5



Similarly, the conditional density of (v, τ ) is
π(v, τ |η, y) ∝

n
Y

i=1
r
Y

·

iyi h



h

τj

1
· exp − τj uTj uj · π(τ ),
2

φ vi − wiT η; 0, 1
qj
2

1(0,∞) (vi )



j=1

1(−∞,0] (vi )

i1−yi



where wiT is the ith row of W for i = 1, . . . , n. Thus, conditional on (η, y), vi , i =
1, . . . , n and τ are independent. We have
ind

vi |η, y ∼ TN(wiT η, 1, yi ), i = 1, . . . , n,

(9)

where TN(µ, σ 2 , ω) denotes the distribution of the normal random variable with mean
µ and variance σ 2 , that is truncated to have only positive values if ω = 1, and, it has
only negative values if ω = 0.
The conditional density of τ is
π(τ |η, y) ∝

r
Y

qj
1
τj 2 exp − τj uTj uj · π(τ ).
2
j=1





(10)

In the special case, if π(τ ) = rj=1 π(τj ), and for each j, π(τj ) is the density of the
gamma random variable with parameters aj and bj , that is
Q

a −1 −bj τj

π(τj ) ∝ τj j

e

,

(11)

then conditional on η, y, τj ’s are independent and τj follows Gamma distribution with
parameter aj + qj /2 and bj + uTj uj /2.
Thus, one single iteration of the block Gibbs sampler {η (m) , v (m) , τ (m) }∞
m=0 has the
following two steps:
Algorithm 1 The (m + 1)st iteration for the two block Gibbs sampler
(m+1)

Draw τ (m+1) from (10) with u = u(m) , and independently draw vi
TN(wiT η (m) , 1, yi ) for i = 1, . . . , n.
2: Draw η (m+1) from (7), that is,
1:

η (m+1) ∼ N

2.2



−1 

W T W + A(τ (m+1) )

 

ind

|η (m) , y ∼

−1 

W T v (m+1) + l , W T W + A(τ (m+1) )

Geometric ergodicity for the block Gibbs sampler

In this section, we establish the geometric convergence of the block Gibbs sampler
Since it is a two-block Gibbs sampler, it has the same rate of

{η (m) , v (m) , τ (m) }∞
m=0 .

6

.

convergence as the η-marginal Markov chain {η (m) }∞
m=0 (Roberts and Rosenthal, 2001).
chain.
Below we analyze the Ψ ≡ {η (m) }∞
m=0
Let η ′ be the current state and η be the next state, then the Markov transition
density (Mtd) of Ψ is
′

k(η|η ) =

Z

Rr+

Z

Rn

π(η|v, τ , y)π(v, τ |η ′ , y)dvdτ ,

(12)

where π(·|·, y)’s are the conditional densities from section 2.1. Routine calculations show
that k(η|η ′ ) is reversible
and thus invariant with respect to the marginal density of η deR R
noted as π(η|y) ≡ Rr Rn π(η, v, τ |y)dvdτ . Let h : Rp+q → R and suppose our interest
+
R
is to estimate the (posterior) mean E(h(η)|y) = Rp+q h(η|π(η|y)dη. Since k(η|η ′ ) is
strictly positive and thus Harris ergodic (Meyn and Tweedie, 1993)), if E(|h(η)||y) < ∞,
then E(h(η)|y) can be consistently estimated by
h̄m =

X
1 m−1
h(η (i) ).
m i=0

As mentioned in the introduction, in order to provide an asymptotically valid confidence
interval for E(h(η)|y) based on h̄m , we need to establish a CLT for it. We say a CLT
exists for h̄m if there exists a constant σh2 ∈ (0, ∞) such that,
√





d





m h̄m − E(h(η)|y) → N 0, σh2 as m → ∞.

(13)

If (13) holds and a consistent estimator σ̂h2 of σh2 is available, then the standard er√
rors σ̂h / m can be used to provide an asymptotic confidence interval for E(h(η)|y)
(Roy and Hobert, 2007). Harris ergodicity does not guarantee (13), although it ensures
consistency of h̄m . One method of proving (13) is to establish the geometric ergodicity of the Markov chain Ψ (Jones and Hobert, 2001). Geometric ergodicity of Ψ also
allows for consistent estimation of σh2 using batch means or spectral variance methods
(Flegal and Jones, 2010).
Let B denote the Borel σ-algebra of Rp+q and K(·, ·) be the Markov transition
function corresponding to the Mtd k(·, ·) in (12), that is, for any set A ∈ B, η ′ ∈ Rp+q
and any j = 0, 1, . . . ,
K(η ′ , A) = Pr(η (j+1) ∈ A|η (j) = η ′ ) =

Z

k(η|η ′ )dη.

(14)

A

Then the m-step Markov transition function is K m (η ′ , A) = Pr(η (m+j) ∈ A|η (j) = η ′ ).
Let Π(·|y) be the probability measure with density π(η|y). The Markov chain Ψ is
geometrically ergodic if there exists a constant 0 < t < 1 and a function G : Rp+q 7→ R+
such that for any η ∈ Rp+q ,
||K m (η, ·) − Π(·|y)|| := sup |K m (η, A) − Π(A|y)| ≤ G(η)tm .
A∈B

7

(15)

Harris ergodicity of Ψ implies that ||K m (η, ·) − Π(·|y)|| ↓ 0 as m → ∞, while (15)
guarantees its exponential rate of convergence. Roberts et al. (1997) show that since Ψ is
reversible, if (15) holds then there exists a CLT, that is (13), for all h with E(h2 (η)|y) <
∞. A proof of Theorem 1 by establishing a drift condition for Ψ is given in the Appendix
A.1 .
Theorem 1. Assume, there exists a constant τ0 > 0 such that τj ≥ τ0 for all j, and
the conditional density π(τ |η, y) defined in (10) is a valid probability density for any u.
Then the Markov chain Ψ is geometrically ergodic.
Remark 1. Theorem 1 holds if the prior distributions for τj ’s are truncated below a
given value τ0 . However, while implementing the block Gibbs sampler in practice, a
number slightly larger than the machine precision zero can be treated as τ0 , practically
avoiding the need to use any rejection sampling algorithms to draw from the truncated
conditional distribution of τ . For example, in R (R Development Core Team, 2011),
this machine precision zero value is 2.220446−16 .
Remark 2. Theorem 1 assumes that π(τ |η, y) is a valid probability density for all u.
a −1
For the important special case when π(τj ) ∝ τj j e−bj τj , this implies that qj + 2aj > 0
and bj > 0 for all j.
Remark 3. Theorem 1 does not make any assumptions on the dimensions p, q or the
sample size n. Thus, it is applicable to the modern high dimensional datasets where p
(and/or q) can be much larger than n.

2.3

Haar PX-DA algorithm

As mentioned in section 1, DA algorithms often suffer from slow convergence and
high autocorrelations. Liu and Wu (1999) proposed parameter expansion for data augmentation (PX-DA) algorithms for speeding up the convergence of DA algorithms.
Hobert and Marchev (2008) compared the performance of PX-DA algorithms based on
Haar measure (called Haar PX-DA algorithms) with PX-DA algorithms based on a
probability measure and DA algorithms. In particular, they showed that, under some
mild conditions, the Haar PX-DA algorithms are better than general PX-DA algorithms and DA algorithms in both efficiency and operator norm ordering. As shown
in Hobert and Marchev (2008), compared to the DA algorithm, in PX-DA, an extra
step is added (sandwiched) between the two steps of the original DA algorithm. In order
to construct this extra step, we derive the marginal density
π (v, τ |y) =
∝
·

Z

Rp+q
n
Yh

π(η, v, τ |y)dη

1(0,∞) (vi )

i=1

iyi h

1(−∞,0] (vi )

(16)
r
i1−yi Y

qj

τj 2

j=1

i
1h
exp − v T M1 v − 2v T M2 − lT Σl π(τ ),
2




8

where




−1

T

M1 = I − W W W + A(τ )

W

T



−1



and M2 = W W T W + A(τ )

l.

Let Z denote the subset of Rn where v lives, that is, Z is the Cartesian product of n
half (positive or negative) lines, where the ith component is (0, ∞) (if yi = 1) or (−∞, 0]
(if yi = 0). Let G be the unimodular multiplicative group on R+ with Haar measure
ν(dg) = dg/g, where dg is Lebesgue measure on R+ . For constructing an efficient
extra step, as in Roy (2014), we let the group G act on Z × Rr+ through a group action
T (v, τ ) = (gv, τ ) = (gv1 , gv2 , . . . , gvn , τ ). With the group action defined this way, it can
be shown that the Lebesgue measure on Z ×Rr+ is relatively left invariant with multiplier
χ(g) = gn (Roy, 2014; Hobert and Marchev, 2008). Following Hobert and Marchev
(2008), consider a probability density function q(g) on G where
q (g) dg ∝ π (gv, τ |y) χ (g) ν(dg)

i
1h
∝ gn−1 exp − g2 v T M1 v − 2gv T M2 dg.
2

(17)

Given (v, τ ), q (g) is a valid density since M1 is a positive definite matrix. From
Hobert and Marchev (2008), it follows that the transition (v, τ ) → (v ′ , τ ) ≡ T (v, τ ) =
(gv, τ ) where g ∼ q(g), is reversible with respect to π(v, τ |y) defined in (16). Below are
the three steps involved in every iterations of the Haar PX-DA algorithm.
Algorithm 2 The (m + 1)st iteration for the Haar PX-DA algorithm
ind

Draw τ from (10) with u = u(m) , and independently draw vi |η, y ∼
TN(wiT η (m) , 1, yi ), i = 1, . . . , n.
2: Draw g from (17).
3: Calculate vi′ = gvi for i = 1, . . . , n, and draw η (m+1) from (7) conditional on v ′ =
(v1′ , . . . , vn′ ) and τ , that is, draw

1:

η

(m+1)



∼N

−1 

T

W W + A(τ )

T

 

′

−1 

T

W v + l , W W + A(τ )

.

The Mtd of the above Haar PX-DA algorithm can be written as
k∗ η|η ′ =


Z

Rn

Z

Rr+

Z

Rn

π η|v ′ , τ , y F v, dv ′ π v, τ |η ′ , y dvdτ dv ′ .






(18)

where F (·, ·) is the Markov transition function corresponding to the move (v, τ ) →
(v ′ , τ ) = T (v, τ ). Let K ∗ and K be the Markov operators associated with the Mtds k∗
and k defined in (18) and (12) respectively. From Hobert and Marchev (2008), we have
||K ∗ || ≤ ||K|| (see also Roy, 2012b). Since the block Gibbs sampler is geometrically
ergodic, we have ||K ∗ || ≤ ||K|| < 1 (Roberts et al., 1997). Thus we have the following
corollary.
9

Corollary 1. Under the conditions of Theorem 1, the Markov chain underlying the Haar
PX-DA algorithm described in Algorithm 2 is geometrically ergodic.
Since q(g) in (17) is log-concave, adaptive rejection sampling algorithm (Gilks and Wild,
1992) is used to efficiently sample from q(g). The extra step of Algorithm 2 is a single
draw from the univariate density q(g), which is easy to sample from. Thus, the computational burden, per iteration, for the Haar PX-DA algorithm is similar to that of
the block Gibbs sampler described in section 2.1. Two other Haar PX-DA algorithms
can be constructed by using group actions T1 (v, τ ) = (v, gτ ) and T2 (v, τ ) = (gv, gτ ).
However, the corresponding q(g)’s are not easy to sample from, thus we do not consider
them here.

3

Geometric ergodicity of the two block Gibbs sampler
under improper priors

In this section we consider an improper prior for β, that is π(β) ∝ 1 and τj ’s are
a −1
apriori independent with (11), that is, π(τj ) ∝ e−bj τj τj j , which can be proper and
improper. Since we are using improper priors, the posterior density is not guaranteed
to exist. In section 3.1, we discuss the conditions under which this posterior density is
proper.

3.1

Geometric ergodicity of the Gibbs sampler

In this case, the joint posterior density (up to a normalizing constant) of β, u, v, τ ,
if it exists, is
"

n
Y

2  h
iyi h
i1−yi
1
T
T
1(0,∞) (vi )
1(−∞,0] (vi )
exp − vi − xi β − zi u
π (β, u, v, τ |y) ∝
2
i=1
∗

×

r
Y



qj
2

τj

+aj −1

j=1

(

exp −τj

uTj uj
bj +
2

!)

.

#

(19)

A similar two block Gibbs sampler as in section 2.1 can be considered in this setting
by alternating draws from the two conditionals η|v, τ , y and v, τ |η, y. From (19), it
follows that


η|v, τ , y ∼ Np+q Σ∗−1 W T v, Σ∗−1 ,
where

∗

Σ =
and
Σ∗−1

XT X
XT Z
T
T
Z X Z Z +D

 
−1
TX
+ RS −1 RT
X
=

−S −1 RT
10

!

,


−RS −1 
,
S −1

with S = Z T (I − PX ) Z + D, R =
Thus



XT X

−1



X T Z and PX = X X T X

−1

XT .

 
−1
h
i 
T I − ZS −1 Z T (I − P ) v
TX
X
X
X
.
Σ∗−1 W T v = 

S −1 Z T (I − PX ) v

The conditional distribution for vi conditional on η, y is the same as in Section
2.1


provided in (9). And conditioned on η, y, τj ’s are independent with τj ∼ Gamma aj + qj /2, bj + uTj uj /2

for j = 1, . . . , r. Thus, one single iteration of the block Gibbs sampler {η (m) , v (m) , τ (m) }∞
m=0
has the following two steps:
Algorithm 3 The (m + 1)st iteration for the two block Gibbs sampler
1:



(m)T

Draw τ (m+1) ∼ Gamma aj + qj /2, bj + uj

(m)

uj



/2

ind
(m+1) (m)
T (m)
vi
|η , y ∼ TN(w
 i η , 1, yi ), i = 1, . . . , n.

2: Draw η (m+1) ∼ Np+q Σ∗(m+1)−1 W T v (m+1) , Σ∗(m+1)−1 ,
ated at τ (m+1) .

and independently draw
where Σ∗(m+1) is evalu-

We now discuss conditions under which the posterior density (19) is proper. Let
∗
ci = 1 if yi = 0 and ci = −1 if yi = 1 for i = 1, . . . , n. Suppose Wn×(p+q)
is a matrix
T
T
with ith row ci (xi , zi ). In the special case when bj = 0, that is when τj has the power
a −1
prior π(τj ) ∝ τj j , Chen et al. (2002) (Theorem 4.2) shown that the corresponding
posterior distribution is proper under the following conditions:
(A1) W = (X, Z) is a full rank matrix;
(A2) There exists a n × 1 positive vector e > 0 such that eT W ∗ = 0;
(A3) min{2a1 + q1 , . . . , 2ar + qr } > 0;
(A4) aj < 0 for j = 1, . . . , r.
Roy and Hobert (2007) provide a simple method of checking the condition A2 using
publicly available softwares.
Here we consider a general form of the prior distribution of τj , that is π(τj ) ∝
a −1
−b
j
e τj τj j . The following theorem gives the conditions for geometric ergodicity of the
Markov chain in algorithm 3 and the proof is given in the appendix A.2.
Theorem 2. The Markov chain underlying the block Gibbs sampler is geometrically
ergodic if the following conditions hold:
(1) aj < bj = 0 or bj > 0 for j = 1, . . . , r;
(2) (A1) − (A3) hold.
11

The condition A1 assumes W is a full rank matrix. Unfortunately, when Z is a
design matrix with elements 1 and 0, which is pretty common in practice, such as in
two way random effect models, this assumption is not reasonable. In this case, we
provide Theorems 3 and 4 showing the posterior propriety and geometric convergence
of the Gibbs sampler without the assumption A1. Assume X is a matrix with 1 as the
first column corresponding to the intercept term and Z = (Z1 , . . . , Zr ), where Zj is a
n × qj matrix such that the ikth element is 1 if observation i = 1, . . . , n is observed at
level k = 1, . . . , qj of random effect uj , 0 otherwise. We have β = (β0 , β1 , . . . , βp−1 )T
and uj =



T

uj1 , . . . , ujqj

. Let µ0 = β0 +

Pr

j=1 uj1 ,

djk = uj,k+1 − uj1 , for k =

1, . . . , qj −1, j = 1, . . . , r and η̃ = (µ0 , β1, , . . . , βp−1 , d11 , . . . , d1,q1 −1 , . . . dr1 , . . . , dr,qr −1 )T .
Define Z̃ = (Z̃1 , . . . , Z̃r ), where Z̃j contains the columns

 of Zj without the first column.
Thus, W η can be written as W̃ η̃, where W̃ = X, Z̃ . Let W̃ ∗ be a matrix with ith

row w̃iT = ci (xTi , z̃iT ), where z̃iT is the ith row of Z̃.
Theorem 3. Assume the following conditions hold,
(B1) aj < bj = 0, qj ≥ 2 or bj > 0 for j = 1, . . . , r;
(B2) 2aj + qj − 1 > 0 for j = 1, . . . , r;
(B3) W̃ is a full rank matrix;

(B4) There exists a n × 1 positive vector e > 0 such that eT W̃ ∗ = 0.
Then the joint posterior distribution (19) is proper, i.e,
Z

Rr+

·

Z

r
Y

Z

Rq

Rp i=1

qj
2

τj

j=1

n
Y

+aj −1



Φ xTi β + ziT u


exp −τj

yi h



1 − Φ xTi β + ziT u

1
bj + uTj uj
2





i1−yi

dβdudτ < ∞.

(20)

The following theorem provides conditions for geometric convergence of the Markov
chain underlying the Gibbs sampler.
Theorem 4. The block Gibbs sampler is geometrically ergodic under the following conditions:
(1) (B1) − (B4) hold;
(2) There exists an s ∈ (0, 1] ∩ (0, s̃) such that
2−s

r
X
Γ (qj /2 + aj − s) h

j=1

Γ (qj /2 + aj )







tr Rj I − PZ T (I−PX )Z RjT

is

< 1.

where s̃ = min{a1 + q1 /2, . . . , ar + qr /2}, Rj is a qj × q matrix with 0 and 1
such that Rj u = uj and PZ T (I−PX )Z is the projection on the column space of
Z T (I − PX ) Z.
12

The proofs of Theorems 3 and 4 are given in Appendix A.3 and Appendix A.4
respectively.

3.2

Haar PX-DA algorithm

Under improper priors, we can also construct a Haar PX-DA algorithm as in Section
2.3. From (19), we derive the marginal density of v, τ given by
∗

π (v, τ |y) =
∝

Z

Rp+q
n h
Y

π ∗ (η, v, τ |y)dη
r
i1−yi Y

iyi h

1(0,∞) (vi )

i=1

(21)

1(−∞,0] (vi )

+aj −bj τj

e

j=1

1
· exp − v T M4 v ,
2


qj

τj 2



where M4 = I − W Σ∗−1 W T . As in section 2.3, we consider a probability density
function q ∗ (g) on G where
q ∗ (g) dg ∝ π ∗ (gv, τ |y) χ (g) ν(dg)


1
∝ gn−1 exp − g2 v T M4 v dg.
2

(22)

Since q ∗ (·) is a log concave function, we can use adaptive rejection sampling to
effectively sample from it. Thus one step of the Haar PX-DA algorithm has the following
3 steps.
Algorithm 4 The (m + 1)st iteration for the Haar PX-DA algorithm
1:



(m)T

τ ∼ Gamma aj + qj /2, bj + uj

(m)

uj



/2

ind

and independently draw vi |η, y ∼

TN(wiT η (m) , 1, yi ) for i = 1, . . . , n.
2: Draw g from (22).


3: Calculate vi′ = gvi and draw η (m+1) ∼ Np+q Σ∗−1 W T v ′ , Σ∗−1 .
Similar as section 2.3, we have the following corollary.

Corollary 2. Under the conditions of Theorem 2 or Theorem 4, the Markov chain
underlying the Haar PX-DA algorithm described in Algorithm 4 is geometrically ergodic.

4

An example

In this section, we use the bacteria data set from Venables and Ripley (2002) to
compare the performance of Algorithm 1 and Algorithm 2 mentioned in sections 2.1 and
2.3 respectively. The response variable here is a binary variable indicating the presence
or absence of the bacteria H.influenzae. There are three covariates: week of test, three
13

Table 1: Results for the block Gibbs sampler and the Haar PX-DA algorithm
DA algorithm
PX-DA algorithm
Parameter estimate
s.e.
estimate
s.e
ratio
β0
2.2404
0.0015
2.2383
0.0009 2.6487
β1
-0.9795 0.0007 -0.9786
0.0005 1.8586
-0.8502 0.0009 -0.8496
0.0007 1.6401
β2
β3
-0.5183 0.0009 -0.5172
0.0007 1.4602
treatments (a placebo, and an active drug with and without extra effort to ensure that
it is taken), and the subject ID. Venables and Ripley (2002) used the logit linear mixed
model to analyze this data. Instead of the logit link function, we use the probit link.
The fixed week effect enters the model as I(week > 2), where I(·) is a indicator function.
The model we consider is
P (Yij = 1) = Φ (β0 + β1 I(weekij > 2) + β2 x2i + β3 x3i + ui ) ,
for i = 1, . . . , n = 50, j = 1, . . . , ni , where ni is the number of observations for subject
i, x2 and x3 are binary variables corresponding to active drug treatment without and
with extra effort to guarantee the drug is taken and u = (u1 , . . . , un )T is the random
effect with 50 levels with u ∼ N (0, 1/τ I50 ). We consider proper priors for both β and
τ . We assume βk ’s are apriori independent and the prior distribution of βk is normal
with mean 0 and variance 1000 for k = 0, 1, 2, 3. For the prior distributions of τ , we
ind
assume τ ∼ Gamma(a, b). We tried priors with parameters a = b = 1, a = b = 0.1
and a = b = 0.01. As mentioned in Remark 1, we set τ0 as a value which is larger than
the machine precision value in R (such as 10−10 ). We compare the block Gibbs sampler
and the Haar PX-DA algorithm using the standard errors for the parameter estimates.
We use the spectral variance estimator with Tukey-Hanning window (Flegal and Jones,
2010) for estimating the standard errors. Markov chains of length 2 × 106 after an initial
burn-in of 2 × 106 iterations are used to calculate the parameter estimates and standard
error estimates. Table 1 shows estimates of the parameters and corresponding standard
errors based on a = b = 1. The two algorithms have similar parameter estimates. We
also calculate the ratio of variance estimates from the DA algorithm and the Haar PX-DA
algorithm. The asymptotic variance ratio estimates (σ̂h2 /σ̂h∗2 ) for different parameters
can be as high as 2.6. That means, the block Gibbs sampler needs about 2.6 times
as many iterations as the Haar PX-DA algorithm in order to achieve the same level of
precision. Since the two algorithms basically require same amount of time per iteration,
use of the Haar PX-DA algorithm can lead to significant gain in efficiency. Table 2
shows the variance ratio estimates for other two priors. The ratios are still larger than
1, the ratios decrease with smaller values of a and b. Using the method describe in
Roy and Hobert (2007), we see that the condition (B4) does not hold in this example.
Thus we do not use improper flat prior for β here.

14

Table 2: Variance ratio estimates for other two priors
a = b = 0.1 a = b = 0.01
β0
2.2435
1.8655
1.8151
1.6254
β1
1.6516
1.5369
β2
β3
1.4543
1.3312

5

Discussion

We develop two block Gibbs samplers for the Bayesian probit linear mixed models
under both proper priors and improper priors. The block Gibbs algorithm samples
the fixed effects and the random effects jointly. We prove the geometric ergodicity of
the two-block Gibbs samplers, which guarantee the existence of central limit theorem
for MCMC estimators under a finite second moment condition. This is the first ever
theoretical analysis of Gibbs samplers for any Bayesian generalized linear mixed model.
We propose the corresponding Haar PX-DA algorithms for the two cases, which not
only improves the efficiency of the Gibbs samplers, but also inherits their geometric
convergence property. Using an example, we show that the Haar PX-DA algorithm can
lead to much efficiency gain compared to the block Gibbs sampler.
Another popular link function is the logit link function. Polson et al. (2013) propose
a DA algorithm for the logistic regression model. Choi and Hobert (2013) prove the
uniform ergodicity of this DA algorithm. As mentioned in Polson et al. (2013), their DA
algorithm can be extended to the logistic linear mixed model. However, the convergence
properties of the corresponding Markov chain have not been studied, and can be a topic
for future research. Another future project can be deriving similar extensions of the
results in Roy (2012a) for proving geometric convergence of the Gibbs samplers for the
robit linear mixed model.

A
A.1

Appendices
Proof of Theorem 1

Proof. We prove the geometric ergodicity of Ψ by establishing a drift condition. In
particular, we consider the drift function
V (η) ≡
and show that for any

η, η ′

n 
X

xTi β + ziT u

i=1

2

+ β T Qβ + τ0 uT u,

∈ Rp+q , there exist ρ ∈ (0, 1) and L > 0 such that
E V (η)|η ′ ≤ ρV (η ′ ) + L.


Let
M=



XT X + Q
XT Z
T
T
Z X
Z Z + τ0 Iq
15

(23)
!

.

The drift function can be written as V (η) = η T M η. Now, we prove (23). By Fubini’s
theorem, we have


E V (η) |η

′

=
=

Z

p+q

V (η) k η|η ′ dη

ZR Z

Rr+

Rn

Z



Rp+q

V (η) π (η|v, τ , y) π v, τ |η ′ , y dηdτ dv.


So the expectation can be calculated in two steps. The first step is to calculate E (V (η) |v, τ , y).
Note that
Q
0
0 τ0 Iq

T

M = W W + A(τ0 ), where A(τ0 ) =

!

.

Since τj ≥ τ0 for j = 1, . . . , r, we have M ≤ Σ, where Σ is defined in (8). Thus,


E η T M η|v, τ , y





≤ E η T Ση|v, τ , y



= p + q + [E (η|v, τ , y)]T ΣE (η|v, τ , y)
= p + q + µT Σ−1 µ


= p + q + WTv + l
1



T



Σ−1 W T v + l


= p + q + ||Σ− 2 W T v + l ||2 .



Using the inequality (Chakraborty and Khare, 2017)


2

2

||a + b|| ≤ 1 + c



1
||a|| + 1 + 2 ||b||2 ,
c




2

for any a, b ∈ Rp+q and c > 0, we have
||Σ

− 21



T



2

W v + l ||



≤



=

−1

Since Σ−1 = W T W + A(τ )
Bernstein (2009)), we have





T

v WΣ



−1

1
W v + 1 + 2 lT Σ−1 l
c


T



−1

1 + c2 v T W W T W + A(τ )



WTv


−1
1 T T
1 + 2 l W W + A(τ )
l.
c

≤

−1



W W T W + A(τ )
and

1+c



+


2

−1



W T W + A(τ0 )

= M −1 (Proposition 8.6.6,



−1

W T ≤ W W T W + A(τ0 )
−1

lT W T W + A(τ )

16

l ≤ lT M −1 l.

WT,

Let X̃ = W A(τ0 )−1/2 , then
−1



W W T W + A(τ0 )



W T = X̃ X̃ T X̃ + In

−1

X̃ T .

From Chakraborty and Khare (2017)[Proposition A.1], we know that all eigenvalues of


X̃ X̃ T X̃ + In

−1

X̃ T lie in [0, 1). Let d1 , . . . , dt denote the nonzero singular values of

X̃, where t is the rank of W . Then from Chakraborty and Khare (2017)[Proposition
−1



WT.
A.1], d2i /(d2i +1), j = 1, . . . , t are the nonzero eigenvalues of W W T W + A(τ0 )
Let λmax (τ0 ) denote max{d21 /(d21 + 1), . . . , d2t /(d2t + 1)}. Note that λmax (τ0 ) < 1. Then
−1



W W T W + A(τ )
Hence,

−1



v T W W T W + A(τ )

W T ≤ λmax (τ0 )In .

W T v ≤ λmax (τ0 )

n
X

vi2 .

i=1

So


T

E η M η|v, τ , y



2

≤ p + q + (1 + c )λmax (τ0 )

n
X

= A1 (c) + (1 + c2 )λmax (τ0 )

vi2

i=1
n
X

1
+ 1 + 2 lT M −1 l
c




vi2 ,

(24)

i=1

where A1 (c) ≡ (1 + c−2 )lT M −1 l + p + q.
Now, we consider the expectation related to the conditional distribution of v and τ
given η ′ ≡ (β ′T , u′T )T and y. From Roy and Hobert (2007), we have
E(vi2 |η ′ , y) =


(xT β′ +ziT u′ )φ(xTi β′ +ziT u′ )

 1 + (xTi β ′ + ziT u′ )2 + i
T ′
T ′

 1 + (xT β ′ + z T u′ )2 − (
i
i

if yi = 1

)

(

Φ xi β +zi u
T ′
T ′
′
T ′
xT
i β +zi u φ xi β +zi u
T
T
′
′
1−Φ xi β +zi u

) (

)

)

(

.

if yi = 0

The above expectation can be written as,


E(vi2 |η ′ , y) = 1 + wi∗T η ′

2

−





1 − Φ wi∗T η ′

where wi∗ = ci wiT is the ith row of W ∗ . Also,
−







wi∗T η ′ φ wi∗T η ′
1 − Φ wi∗T η ′





≤
≤


∗T ′
∗T ′

 (wi η )φ(wi η )
1−Φ(wi∗T η ′ )


0

sup

u∈(−∞,0]

17



wi∗T η ′ φ wi∗T η ′




,

if wi∗T η ′ ≤ 0
if wi∗T η ′ > 0

uφ (u)
≡ Λ,
1 − Φ (u)

(25)

where Λ ∈ (0, ∞). So E(vi2 |η ′ , y) ≤ 1 + Λ + (xTi β ′ + ziT u′ )2 for i = 1, . . . , n. Since
λmax (τ0 ) < 1, from (24)


E η T M η|η ′ , y



≤ A1 (c) + n(1 + c2 )λmax (τ0 )(1 + Λ)
+ (1 + c2 )λmax (τ0 )

n
X

(xTi β ′ + ziT u′ )2

i=1

′

< L + ρ(c)V (η ),

2
2
where L = A1 (c)+n(1+c
p )λmax (τ0 )(1+Λ) and ρ(c) = λmax (τ0 )(1+c ). Since λmax (τ0 ) <
−1
1, there exists c ∈ (0, λmax (τ0 ) − 1) such that ρ(c) < 1. That means there exist
ρ ∈ (0, 1) and L > 0 such that (23) holds.
Since M is positive definite, V (η) is unbounded off compact sets, that is, for any
a > 0, the set {η : V (η) ≤ a} is compact. We now show that η-chain is a Feller chain,
which means K (η, O) is a lower semi-continuous function on Rp+q for each fixed open
set O. That is,

lim inf K (ηm , O) = lim inf
m→∞

m→∞

= lim inf
m→∞

≥

Z

Z "Z

Rr+

O

Z

Z Z
O

k (η|ηm ) dη

O

Rr+

Rn

Z

Rn

#

π(η|v, τ , y)π(v, τ |ηm , y)dvdτ dη

π(η|v, τ , y) lim inf π(v, τ |ηm , y)dvdτ dη,
m→∞

where the inequality follows from Fatou’s lemma. Since we have π(v, τ |η, y) = π(v|η, y)π(τ |η, y).
And by the assumption, both π(v|η, y) and π(τ |η, y) are continuous functions in η. So
if η (m) → η,
lim inf K (ηm , O) ≥
m→∞

Z Z
O

Rr+

Z

Rn

π(η|v, τ , y)π(v, τ |η, y)dvdτ dη

= K (η, O) .
Thus by Meyn and Tweedie (1993)(chap. 15), (23) implies the block Gibbs sampler
Markov chain Ψ is geometrically ergodic.

A.2

Proof of Theorem 2

Before presenting the proof of Theorem 2, we list some preliminary results.
Lemma 1. The following inequalities hold for all τj ∈ R+ for j = 1, . . . , r:


1. S −1  Z T (I − PX ) Z
2.



−1

+

+



Pr



j=1 1/τj I − PZ T (I−PX )Z .

 (λp + τj ) Iqj , where λp is the largest eigenvalue of Z T (I − PX ) Z
Ri S −1 RiT
and Rj is a qj × q matrix with 0 and 1 such that Rj u = uj .
18

Let Z T (I − PX ) Z = U T ΛU be its spectral decomposition with λi ’s as diagonal elements


of Λ. Also let Z T (I − PX ) Z

+

= U T Λ+ U , and λ+
i = 1/λi if λi 6= 0, 0 otherwise.

The above results are similar to Lemma 1 in Román and Hobert (2012) and we omit
the proof.
Lemma 2. For any τ ∈ Rr+ , we have,
S −1 Z T (I − PX ) v ≤ K̂

n
X
i=1

|vi | ,

where K̂ is a finite number which depends on W .
Proof. Let ZP = Z T (I − PX ) and zPT i be the ith row of ZP . Then


ZPT ZP + D

S −1 Z T (I − PX ) v =
n 
X

=

ZPT ZP + D

i=1

=

n
X
i=1

n
X

zP k zPT k

−1

zP i vi ≤

+D

k=1

!−1

−1

ZPT v

n

X

ZPT ZP + D

i=1

zP i vi ≤

n
X
i=1

−1

zP i vi

|vi | Ki (τ ) ,

where


Ki2 (τ ) = zPT i zP i zPT i +


= zPT i zP i zPT i +


X

k∈{1,...,n}\{i}

k∈{1,...,n}\{i}

≤ sup tTi ti tTi +
ϑ∈Rn+q
+

X

X

−2

zP k zPT k + D 

zP i

zP k zPT k + D − Pr

1

j=1 1/τj

ϑk tk tTk +

n+q
X

k=n+1

k∈{1,...n}\{i}

≡ K̂i2 ,

Iq + Pr

1

j=1 1/τj

−2

ϑk tk tTk + ϑ1 Iq 

−2

Iq 

zP i

ti

where ti = zP i , tk = zP k for k ∈ {1, . . . , n}\ {i} and tk is q × 1 vector with 1 in (k − n)th
element, 0 otherwise for k = n + 1, . . . , n + q. By Lemma 3 in Román and Hobert (2012),
we know that K̂i2 is finite. Let K̂ = maxi K̂i , then
S −1 Z T (I − PX ) v ≤

n
X
i=1

19

|vi | K̂i ≤ K̂

n
X
i=1

|vi | .

The two block Gibbs sampler {η (m) , (v (m) , τ (m) )}∞
m=1 in Algorithm 3 has the same
rate of convergence as its two marginal chains η-chain and (v, τ )-chain. Here we work
with the η-chain, denoted as Ψ∗ = {η}∞
m=1 and establish its geometric rate of convergence. Since the conditional distribution of τj does not exist when ||uj || = 0 for
j ∈ A = {j ∈ {1, . . . , r} : bj = 0}, we define π(τ |η, y) as follows, such that η chain is
well defined on Rp+q ,
π (τ |η, y) =
n



Qr

j=1 fG


Qr



q
τj , 2j

+

uT uj
aj , j2

+ bj

j=1 fG (τj , 1, 1)



if η ∈
/N

,

if η ∈ N

o

where N = η ∈ Rp+q ; j∈A ||uj || = 0 , and fG (x, a, b) = ba /Γ(a)xa−1 e−bx . Note that
N being a zero probability measure set, π(τ |η, y) can be defined arbitrarily on N and
simulation of the Gibbs sampler is not affected by the definition of π(τ |η, y) on N .
The Markov chain defined on Rp+q \N is denoted as Ψ̃∗ , which is Harris ergodic on
p+q
R \N . The proof of geometric ergodicity of Ψ∗ is through that of Ψ̃∗ . The following
proof establishes the geometric ergodicity of Ψ̃∗ .
Q

Proof. The proof of the geometric ergodicity is by constructing a drift function which
has the following form,
V (η) =

n 
X

xTi β

+

ziT u

i=1

2

+

r 
X

uTj uj

j=1

−c

,

(26)

where c ∈ (0, 1/2) is a positive constant. We show that for any η, η ′ ∈ Rp+q \N , there
exists ρ1 ∈ [0, 1) and L1 > 0 such that
E[V (η) |η ′ ] ≤ ρ1 V (η ′ ) + L1 .

(27)

The expectation can be calculated using two steps by Fubini’s theorem. First
E[V (η)|v, τ , y] is calculated with respect to the conditional distribution of η given v, τ
and y.
We have W Σ∗−1 W T = PX +(I − PX ) ZS −1 Z T (I − PX ) and (I − PX ) = (I − PX )2 .
Let P̃ = (I − PX ) ZD −1/2 , then
(I − PX ) ZS −1 Z T (I − PX ) = (I − PX )2 ZS −1 Z T (I − PX )2


= (I − PX ) P̃ P̃ T P̃ + I
Thus, W Σ∗−1 W T ≤ PX + I − PX = I and
E

" n
X
i=1

xTi β

+

ziT u

2

#

h

−1

P̃ T (I − PX ) ≤ I − PX .

|v, τ | ≤ E η T Σ∗ η|τ , v

= p + q + vW Σ∗−1 W T v ≤ p + q + v T v.
20

i

According to Román and Hobert (2012) for c ∈ (0, 1/2)
E



uTj uj

−c



|v, τ ≤ 2−c

i
Γ (qj /2 − c) h c
λp + τjc .
Γ (qj /2)

where λp is the largest eigenvalue of Z T (I − PX )Z. Then E[V (η)|v, τ , y] becomes,
E [V (η) |τ , v, y] ≤ v T v + 2−c

r
X
Γ (qj /2 − c)

Γ (qj /2)

j=1

τjc + 2−c

r
X
Γ (qj /2 − c)

Γ (qj /2)

j=1

λcp + p + q.

Now we consider the expectation corresponding to the conditional distribution of v
and τ given η ′ and y. The following procedure is almost the same as in Roy and Hobert
(2007), except wehave η ′ instead
of β.

Pn
2 |η ′ = n. If η ′ ∈ Rp+q \ {0}, denote A , . . . , A n as all the
If η ′ = 0, E
v
1
2
i=1 i
subsets of Nn = {1, 2 . . . , n} and
o

n

Sj = η ′ ∈ Rp+q \ {0} : wiT η ′ ≤ 0 for all i ∈ Aj and wiT η ′ > 0 for all i ∈ Āj ,
where Āj is the complement of Aj . Here Sj has the same properties as in Roy and Hobert
(2007). For j ∈ C = {j ∈ N2n : Sj 6= ∅}, define
′

P

i∈Aj

Rj η = Pn

i=1



wi∗T η ′

wi∗T η ′

2



P

i∈Aj
2
wi∗T η ′

2 = P

i∈Aj

wi∗T η ′

+

P

2

i∈A¯j

wi∗T η ′

2 .

where wi∗ ’s are defined in the proof of Theorem 1. For η ∈ Sj , j ∈ C, we have
E

" n
X
i=1

#

vi2 |η ′ = n +
=n+

n 
X

wi∗T η ′

i=1

n 
X

wi∗T η ′

i=1

≤n+

n 
X

wi∗T η ′

i=1

2

−



i∈Aj

2

+

2

+ nΛ −

= n (1 + Λ) + Rj η ′

X

i∈Aj



n 
X





wi∗T η ′ φ wi∗T η ′

X

1 − Φ wi∗T η ′











wi∗T η ′ φ wi∗T η ′
1 − Φ wi∗T η ′

X

wi∗T η ′

i∈A¯j

wi∗T η ′

i=1

2

2





−

X

i∈A¯j

−



X

i∈A¯j





wi∗T η ′ φ wi∗T η ′
1 − Φ wi∗T η ′










wi∗T η ′ φ wi∗T η ′
1 − Φ wi∗T η ′



,

where the inequality is based on the fact that uφ (u) / [1 − Φ (u)] ≥ u2 for u ≥ 0. Define
λj = supη′ ∈Sj {Rj (η ′ )} ∈ [0, 1] and
λ0 = max λj .
j∈C

Since conditions A1 and A2 are in force, using the techniques as in Roy and Hobert
(2007), it can be shown that λ0 < 1 and for all η ′ ∈ Rp+q ,
21



E

" n
X
i=1

vi2 |η ′

#

≤ λ0

n 
X

xTi β ′ + ziT u′

i=1

2

+ n (1 + Λ) .

Define Gj (−c) = 2c Γ(qj /2 + aj + c)/Γ(qj /2 + aj ) for j = 1, . . . , r, then
E

h

τjc |u′j

i

′
u′T
j uj
= 2 Gj (−c) bj +
2

"

−c

#−c





′
≤ Gj (−c) (2bj )−c I(0,∞) (bj ) + u′T
j uj

−c

We consider two cases, A is empty and A is not empty.
Case 1: A is not empty.
Then
E V (η) |η ′ ≤ λ0




n 
X

xTi β ′ + ziT u′

i=1

2

+ δ1 (c)

X

j∈A

′
u′T
j uj

−c

+ L (c) ,

where
δ1 (c) = 2−c max Gj (−c)
j∈A

Γ (qj /2 − c)
,
Γ (qj /2)

L (c) = n(1 + Λ) + p + q + 2−c λcp

r
X
Γ (qj /2 − c)

j=1

Γ (qj /2)

+ 2−c

X

j ∈A
/

Gj (−c)

Γ (qj /2 − c)
(2bj )−c .
Γ (qj /2)

By Román and Hobert (2012), there exists c ∈ C1 = (0, 1/2) ∩ (0, − max j∈A aj ) such
that δ1 (c) < 1. Thus, taking ρ1 = max(λ0 , δ1 (c)), and L1 = L(c), we have
E V (η) |η ′ ≤ ρ1 V η ′ + L1 .






Case 2: A is empty.
In this case, the conditional expectation of τjc can be bounded by a constant, that is,
E

h

τjc |u′j

i

′
u′T
j uj
= 2 Gj (−c) bj +
2

"

−c

#−c

≤ Gj (−c) (2bj )−c .

Then
E V (η) |η ′ ≤ λ0




n 
X

xTi β ′ + ziT u′

i=1

2

+ L(c) ≤ λ0 V η ′ + L(c).


Thus in both cases, (27) holds. Since by condition A3, on Rp+q \N , τj has full
conditional distribution of a Gamma distribution for any given value of uj , the η-chain
is a Feller chain on Rp+q \N . By condition A1, V (η) is unbounded off compact on
Rp+q \N . Thus η-chain is geometrically ergodic on Rp+q \N .
Then we need to show that the original Markov chain Ψ∗ is geometrically ergodic.
The technique of Lemma 12 in Román (2012) can be applied here for this purpose.

22



I{0} (bj ) .

Let M and M̃ be the corresponding Mtfs of Ψ∗ and Ψ̃∗ , M n and M̃ n be the n-step
Mtfs, and X ≡ Rp+q , X̃ ≡ Rp+q \N . Since the Lebesgue measure of N is 0, then for any
x ∈ X̃ and B ∈ BX̃ = {X̃ ∩ A : A ∈ B}
M̃ (x, B) = M (x, B).
Let µ an µ̃ be the Lebesgue measures on X and X̃. Then Ψ∗ and Ψ̃∗ are µ-irreducible
and µ̃-irreducible respectively. Also, µ and µ̃ are the corresponding maximal irreducibility measures. These two Markov chains are also aperiodic. According to Theorem 15.0.1
in Meyn and Tweedie (1993), there exists a ν-petite set C ∈ BX̃ and a number M̃ ∞ (C)
such that µ̃(C) > 0 and
|M̃ n (x, C) − M̃ ∞ (C)| < MC ρnC ,
for all x ∈ C. The set C is a ν-petite set for X̃ means that ν is a nontrivial measure on
BX̃ such that,
∞
X

n=0

M̃ n (x, B)ã(n) ≥ ν(B)

for all x ∈ C and B ∈ BX̃ , where ã(n) is a mass function on {0, 1, 2, . . . , }.
Since M̃ n (x, B) = M n (x, B) for any x ∈ X̃ and B ∈ BX̃ , then M n (x, C) = M̃ n (x, C).
So
|M n (x, C) − M̃ ∞ (C)| < MC ρnC .
Also, since µ(N ) = 0, we know that µ(C) > 0. It can be checked that C is also petite
for the original Markov chain Ψ∗ . Thus from Theorem 15.0.1 of Meyn and Tweedie
(1993), it follows that Ψ∗ is geometrically erdogic.

A.3

Proof of Theorem 3

Proof. The integral in (20) can be written as,
Z

Rp+q−r

r
Y

j=1

qj
2

τj

Z

Rr+

Z

n
Y

Rr i=1





Φ w̃iT η̃


 yi h



1 − Φ w̃iT η̃

i1−yi

(28)


qj −1
X
τj
+aj −1
(djk + uj1 )2 + 2bj  du11 · · · dur1 dτ dη̃
exp − u2j1 +
2
k=1

Pqj −1
Let d¯j = k=1
djk , (28) becomes,

23

Z

Rp+q−r

· (2π)
=

r
2

Rr+ i=1

r
Y

j=1
n
Y

Z

Rp+q−r i=1

· (2π)
≤K1

r
2

n
Y

Z

r
Y

−1/2

qj

qj
2

τj



Φ w̃iT η̃





1 − Φ w̃iT η̃

yi h



Φ w̃iT η̃



i1−yi



qj −1
2
qj − 1 ¯2
τj X 
djk − d¯j +
dj + 2bj  dτ dη̃
exp − 
2 k=1
qj



1 − Φ w̃iT η̃

j=1

Rp+q−r i=1

yi h

+aj − 32

−1/2
qj
Γ (qj /2

n
Y

Z



Φ w̃iT η̃

i1−yi

+ aj − 1/2) 2

 yi h



qj
2

1 − Φ w̃iT η̃

+aj − 12



qj −1 



djk − d¯j

X

k=1

r
i1−yi Y

j=1



2

qj −1 



qj − 1 ¯2
+
dj + 2bj 
qj

djk − d¯j

X

k=1

− qj −aj + 1

2

+

2

2

dη̃

− qj −aj + 1

qj − 1 ¯2
dj + 2bj 
qj

2

2

where K1 is a constant depending on r, qj and aj , j = 1, . . . , r.
Let δi , i = 1, . . . , n be i.i.d standard normal random variables and δ∗ = (c1 δ1 , . . . , cn δn )T .
h n

We have E 1 ci w̃iT η̃ ≤ ci δi
n
Y



Φ w̃iT η̃

i=1

oi

= Φ w̃iT η̃

iyi h



i1−yi

h 

 yi h

1 − Φ w̃iT η̃



1 − Φ w̃iT η̃

i1−yi

h n

. Then

= E 1 W̃ ∗ η̃ ≤ δ∗

oi

.

Since conditions B3 and B4 are in force, according to Chen and Shaon(2001) (Lemma
o
4.1), there exists a constant K depending on W̃ and y, such that 1 W̃ ∗ η̃ ≤ δ∗ ≤
1 {kη̃k ≤ K kδ∗ k}. Then (28) is bounded by


K1 E 


Z

Rp+q−r

1 {kη̃k ≤ K kδ∗ k}



≤ 2p K p K1 E kδ∗ kp

Z

≤ 2p K p K1 E kδ∗ kp

Z






r
Y

j=1




r
Y

|djk |≤Kkδ ∗ k j=1
|djk |

≤Kkδ ∗ k

r
Y

j=1

qj −1 

X

k=1



djk − d¯j

qj −1 



X

k=1

2

djk − d¯j

where dj = (dj1 , . . . , dj,qj −1 )T .
We consider two cases.
Case 1: aj < bj = 0.

24

qj − 1 ¯2
dj + 2bj 
qj

+

2

qj − 1 ¯2
dj + 2bj
qj

− qj −aj + 1

+

2

− qj −aj + 1

qj − 1 ¯2
dj + 2bj 
qj

!− qj −aj + 1
2

2

2

2



dd1 · · · , ddr  .


2



dη̃ 




dd1 · · · , ddr 


(29)

dη̃,

We know that,
Z

|dj1 |≤1

=

−

q

2 − 2j −aj + 12
qj −1

 X
ddj1 =
djk  


k=1

1
2 − qj − 2aj

1
2 − qj − 2aj

2  −
X


djk  



qj −1

qj
2

−aj + 12

k=1

1

 −1 +
2 − qj − 2aj

X

k=1

qj −1

X

k=2

2 

djk  

q
− 2j

−aj + 12



qj −1

X

−1 +

k=2



djk 



q
q
2 − 2j −aj +1 
2 − 2j −aj +1 




qj −1
qj −1


X
X
1




djk  
+ −1 +
. (30)
≤
djk  
 1+

2 − qj − 2aj 




k=2
k=2


Note that, (30) is a linear combination of terms α0 + dj,qj −1

−2aj

1

, which is finite

−1

since aj < 0, where α0 is a constant value of linear combination of -1 and 1. So,
qj − 1 ¯2
dj + bj
qj

Z

|djk |≤Kkδ ∗ k
qj

!− qj −aj + 1
2

2

1

= [qj (qj − 1)] 2 +aj − 2 (K kδ∗ k)−2aj

ddj

Z

|djk |≤1

≤K2j kδ∗ k−2aj .

q

2 − 2j −aj + 12
qj −1

 X
ddj
djk  


k=1

Case 2: bj > 0.
We have
Z

|djk |≤Kkδ ∗ k
q
− 2j

≤(2bj )

−aj + 12

qj − 1 ¯2
dj + 2bj
qj

!− qj −aj + 1
2

2

ddj ≤

Z

qj

|djk |≤Kkδ ∗ k

2qj −1 K qj −1 kδ∗ kqj −1 ≤ K3j kδ∗ kqj −1 .

So (28) can be bounded by


2p K p K1 E kδ∗ kp


≤2p K p K1 

r
Y

j=1

r n
Y

j=1

o



K2j kδ∗ k−2aj + K3j kδ∗ kqj −1 

K2j E kδ∗ k

p−2

Pr

j=1

25

aj

+

r
Y

j=1

1

(2bj )− 2 −aj + 2 ddj



K3j E kδ∗ kp+q−r  < ∞.



djk 

k=2

k=2



qj −1



q

2 − 2j −aj + 12 
qj −1
qj −1
X
X


1 +
djk 
djk  
 1+



1



−1

Remark 4. If qj = 1 and bj = 0, (28) is ∞ since
can be 1.

A.4

qj /2+aj −3/2
dτj
R τj

R

= ∞. If bj > 0, qj

Proof of Theorem 4

Proof. We define the same conditional distribution for τj if η ∈ N the same way as in
A.2.
Consider the following drift function on Rp+q \N ,
V (η) = α

n 
X

xTi β + ziT u

i=1

2

r
X

+



Gj (s) uTj uj

j=1

s

+

r 
X

j=1

uTj uj

−c

.

where α, s ∈ S1 = (0, 1] ∩ (0, s̃) and c ∈ C1 = (0, 1/2) ∩ (0, − max j∈A aj ) are positive
constants to be chosen later. We need to show that for any η, η ′ ∈ Rp+q \N , there exists
a constant ρ2 ∈ [0, 1) and L2 > 0 such that
E[V (η) |η ′ ] ≤ ρ2 V (η ′ ) + L2 .

(31)

Again, we calculate the expectation with respect to η conditional on v and τ first.
Same calculations as in the proof of Theorem 2 show that,
E

"

n 
X

xTi β + ziT u

i=1

2

#

|τ , v ≤ p + q + v T v.

For s ∈ (0, 1], by Jensen inequality,
E
Also,


h



uTj uj

s

i

h



|v, τ ≤ E uTj uj |v, τ

is

.





E uTj uj |v, τ = tr Rj S −1 Rj + [E (Rj u|v, τ )]T [E (Rj u|v, τ )] .
For the trace part, we have


tr Rj S

−1

RjT







T

= tr Rj Z (I − PX ) Z
= ξj + ςj

r
X

+

RjT





h



+ tr Rj I − PZ T (I−PX )Z RjT

r
iX

τl−1

l=1

τl−1 ,

l=1



where ξj = tr Rj



ZT

(I − PX ) Z

+

RjT



h





i

and ςj = tr Rj I − PZ T (I−PX )Z RjT . For

26

the second part,
[E (Rj u|v, τ )]T [E (Rj u|v, τ )] = v T (I − PX ) ZS −1 RjT Rj S −1 Z T (I − PX ) v
≤ v T (I − PX ) ZS −1 S −1 Z T (I − PX ) v

= S −1 Z T (I − PX ) v
≤

K̂

n
X
i=1

!2

2

2

≤ K̂ n

|vi |

n
X

vi2 .

i=1

The second inequality holds because of Lemma 2. Then
h



E uTj uj |v, τ

is



≤  ξ j + ςj

r
X

τj−1 + K̂ 2 n

j=1

n
X
i=1

s

vi2  ≤ ξjs + ςjs

r
X

τl−s + K̂ 2s ns

n
X

vi2s .

i=1

l=1

For vis , if vi2 ≤ 1, then vi2s ≤ 1, and if vi2s > 1, then vi2s < vi2 , so vi2s ≤ 1 + vi2 . Thus,
h



E uTj uj |v, τ

is

≤ ςjs

We also have,
E



uTj uj

−c

r
X

τl−s + K̂ 2s ns

n
X

vi2 + K̂ 2s n1+s + ξjs .

i=1

l=1



|v, τ ≤ 2−c

i
Γ (qj /2 − c) h c
λp + τjc .
Γ (qj /2)

Thus,
E [V (η) |v, τ ] ≤ (α + δ1 (s))

n
X

vi2 + δ2 (s)

i=1

r
X

τj−s + 2−c

r
X
Γ (qj /2 − c)

j=1

j=1

Γ (qj /2)

τjc + κ1 (α, s, c) ,

where
δ1 (s) = K̂ 2s ns

r
X

Gj (s) ,

j=1

δ2 (s) =

r
X

Gj (s) ςjs ,

j=1

κ1 (α, s, c) = α (p + q) +

r
X





Gj (s) K̂ 2s n1+s + ξjs + 2−c λcp

j=1

r
X
Γ (qj /2 − c)

j=1

Γ (qj /2)

.

Then we calculate the expectation of the drift function with respect to the distribution of v and τ conditional on η ′ .
P
When calculating the upper bound of E( ni=1 vi2 |η ′ ), we need to take into account
P
the fact that W is not full rank in this setting. But, E( ni=1 vi2 |η) can be written as,
E

" n
X
i=1

#

vi2 |η ′ = n +

n 
X

w̃iT η̃ ′

i=1

27

2


 

n
w̃iT η̃ ′ φ w̃iT η̃ ′
X
 .
−
T
i=1

1 − Φ w̃i η̃ ′

And we know that W̃ is a full rank matrix. Then the same technique as in the proof
of Theorem 2 can be applied. That is, there exists λ0 ∈ [0, 1) such that
E

" n
X

vi2 |η ′

i=1

#

≤ λ0
= λ0
= λ0

n 
X

i=1
n 
X

i=1
n 
X

w̃iT η̃ ′

2

+ n (1 + Λ)

wiT η ′

2

+ n (1 + Λ)

xTi β ′ + ziT u′

i=1

2

+ n (1 + Λ) .

For s ∈ S1 , we have
E

h

τj−s |η ′

i

′
u′T
j uj
= 2 Gj (s) bj +
2
s

!s



′
≤ Gj (s) u′T
j uj

s

+ 2s Gj (s) bsj ,

and c ∈ C1 ,
′
u′T
j uj
E τjc |η ′ = 2−c Gj (−c) bj +
2

h

"

i

#−c



−c

≤ Gj (−c) (2bj )

I(0,∞) (bj ) +



′
u′T
j uj

−c

I{0} (bj ) .

Again, we consider it in two cases.
Case 1: A is not empty.
′



E V (η) |η = αλ0

 n 
r
2

s
X
δ1 (s) X
′
xTi β ′ + ziT u′ + δ2 (s)
Gj (s) u′T
u
1+
j
j
α
i=1
j=1



+ δ3 (c)

X

′
u′T
j uj

j∈A

−c

+ L2 (α, s, c) ,

where
δ3 (c) = 2−c max Gj (−c)
j∈A

Γ (qj /2 − c)
,
Γ (qj /2)

L2 (α, s, c) = κ1 (α, s, c) + n (1 + Λ) (α + δ1 (s)) + δ2 (s) 2s

r
X

Gj (s) bsj

j=1

+2−c

X Γ (qj /2 − c)

j ∈A
/

Γ (qj /2)

Gj (−c) (2bj )−c .

And we know that for c ∈ C1 , δ3 (c) < 1. By condition 2 of Theorem 4, then
0 δ1 (s)
. So there exists a ρ2 (α, s, c) =
δ2 (s) < 1. For a fixed s, λ0 1 + δ1α(s) < 1 if α > λ1−λ
0
n





o

max λ0 1 + δ1α(s) , δ2 (s), δ3 (c) < 1 and L2 (α, s, c) > 0 such that (31) holds.
Case 2: A is empty

28





′

E V (η) |η = αλ0

 n 
r

2
s
X
δ1 (s) X
′
xTi β ′ + ziT u′ + δ2 (s)
Gj (s) u′T
1+
j uj
α
i=1
j=1



+ L2 (α, s, c) .
n





o

Similarly, So there exists ρ2 (α, s, c) = max λ0 1 + δ1α(s) , δ2 (s) < 1 and L2 (α, s, c) >
0 such that (31) holds.
Since η-chain is a Feller chain on Rp+q \N , and V (η) is unbounded off compact on
p+q
R \N , η-chain is geometrically ergodic on Rp+q \N . Using the same techniques as in
A.2, it can be shown that the original Markov chain defined on Rp+q is also geometrically
ergodic.

References
Albert, J. H. and Chib, S. (1993). Bayesian analysis of binary and polychotomous
response data. Journal of the American statistical Association, 88(422):669–679.
Asmussen, S. and Glynn, P. W. (2011). A new proof of convergence of MCMC via the
ergodic theorem. Statistics and Probability Letters, 81:1482–1485.
Baragatti, M. (2011). Bayesian variable selection for probit mixed models applied to
gene selection. Bayesian Analysis, 6(2):209–229.
Bernstein, D. S. (2009). Matrix mathematics: theory, facts, and formulas. Princeton
University Press.
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, M.-H. and Shao, Q.-M. (2001). Propriety of posterior distribution for dichotomous quantal response models. Proceedings of the American Mathematical Society,
129(1):293–302.
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. Sankhyā:
The Indian Journal of Statistics, Series A, pages 57–85.
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.
Doss, H. and Hobert, J. P. (2010). Estimation of Bayes factors in a class of hierarchical
random effects models using a geometrically ergodic MCMC algorithm. Journal of
Computational and Graphical Statistics, 19:295–312.

29

Flegal, J. M. and Jones, G. L. (2010). Batch means and spectral variance estimators in
Markov chain Monte Carlo. The Annals of Statistics, 38(2):1034–1070.
Gilks, W. R. and Wild, P. (1992). Adaptive rejection sampling for Gibbs sampling.
Applied Statistics, pages 337–348.
Hobert, J. P. and Marchev, D. (2008). A theoretical comparison of the data augmentation, marginal augmentation and PX-DA algorithms. The Annals of Statistics,
36(2):532–554.
Johnson, A. A. and Jones, G. L. (2010). Gibbs sampling for a bayesian hierarchical
general linear model. Electronic Journal of Statistics, 4:313–333.
Jones, G. L. and Hobert, J. P. (2001). Honest exploration of intractable probability
distributions via Markov chain Monte Carlo. Statistical Science, pages 312–334.
Jones, G. L. and Hobert, J. P. (2004). Sufficient burn-in for gibbs samplers for a hierarchical random effects model. Annals of statistics, 32(2):784–817.
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, pages 27–40.
Liu, J. S. and Wu, Y. N. (1999). Parameter expansion for data augmentation. Journal
of the American Statistical Association, 94(448):1264–1274.
McCulloch, C. E., Searle, S. R., and Neuhaus, J. M. (2011). Generalized, Linear, and
Mixed Models. John Wiley & Sons.
Meng, X.-L. and van Dyk, D. A. (1999). Seeking efficient data augmentation schemes
via conditional and marginal augmentation. Biometrika, 86:301–320.
Meyn, S. P. and Tweedie, R. L. (1993). Markov chains and stochastic stability. Springer.
Polson, N. G., Scott, J. G., and Windle, J. (2013). Bayesian inference for logistic models
using Pólya–Gamma latent variables. Journal of the American statistical Association,
108(504):1339–1349.
R Development Core Team (2011). R: A Language and Environment for Statistical
Computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-90005107-0.
Roberts, G. O. and Rosenthal, J. S. (2001). Markov chains and de-initializing processes.
Scandinavian Journal of Statistics, 28:489–504.
Roberts, G. O., Rosenthal, J. S., et al. (1997). Geometric ergodicity and hybrid Markov
chains. Electron. Comm. Probab, 2(2):13–25.

30

Román, J. C. (2012). Convergence analysis of block Gibbs samplers for Bayesian general
linear mixed models. PhD thesis, University of Florida.
Romá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(6):2823–2849.
Román, J. C. and Hobert, J. P. (2015). Geometric ergodicity of Gibbs samplers for
Bayesian general linear mixed models with proper priors. Linear Algebra and its
Applications, 473:54–77.
Roy, V. (2012a). Convergence rates for MCMC algorithms for a robust Bayesian binary
regression model. Electronic Journal of Statistics, 6:2463–2485.
Roy, V. (2012b). Spectral analytic comparisons for data augmentation. Stat. and Prob.
Letters, 82:103–108.
Roy, V. (2014). Efficient estimation of the link function parameter in a robust Bayesian
binary regression model. Computational Statistics & Data Analysis, 73:87–102.
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(4):607–623.
Tan, A. and Hobert, J. P. (2009). Block Gibbs sampling for Bayesian random effects
models with improper priors: Convergence and regeneration. Journal of Computational and Graphical Statistics, 18(4):861–878.
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.
Venables, W. N. and Ripley, B. D. (2002). Modern applied statistics with S. Springer.

31

