Biostatistics (2007), 8, 1, pp. 1–27
doi:10.1093/biostatistics/kxl000
Advance Access publication on May 15, 2006

Bayesian Inference for Generalized Linear Mixed Models
YOUYI FONG
Department of Biostatistics, University of Washington, Seattle, WA 98112, USA
HÅVARD RUE
The Norwegian University for Science and Technology, Trondheim, Norway
JON WAKEFIELD∗
Departments of Statistics and Biostatistics, University of Washington, Seattle, WA 98112, USA
jonno@u.washington.edu

S UMMARY
Generalized linear mixed models (GLMMs) continue to grow in popularity due to their ability to directly
acknowledge multiple levels of dependency, and model different data types. For small sample sizes especially, likelihood-based inference can be unreliable with variance components being particularly difficult
to estimate. A Bayesian approach is appealing, but has been hampered by the lack of a fast implementation, and the difficulty in specifying prior distributions with variance components again being particularly
problematic. Here we briefly review previous approaches to computation in Bayesian implementationsof GLMMs, and illustrate in detail the use of integrated nested Laplace approximations in this context.
We consider a number of examples, carefully specifying prior distributions on meaningful quantities ineach case. The examples cover a wide range of data types including those requiring smoothing over time,and a relatively complicated spline model for which we examine our prior specification in terms of the∗ To whom correspondence should be addressed.

c The Author 2006. Published by Oxford University Press. All rights reserved. For permissions, please e-mail: journals.permissions@oxfordjournals.org.

2

Y. F ONG , H. RUE AND J. WAKEFIELD

implied degrees of freedom. We conclude that Bayesian inference is now practically feasible for GLMMs,
and provides an attractive alternative to likelihood-based approaches such as penalized quasi-likelihood.
As with likelihood-based approaches, great care is required in the analysis of clustered binary data, sinceapproximation strategies may be less accurate for such data.

Keywords: Integrated nested Laplace approximations; Longitudinal data; Penalized quasi-likelihood; Prior specification; Spline models

1. I NTRODUCTION
Generalized linear mixed models (GLMMs) combine a generalized linear model with normal random effects on the linear predictor scale, to give a rich family of models that have been used in a wide variety
of applications, see for example Diggle et al. (2002), Verbeke and Molenberghs (2000, 2005), and McCulloch et al. (2008). This flexibility comes at a price, however, in terms of analytical tractability, which
has a number of implications including computational complexity, and an unknown degree to which inference is dependent on modeling assumptions. Likelihood-based inference may be carried out relatively
easily within many software platforms (except perhaps for binary responses), but inference is dependent
on asymptotic sampling distributions of estimators, with few guidelines available as to when such theory
will produce accurate inference. A Bayesian approach is attractive, but requires the specification of priordistributions which is not straightforward, in particular for variance components. Computation is also an
issue since the usual implementation is via Markov chain Monte Carlo (MCMC), which carries a large
computational overhead. The seminal article of Breslow and Clayton (1993) helped to popularize GLMMs
and placed an emphasis on likelihood-based inference via penalized quasi-likelihood (PQL). It is the aimof this article to describe, through a series of examples (including all of those considered in Breslow andBayesian Generalized Linear Mixed Modeling

3

Clayton, 1993), how Bayesian inference may be performed with computation via a fast implementation,
and with guidance on prior specification.
The structure of this article is as follows. In Section 2 we define notation for the GLMM, and in Section
3 describe the integrated nested Laplace approximation (INLA) which has recently been proposed as acomputationally convenient alternative to MCMC. Section 4 gives a number of prescriptions for priorspecification. Three examples are considered in Section 5 (with additional examples being reported in the
supplementary material, along with a simulation study that reports the performance of INLA in the binary
response situation). We conclude the paper with a discussion in Section 6.

2. T HE G ENERALIZED L INEAR M IXED M ODEL
GLMMs extend the generalized linear model, as proposed by Nelder and Wedderburn (1972) and comprehensively described in McCullagh and Nelder (1989), by adding normally distributed random effects
on the linear predictor scale. Suppose Yij is of exponential family form: Yij |θij , φ1 ∼ p(·) where p(·) isa member of the exponential family, that isp(yij |θij , φ1 ) = exp
yij θij − b(θij )+ c(yij , φ1 ) ,a(φ1 )

for i = 1, ..., m units (clusters), and j = 1, ..., ni , measurements per unit, and where θij is the (scalar)canonical parameter. Let µij = E[Yij |β, bi , φ1 ] = b0 (θij ) withg(µij ) = ηij = xij β + z ij bi ,where g(·) is a monotonic link function, xij is 1 × p and z ij is 1 × q, with β a p × 1 vector of fixed effects
and bi a q × 1 vector of random effects, hence θij = θij (β, bi ). Assume bi |Q ∼ N (0, Q−1 ), where
the precision matrix Q = Q(φ2 ) depends on parameters φ2 . For some choices of model the matrix Q is

4

Y. F ONG , H. RUE AND J. WAKEFIELD

singular; examples include random walk models (as considered in Section 5.2) and intrinsic conditional
autoregressive (ICAR) models. We further assume that β is assigned a normal prior distribution. Let
γ = (β, b) denote the G × 1 vector of parameters assigned Gaussian priors. We also require priors forφ1 (if not a constant) and for φ2 . Let φ = (φ1 , φ2 ) be the variance components for which non-Gaussianpriors are assigned, with V = dim(φ).

3. I NTEGRATED N ESTED L APLACE A PPROXIMATION
Before the MCMC revolution, there were few examples of the applications of Bayesian GLMMs since,outside of the linear mixed model, the models are analytically intractable. Kass and Steffey (1989) describe the use of Laplace approximations in Bayesian hierarchical models, while Skene and Wakefield(1990) used numerical integration in the context of a binary GLMM. The use of MCMC for GLMMs
is particularly appealing since the conditional independencies of the model may be exploited when therequired conditional distributions are calculated. Zeger and Karim (1991) described approximate Gibbs
sampling for GLMMs, with non-standard conditional distributions being approximated by normal distributions. More general Metropolis-Hastings algorithms are straightforward to construct, see for example
Clayton (1996) and Gamerman (1997). The WinBUGs (Spiegelhalter et al., 1998) software example manuals contain many GLMM examples. There are now a variety of additional software platforms for fitting
GLMMs via MCMC including JAGS (Plummer, 2009) and BayesX (Fahrmeir et al., 2004). A large
practical impediment to data analysis using MCMC is the large computational burden. For this reason we
now briefly review the INLA computational approach upon which we concentrate. The method combines
Laplace approximations and numerical integration in a very efficient manner, see Rue et al. (2009) for aBayesian Generalized Linear Mixed Modeling5

more extensive treatment. For the GLMM described in Section 2, the posterior is given byπ(γ, φ|y) ∝ π(γ|φ)π(φ)

mYp(y i | γ, φ)

i=1
1/2

∝ π(φ)π(β) | Q(φ2 ) |

(

)
m
X
1 T
exp − b Q(φ2 )b +log p(y i | γ, φ1 )2i=1where y i = (yi1 , ..., yini ) is the vector of observations on unit/cluster i. We wish to obtain the posterior
marginals π(γg | y), g = 1, ..., G, and π(φv |y), v = 1, ..., V . The number of variance components, V ,should not be too large for accurate inference (since these components are integrated out via Cartesian
product numerical integration, which does not scale well with dimension). We writeπ(γg |y) =Zπ(γg |φ, y) × π(φ|y)dφwhich may be evaluated via the approximationπe(γg |y) =≈Z

πe(γg |φ, y) × πe(φ|y)dφ

K
X

k=1

π
e(γg | φk , y) × πe(φk | y) × ∆k

(3.1)

where Laplace (or other related analytical approximations) are applied to carry out the integrations required for evaluation of πe(γg |φ, y). To produce the grid of points {φk , k = 1, ..., K} over which numerical integration is performed, the mode of πe(φ|y) is located, and the Hessian is approximated, from which

the grid is created and exploited in (3.1). The output of INLA consists of posterior marginal distributions,which can be summarized via means, variances and quantiles. Importantly for model comparison, the normalizing constant p(y) is calculated. The evaluation of this quantity is not straightforward using MCMC
(DiCiccio et al., 1997; Meng and Wong, 1996). The deviance information criterion (Spiegelhalter et al.,1998) is popular as a model selection tool, but in random effects models the implicit approximation in its

6

Y. F ONG , H. RUE AND J. WAKEFIELD

use is valid only when the effective number of parameters is much smaller than the number of independent
observations, see Plummer (2008).

4. P RIOR D ISTRIBUTIONS
4.1 Fixed Effects
Recall that we assume β is normally distributed. Often there will be sufficient information in the data for β
to be well estimated with a normal prior with a large variance (of course there will be circumstances underwhich we would like to specify more informative priors, for example, when there are many correlated
covariates). The use of an improper prior for β will often lead to a proper posterior though care should
be taken. For example, Wakefield (2007) shows that a Poisson likelihood with a linear link can lead to an
improper posterior if an improper prior is used. Hobert and Casella (1996) discuss the use of improper
priors in linear mixed effects models.
If we wish to use informative priors we may specify independent normal priors with the parameters
for each component being obtained via specification of two quantiles with associated probabilities. For
logistic and log-linear models these quantiles may be given on the exponentiated scale since these are moreinterpretable (as the odds ratio and rate ratio, respectively). If θ1 , θ2 are the quantiles on the exponentiatedscale, and p1 , p2 are the associated probabilities, then the parameters of the normal prior are given by:µ=z2 log(θ1 ) − z1 log(θ2 )z 2 − z1σ=log(θ2 ) − log(θ1 )z 2 − z1

where z1 , z2 are the p1 , p2 quantiles of a standard normal random variable. For example, in an epidemiological context we may wish to specify a prior on a relative risk parameter, exp(β1 ), which has a medianBayesian Generalized Linear Mixed Modeling7of 1 and a 95% point of 3 (if we think it is unlikely that the relative risk associated with a unit increase in
exposure exceeds 3). These specifications lead to β1 ∼ N (0, 0.6682 ).

4.2

Variance Components

We begin by describing an approach for choosing a prior for a single random effect, based on Wakefield(2009). The basic idea is to specify a range for the more interpretable marginal distribution of bi , anduse this to drive specification of prior parameters. We state a trivial lemma upon which prior specification
is based, but first define some notation. We write τ ∼ Ga(a1 , a2 ) for the gamma distribution with unnormalized density τ a1 −1 exp(−a2 τ ). For q-dimensional x we write x ∼ Tq (µ, Ω, d) for the Student’s
t distribution with unnormalized density [1 + (x − µ)T Ω−1 (x − µ)/d]−(d+q)/2 . This distribution has
location µ, scale matrix Ω and degrees of freedom d.
Lemma 1: Let b|τ ∼ N (0, τ −1 ) and τ ∼ Ga(a1 , a2 ). Integration over τ gives the marginal distributionof b as T1 (0, a2 /a1 , 2a1 ).
To decide upon a prior we give a range for a generic random effect b, and specify the degrees of freedom, d, and then solve for a1 , a2 . For the range (−R, R) we use the relationship ±td1−(1−q)/2pa2 /a1 =±R, where tdq is the 100 × q-th quantile of a Student t random variable with d degrees of freedom, to givea1 = d/2, a2 = R2 d/2(td1−(1−q)/2 )2 . In the linear mixed effects model, b is directly interpretable, while
for binomial or Poisson models it is more appropriate to think in terms of the marginal distribution ofexp(b), the residual odds and rate ratio, respectively, and this distribution is log Student’s t. For example,
if we choose d = 1 (to give a Cauchy marginal), and a 95% range of [0.1, 10] we take R = log 10 andobtain a = 0.5, b = 0.0164.
Another convenient choice is d = 2 to give the exponential distribution Exp(a2 ) (where the mean is

8

Y. F ONG , H. RUE AND J. WAKEFIELD−2
a−1. This leads to closed form expressions for the more interpretable quantiles of σ so that, for2 ) for σ2
example, if we specify the median for σ as σm we obtain a2 = σmlog 2.

Unfortunately, the use of Ga(, ) priors has become popular as a prior for σ −2 in a GLMM context,
arising from their use in the WinBUGS examples manual. As has been pointed out many times (e.g. Kelsall
and Wakefield, 1999, Gelman 2006, Crainiceanu et al. 2008) this choice places the majority of the priormass away from zero, and leads to a marginal prior for the random effects which is Student’s t with 2
degrees of freedom (so that the tails are much heavier than even a Cauchy), and difficult to justify in any
practical setting.
We now specify another trivial lemma, but first establish notation for the Wishart distribution. For theq × q non-singular matrix z we write z ∼ Wishartq (r, S) for the Wishart distribution with unnormalized
density |z|(r−q−1)/2 exp − 21 tr(zS −1 ) . This distribution has E[z] = rS and E[z −1 ] = S −1 /(r−q−1),

and we require r > q − 1 for a proper distribution.

Lemma: Let b = (b1 , ..., bq ), with b|Q ∼iid Nq (0, Q−1 ), Q ∼ Wishartq (r, S). Integration over Q gives
the marginal distribution of b as Tq (0, [(r − q + 1)S]−1 , r − q + 1).
The margins of a multivariate Student’s t are t also, which allows r and S to be chosen as in theunivariate case. Specifically, the k-th element of a generic random effect, bk , follows a univariate Student
t distribution with location 0, scale S kk /(r − q + 1), and degrees of freedom d = r − q + 1 where S kk iselement (k, k) of the inverse of S. We obtain r = d + q − 1 and S kk = (td1−(1−q)/2 )2 /(dR2 ). If a priori
we have no reason to believe that elements of b are correlated we may specify Sjk = 0 for j 6= k andS kk = 1/Skk , to recover the univariate specification, recognizing that with q = 1 the univariate Wishart
has parameters a1 = r/2, a2 = 1/(2S). If we believe that elements of b are dependent then we may
specify the correlations and solve for the off-diagonal elements of S. To ensure propriety of the posterior,

Bayesian Generalized Linear Mixed Modeling9

proper priors are required for Σ; Zeger and Karim (1991) use an improper prior for Σ, so that the posterior
is improper also.

4.3

Effective Degrees of Freedom Variance Components Prior

In Section 5.3 we describe the GLMM representation of a spline model. A generic linear spline model isgiven byy i = xi β +KXzik bk + i

k=1

where xi is a p × 1 vector of covariates with p × 1 associated fixed effects β, zik denote the spline
basis, bk ∼iid N (0, σb2 ) and i ∼iid N (0, σ2 ), with bk and i independent. Specification of a priorfor σb2 is not straightforward, but may be of great importance since it contributes to determining theamount of smoothing that is applied. Ruppert et al. (2003, p. 177) raise concerns, “about the instability of
automatic smoothing parameter selection even for single predictor models”, and continue, “Although weare attracted by the automatic nature of the mixed model-REML approach to fitting additive models, we
discourage blind acceptance of whatever answer it provides and recommend looking at other amounts ofsmoothing”. While we would echo this general advice, we believe that a Bayesian mixed model approach,with carefully chosen priors, can increase the stability of the mixed model representation. There has beensome discussion of choice of prior for σb2 in a spline context (Crainiceanu et al., 2005, 2008). More general
discussion can be found in Natarajan and Kass (2000) and Gelman (2006).
In practice (e.g. Hastie and Tibshirani, 1990) smoothers are often applied with a fixed degrees of
freedom. We extend this rationale by examining the prior degrees of freedom that is implied by the choiceσb−2 ∼ Ga(a1 , a2 ). For the general linear mixed modely = xβ + zb + 

10

Y. F ONG , H. RUE AND J. WAKEFIELD

we haveb + zbb = xβyb = C(C T C + Λ)−1 C T ywhere C = [x|z] is n × (p + K) and

Λ=



0p×p
0K×p

0p×K
σ2 cov(b)−1see, for example, Ruppert et al. (2003, Section 8.3). The total degrees of freedom associated with themodel isdf = tr{(C T C + Λ)−1 C T C}
which may be decomposed into the degrees of freedom associated with β and b, and extends easily tosituations in which we have additional random effects, beyond those associated with the spline basis (suchan example is considered in Section 5.3). In each of these situations the degrees of freedom associatedwith the respective parameter is obtained by summing the appropriate diagonal elements of (C T C +Λ)−1 C T C. Specifically, if we have j = 1, ..., d sets of random effect parameters (there are d = 2 inthe model considered in Section 5.3) then let E j be the (p + K) × (p + K) diagonal matrix with onesin the diagonal positions corresponding to set j. Then the degrees of freedom associated with this set is
dfj = tr{E j (C T C + Λ)−1 C T C. Note that the effective degrees of freedom changes as a function of K,
as expected. To evaluate Λ, σ2 is required. If we specify a proper prior for σ2 then we may specify the joint
prior as π(σb2 , σ2 ) = π(σ2 )π(σb2 |σ2 ). Usually, however, we assume the improper prior π(σ2 ) ∝ 1/σ2 ,since the data provide sufficient information with respect to σ2 . Hence we have found the substitution ofan estimate for σ2 (for example, from the fitting of a spline model in a likelihood implementation) to be apractically reasonable strategy.
As a simple non-spline demonstration of the derived effective degrees of freedom, consider a one-way

Bayesian Generalized Linear Mixed Modeling

11

ANOVA model
Yij = β0 + bi + ij
with bi ∼iid N (0, σb2 ), ij ∼iid N (0, σ2 ) for i = 1, ..., m = 10 groups and j = 1, ..., n = 5 observations
per group. For illustration we assume σb−2 ∼ Ga(0.5, 0.005). Figure 1 displays the prior distribution
for σ, the implied prior distribution on the effective degrees of freedom, and the bivariate plot of thesequantities. For clarity of plotting we exclude a small number of points beyond σ > 2.5 (4% of points).
In panel (c) we have placed dashed horizontal lines at effective degrees of freedom equal to 1 (completesmoothing) and 10 (no smoothing). From panel (b) we conclude that here the prior choice favors quite
strong smoothing. This may be contrasted with the Gamma prior with parameters (0.001, 0.001), which,in this example, gives greater than 99% of the prior mass on an effective degrees of freedom greater than
9.9, again showing the inappropriateness of this prior.
It is appealing to extend the above argument to non-linear models, but unfortunately this is not straightforward. For a non-linear model the degrees of freedom may be approximated by
df = tr{(C T W C + Λ)−1 C T W C}
 2 
−1 dµiand h = g −1 denotes the inverse link function. Unfortunately thiswhere W = diag Vi
dhquantity depends on β and b, which means that in practice we would have to use prior estimates for allof the parameters, which may not be practically possible. Fitting the model using likelihood and thensubstituting in estimates for β and b seems philosophically dubious.

12

Y. F ONG , H. RUE AND J. WAKEFIELD

(b)

2000

Frequency

1500

1500

0

0

500

500

1000

1000

Frequency

2000

2500

2500

3000

3000

3500

(a)

0.0

0.5

1.0

1.5

2.0

2.5

2

σ

4

6

8

10

Effective Degrees of Freedom

(a)

(b)

6
4
2

Effective Degrees of Freedom

8

10

(c)

0.0

0.5

1.0

1.5

2.0

2.5

σ

(c)

Fig. 1. Gamma prior for σ −2 with parameters 0.5 and 0.005, (a) implied prior for σ, (b) implied prior for the effective
degrees of freedom, (c) effective degrees of freedom versus σ.

Bayesian Generalized Linear Mixed Modeling13

4.4 Random Walk Models
Conditionally represented smoothing models are popular for random effects in both temporal and spatial
applications, see for example Besag et al. (1995) and Rue and Held (2005). For illustration, consider
models of the form1p(u|σu2 ) = (2π)−(m−r)/2 |Q? |1/2 σu−(m−r) exp − 2 uT Qu
2σu

(4.1)

where u = (u1 , ..., um ) is the collection of random effects, Q is a (scaled) “precision” matrix of rankm−r, whose form is determined by the application at hand, and |Q? | is a generalized determinant which isthe product over the m−r non-zero eigenvalues of Q. Picking a prior for σu is not straightforward because
σu has an interpretation as the conditional standard deviation, where the elements that are conditioned
upon depends on the application. We may simulate realizations from (4.1) to examine candidate priordistributions. Due to the rank deficiency, (4.1) does not define a probability density, and so we cannot
directly simulate from this prior. However, Rue and Held (2005) give an algorithm for generating samples
from (4.1):
1. Simulate zj ∼ N (0, λ−1j ), for j = m − r + 1, ..., m, where λj are the eigenvalues of Q (there arem − r non-zero eigenvalues as Q has rank m − r).
2. Return u = zm−r+1 en−r+1 +z3 e3 +...+zn em = Ez where ej are the corresponding eigenvectorsof Q, E is the m × (m − r) matrix with these eigenvectors as columns, and z is the (m − r) × 1
vector containing zj , j = m − r + 1, ..., m.
The simulation algorithm is conditioned so that samples are zero in the null-space of Q; if u is a sample
and the null-space is spanned by v 1 and v 2 , then uT v 1 = uT v 2 = 0. For example, suppose Q1 = 0 sothat the null space is spanned by 1, and the rank deficiency is 1. Then Q is improper since the eigenvalue

14

Y. F ONG , H. RUE AND J. WAKEFIELD

corresponding to 1 is zero, and samples u produced by the algorithm are such that uT 1 = 0. In Section
5.2 we use this algorithm to evaluate different priors, via simulation. It is also useful to note that if we wish
to compute the marginal variances only, simulation is not required, as they are available as the diagonal
elements of the matrixP−1Tj λj ej ej .

5. E XAMPLES
Here we report three examples, with four others described in the supplementary material. Together thesecover all of the examples in Breslow and Clayton (1993), along with an additional spline example. In
the first example results using the INLA numerical/analytical approximation described in Section 3 werecompared with MCMC as implemented in the JAGS software (Plummer, 2009) and, found to be accurate.
For the models considered in the second and third examples the approximation was compared with the
MCMC implementation contained in the INLA software.

5.1

Longitudinal Data

We consider the much analyzed epilepsy dataset of Thall and Vail (1990). These data concern the numberof seizures, Yij for patient i on visit j, with Yij |β, bi ∼ind Poisson(µij ), i = 1, ..., 59, j = 1, ..., 4. We
concentrate on the three random effects models fitted by Breslow and Clayton (1993):

log µij = xij β + b1i(5.1)

log µij = xij β + b1i + b0ij

(5.2)

log µij = xij β + b1i + b2i Vj /10(5.3)

Bayesian Generalized Linear Mixed Modeling15

where xij is a 1 × 6 vector containing a 1 (representing the intercept), an indicator for baseline measurement, a treatment indicator, the baseline by treatment interaction, which is the parameter of interest, age,and either an indicator of the fourth visit (models (5.1) and (5.2) and denoted V4 ), or visit number coded−3, −1, +1, +3 (model (5.3), and denoted Vj /10); β is the associated fixed effect. All three models
include patient-specific random effects b1i ∼ N 0, σ12 , while in model (5.2) we introduce independent“measurement errors”, b0ij ∼ N (0, σ02 ). Model (5.3) includes random effects on the slope associated withvisit, b2i with
We assume Q ∼ Wishart(r, S) with S =

b1i
b2i



∼ N (0, Q−1 ).

(5.4)

 S11 S12 
S21 S22 . For prior specification we begin with the bivariate

model and assume that S is diagonal. We assume the upper 95% point of the priors for exp(b1i ) andexp(b2i ) are 5 and 4, respectively, and that the marginal distributions are t with 4 degrees of freedom.
Following the procedure outlined in Section 4.2 we obtain r = 5 and S = diag(0.439, 0.591). We takethe prior for σ1−2 in model (5.1) to be Ga(a1 , a2 ) with a1 = (r − 1)/2 = 2 and a2 = 1/2S11 = 1.140 (sothat this prior coincides with the marginal prior obtained from the bivariate specification). In model (5.2)
we assume b1i and b0ij are independent, and that σ0−2 follows the same prior as σ1−2 , i.e. Ga(2, 1.140).
We assume a flat prior on the intercept, and assume that the rate ratios, exp(βj ), j = 1, ..., 5, lie between
0.1 and 10 with probability 0.95 which gives, using the approach described in Section 4.1, a normal priorwith mean 0 and variance 1.172 .
Table 1 gives PQL and INLA summaries for models (5.1)–(5.3). There are some differences between
the PQL and Bayesian analyses, with slightly larger standard deviations under the latter, which probably
reflects that with m = 59 clusters a little accuracy is lost when using asymptotic inference. There aresome differences in the point estimates which is at least partly due to the non-flat priors used — the priors

16

Y. F ONG , H. RUE AND J. WAKEFIELD

have relatively large variances, but here the data are not so abundant so there is sensitivity to the prior.
Reassuringly under all three models inference for the baseline-treatment interaction of interest is virtually
identical, and suggests no significant treatment effect. We may compare models using log p(y): for thethree models we obtain values of -674.8, -638.9 and -665.5, so that the second model is strongly preferred.

Variable
Base
Trt
Base × Trt
Age
V4 or V/10
σ0
σ1
σ2

Model (5.1)
PQL
INLA
0.87 ± 0.14
0.88 ± 0.15
−0.91 ± 0.41 −0.94 ± 0.44
0.33 ± 0.21
0.34 ± 0.22
0.47 ± 0.36
0.47 ± 0.38
−0.16 ± 0.05 −0.16 ± 0.05——
0.53 ± 0.06
0.56 ± 0.08——

Model (5.2)
PQL
INLA
0.86 ± 0.13
0.88 ± 0.15
−0.93 ± 0.40 −0.96 ± 0.44
0.34 ± 0.21
0.35 ± 0.23
0.47 ± 0.35
0.48 ± 0.39
−0.10 ± 0.09 −0.10 ± 0.09
0.36 ± 0.04
0.41 ± 0.04
0.48 ± 0.06
0.53 ± 0.07——

Model (5.3)
PQL
INLA
0.87 ± 0.14
0.88 ± 0.14
−0.91 ± 0.41 −0.94 ± 0.44
0.33 ± 0.21
0.34 ± 0.22
0.46 ± 0.36
0.47 ± 0.38
−0.26 ± 0.16 −0.27 ± 0.16——
0.52 ± 0.06
0.56 ± 0.06
0.74 ± 0.16
0.70 ± 0.14

Table 1. PQL and INLA summaries for the epilepsy data.

5.2 Smoothing of Birth Cohort Effects in an Age-Cohort Model
We analyze data from Breslow and Day (1975) on breast cancer rates in Iceland. Let Yjk be the numberof breast cancer of cases in age group j (20–24, ..., 80–84) and birth cohort k (1840–1849, ...,1940–1949)
with j = 1, ..., J = 13 and k = 1, ..., K = 11. Following Breslow and Clayton (1993) we assume
Yjk |µjk ∼ind Poisson(µjk ) withlog µjk = log njk + βj + βk + vk + uk

(5.5)

and where njk is the person-years denominator, exp(βj ), j = 1, ..., J, represent fixed effects for age
relative risks, exp(β) is the relative risk associated with a one group increase in cohort group, vk ∼iid
N (0, σv2 ) represent unstructured random effects associated with cohort k, with smooth cohort terms uk
following a second-order random effects model with E[uk |{ui : i < k}] = 2uk−1 −uk−2 and Var(uk |{ui :Bayesian Generalized Linear Mixed Modeling17

i < k}) = σu2 . This latter model is to allow the rates to vary smoothly with cohort. An equivalent
representation of this model is, for 2 < k < K − 1,
1
E[uk |{ul : l 6= k}] = (4uk−1 + 4uk+1 − uk−2 − uk+2 )6
Var(uk |{ul : l 6= k}) =σu2
6The rank of Q in the (4.1) representation of this model is K − 2 reflecting that both the overall level andthe overall trend are aliased (hence the appearance of β in (5.5)).
The term exp(vk ) reflects the unstructured residual relative risk and, following the argument in Section
4.2, we specify that this quantity should lie in [0.5,2.0] with probability 0.95, with a marginal log Cauchy
distribution, to obtain the gamma prior σv−2 ∼ Ga(0.5, 0.00149).
The term exp(uk ) reflects the smooth component of the residual relative risk, and the specificationof a prior for the associated variance component σu2 is more difficult, given its conditional interpretation.
Using the algorithm described in Section 4.2 we examined simulations of u for different choices of gamma
hyperparameters, and decided on the choice σu−2 ∼ Ga(0.5, 0.001); Figure 2 shows 10 realizations fromthe prior. The rationale here is to examine realizations to see if they conform to our prior expectations
and in particular exhibit the required amount of smoothing. All but one of the realizations vary smoothly
across the 11 cohorts, as is desirable. Due to the tail of the gamma distribution, we will always have someextreme realizations.
The INLA results, summarized in graphical form, are presented in Figure 2(b), alongside likelihood
fits in which the birth cohort effect is incorporated as a linear term and as a factor. We see that thesmoothing model provides a smooth fit in birth cohort, as we would hope.

Y. F ONG , H. RUE AND J. WAKEFIELD

2.0
1.0
0.5
0.1

0.2

Residual Relative Risk

5.0

10.0

18

2

4

6

8

10

Cohort Number

5.0

10.0

(a)

●

+

+
●

1.0

+
●
●
+

+
●

●
+

+
●

0.5

Relative Risk

2.0

●
+
●
+

+
●

+

0.1

0.2

●

2

4

6

8

10

Cohort Number(b)Fig. 2. (a) Ten realizations (on the relative risk scale) from the random effects second-order random walk modelin which the prior on the random effects precision is Ga(0.5,0.001). (b) Summaries of fitted models: the solid line
corresponds to a loglinear model in birth cohort, the circles to birth cohort as a factor, and “+” to the Bayesiansmoothing model.

Bayesian Generalized Linear Mixed Modeling
5.3

19

B-Spline Nonparametric Regression

We demonstrate the use of INLA for nonparametric smoothing using O’Sullivan splines, which are basedon a B-spline basis. We illustrate using data from Bachrach et al. (1999) that concerns longitudinal measurements of spinal bone mineral density (SBMD) on 230 female subjects aged between 8 and 27, andof one of four ethnic groups: Asian, Black, Hispanic and White. Let yij denote the SBMD measure forsubject i at occasion j, for i = 1, ..., 230 and j = 1, ..., ni with ni being between 1 and 4. Figure 3 showsthese data, with the grey lines indicating measurements on the same woman.
We assume the model

Yij = xi β 1 + ageij β2 +KXzijk b1k + b2i + ij

k=1

where xi is a 1 × 4 vector containing an indicator for the ethnicity of individual i, with β 1 the associated
4 × 1 vector of fixed effects, zijk is the k-th basis associated with age, with associated parameter b1k ∼N (0, σ12 ), and b2i ∼ N (0, σ22 ) are woman-specific random effects, finally, ij ∼iid N (0, σ2 ). All random
terms are assumed independent. Note that the spline model is assumed common to all ethnic groups andall women, though it would be straightforward to allow a different spline for each ethnicity. Writing thismodel in the formy = xβ + z 1 b1 + z 2 b2 +  = Cγ + 
we use the method described in Section 4.3 to examine the effective number of parameters implied bydifferent priors σ1−2 ∼ Ga(a1 , a2 ), σ2−2 ∼ Ga(a3 , a4 ).
To fit the model we first use the R code provided in Wand and Ormerod (2008) to construct the basis
functions, which are then input to the INLA program. Running the REML version of the model we obtainσb = 0.033 which we use to evaluate the effective degrees of freedoms associated with priors for σ12 and20

Y. F ONG , H. RUE AND J. WAKEFIELD

10

Hispanic

15

20

25

White
1.4
1.2
1.0

Spinal BMD

0.8
0.6

Asian

Black

1.4
1.2
1.0
0.8
0.6
10

15

20

25

age

Fig. 3. Spinal bone mineral density versus age by ethnicity. Measurements on the same woman are joined with greylines. The solid curve corresponds to the fitted spline, and the dashed lines to the individual fits.

21

Bayesian Generalized Linear Mixed Modeling

10

0

5

Effective Degrees of Freedom

400
100

200

Frequency

300

400
300
200
0

100

Frequency

(c)

15

(b)

500

(a)

0.04

0.08

0

5

10

15

0.00

0.04

σ1

Effective Degrees of Freedom

σ1

(d)

(e)

(f)

0.08

0.2

0.4

0.6

0.8

σ2

200
150
100
0

0
0.0

50

150
100
50

Frequency

200

Effective Degrees of Freedom

250

100 150 200 250 300
050Frequency

300

0.00

100

150

200

Effective Degrees of Freedom

0.0

0.2

0.4

0.6

σ2

Fig. 4. Prior summaries: (a) σ1 , the standard deviation of the spline coefficients, (b) effective degrees of freedom
associated with the prior for the spline coefficients, (c) effective degrees of freedom versus σ1 , (d) σ2 , the standard
deviation of the between-individual random effects, (e) effective degrees of freedom associated with the individual
random effects, (f) effective degrees of freedom versus σ2 . The vertical dashed lines on panels (a), (b), (d) and (e)correspond to the posterior medians.

0.8

22

Y. F ONG , H. RUE AND J. WAKEFIELD

σ22 . We assume the usual improper prior, π(σ2 ) ∝ 1/σ2 for σ2 . After some experimentation we settled on
the prior σ1−2 ∼ Ga(0.5, 5 × 10−6 ). For σ22 we wished to have a 90% interval for b2i of ±0.3 which, with
1 degree of freedom for the marginal distribution, leads to σ2−2 ∼ Ga(0.5, 0.00113). Figure 4 shows thepriors for σ1 and σ2 , along with the implied effective degrees of freedom under the assumed priors. For
the spline component the 90% prior interval for the effective degrees of freedom is [2.4,10].
Table 2 compares estimates from REML and INLA implementations of the model, and we see close
correspondence between the two. Figure 4 also shows the posterior medians for σ1 , σ2 and the two effective degrees of freedom. For the spline and random effects these correspond to 8 and 214, respectively. The
latter figure shows that there is considerable variability between the 230 women here. This is confirmed
in Figure 3 where we observe large vertical differences between the profiles. This figure also shows thefitted spline, which appears to mimic the trend in the data well.
VariableIntercept
Black
Hispanic
White
Age
σ1
σ2
σ

REML
0.560 ± 0.029
0.106 ± 0.021
0.013 ± 0.022
0.026 ± 0.022
0.021 ± 0.002
0.018?
0.109?
0.033?

INLA
0.563 ± 0.031
0.106 ± 0.021
0.013 ± 0.022
0.026 ± 0.022
0.021 ± 0.002
0.024 ± 0.006
0.109 ± 0.006
0.033 ± 0.002

Table 2. REML and INLA summaries for spinal bone data. Intercept corresponds to Asian group. For the entries
marked with a ? standard errors were unavailable.

5.4

Timings

For the three models in the longitudinal data example, INLA takes 1 to 2 seconds to run, using a single
CPU. To get estimates with similar precision with MCMC, we ran JAGS for 100,000 iterations, which

Bayesian Generalized Linear Mixed Modeling23

took 4 to 6 minutes. For the model in the temporal smoothing example, INLA takes 45 seconds to run,using one CPU. Because part of the INLA procedure can be executed in a parallel manner, if there are twoCPUs available, as is the case with today’s prevalent INTEL Core 2 Duo processors, INLA only takes 27
seconds to run. It is not currently possible to implement this model in JAGS. We ran the MCMC utility
built into the INLA software for 3.6 million iterations, to obtain estimates of comparable accuracy, which
took 15 hours. For the model in the B-spline nonparametric regression example, INLA took 5 seconds torun, using a single CPU. We ran the MCMC utility built into the INLA software for 2.5 million iterations
to obtain estimates of comparable accuracy, the analysis taking 40 hours.

6. D ISCUSSION
In this paper we have demonstrated the use of the INLA computational method for GLMMs. We havefound that the approximation strategy employed by INLA is accurate in general, but less accurate forbinomial data with small denominators. The on-line supplementary material contains an extensive simulation study, replicating that presented in Breslow and Clayton (1993). There are some suggestions in thediscussion of Rue et al. (2009) on how to construct an improved Gaussian approximation that does not
use the mode and the curvature at the mode. It is likely that these suggestions will improve the results forbinomial data with small denominators. There is an urgent need for diagnosis tools to flag when INLA is
inaccurate. Conceptually, computation for non-linear mixed effects models (Davidian and Giltinan, 1995;
Pinheiro and Bates, 2000) can also be handled by INLA, but this capability is not currently available.
The website www.r-inla.org contains all of the data and R scripts to perform the analyses andsimulations reported in the paper. The latest release of software to implement INLA can also be foundat this site. Recently, Breslow (2005) revisited PQL and concluded that, “PQL still performs remarkably

24

Y. F ONG , H. RUE AND J. WAKEFIELD

well in comparison with more elaborate procedures in many practical situations”. We believe that INLA
provides an attractive alternative to PQL for GLMMs, and we hope that this paper stimulates the greateruse of Bayesian methods for this class.

ACKNOWLEDGMENTS
JW was supported by grant R01 CA095994 from the National Institutes of Health.

R EFERENCES
Bachrach, L.K., Hastie, T., Wang, M.C., Narasimhan, B., and Marcus, R. (1999). Bone mineral acquisition in
healthy Asian, Hispanic, Black and Caucasian youth. A longitudinal study. Journal of Clinical Endocrinology
and Metabolism 84, 4702–4712.
Besag, J., Green, P.J., Higdon, D., and Mengersen, K. (1995). Bayesian computation and stochastic systems (withdiscussion). Statistical Science 10, 3–66.
Breslow, N.E. (2005). Whither PQL? In D. Lin and P.J. Heagerty (Eds.), Proceedings of the Second Seattle Symposium, New York, pp. 1–22. Springer-Verlag.
Breslow, N.E. and Clayton, D.G. (1993). Approximate inference in generalized linear mixed models. Journal of theAmerican Statistical Association 88, 9–25.
Breslow, N.E. and Day, N.E. (1975). Indirect standardization and multiplicative models for rates, with reference tothe age adjustment of cancer incidence and relative frequency data. Journal of Chronic Diseases 28, 289–301.
Clayton, D.G. (1996). Generalized linear mixed models. In W.R. Gilks, S. Richardson, and D.J. Spiegelhalter (Eds.),
Markov Chain Monte Carlo in Practice, pp. 275–301. Chapman and Hall.
Crainiceanu, C.M., Diggle, P.J., and Rowlingson, B. (2008). Bayesian analysis for penalized spline regression using
WinBUGS. Journal of the American Statistical Association 102, 21–37.
Crainiceanu, C.M., Ruppert, D., and Wand, M.P. (2005). Bayesian analysis for penalized spline regression using

Bayesian Generalized Linear Mixed Modeling

25

WinBUGS. Journal of Statistical Software 14.
Davidian, M. and Giltinan, D.M. (1995). Nonlinear Models for Repeated Measurement Data. Boca Raton: CRC
Press.
DiCiccio, T.J., Kass, R.E., Raftery, A., and Wasserman, L. (1997). Computing Bayes Factors by Combining Simulation and Asymptotic Approximations. Journal of the American Statistical Association 92, 903–915.
Diggle, P., Heagerty, P., Liang, K.Y., and Zeger, S. (2002). Analysis of Longitudinal Data, Second Edition. Oxford
University Press.
Fahrmeir, L., Kneib, T., and Lang, S. (2004). Penalized structured additive regression for space-time data: A Bayesian
perspective. Statistica Sinica 14, 715–745.
Gamerman, D. (1997). Sampling from the posterior distribution in generalized linear mixed models. Statistics andComputing 7, 57–68.
Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian Analysis 1, 515–534.
Hastie, T.J. and Tibshirani, R.J. (1990). Generalized additive models. Chapman and Hall, London.
Hobert, J.P. and Casella, G. (1996). The effect of improper priors on Gibbs sampling in hierarchical linear mixed
models. Journal of the American Statistical Association 91, 1461–1473.
Kass, R.E. and Steffey, D. (1989). Approximate Bayesian inference in conditionally independent hierarchical models(parametric empirical Bayes models). Journal of the American Statistical Association 84, 717–726.
Kelsall, J.E. and Wakefield, J.C. (1999). Discussion of “Bayesian models for spatially correlated disease and exposuredata” by N. Best, L. Waller, A. Thomas, E. Conlon and R. Arnold. In J.M. Bernardo, J.O. Berger, A.P. Dawid,
and A.F.M. Smith (Eds.), Sixth Valencia international meeting on Bayesian statistics, London. Oxford University
Press.
McCullagh, P. and Nelder, J.A. (1989). Generalized Linear Models, Second Edition. London: Chapman and Hall.
McCulloch, C.E., Searle, S.R., and Neuhaus, J.M. (2008). Generalized, Linear, and Mixed Models (Second ed.). John
Wiley and Sons.

26

Y. F ONG , H. RUE AND J. WAKEFIELD

Meng, X. and Wong, W. (1996). Simulating ratios of normalizing constants via a simple identity. Statistical Sinica 6,
831–860.
Natarajan, R. and Kass, R.E. (2000). Reference Bayesian methods for generalized linear mixed models. Journal of
the American Statistical Association 95, 227–237.
Nelder, J. and Wedderburn, R. (1972). Generalized linear models. Journal of the Royal Statistical Society, Series
A 135, 370–384.
Pinheiro, J.C. and Bates, D.M. (2000). Mixed-Effects Models in S and S-plus. New York: Springer.
Plummer, M. (2008). Penalized loss functions for Bayesian model comparison. Biostatistics 9, 523–539.
Plummer, M. (2009). JAGS Version 1.0.3 Manual. Technical report.
Rue, H. and Held, L. (2005). Gaussian Markov Random Fields: Thoery and Application. Chapman and Hall/CRC.
Rue, H., Martino, S., and Chopin, N. (2009). Approximte Bayesian inference for latent Gaussian models using
integrated nested Laplace approximations (with discussion). Journal of the Royal Statistical Society, Series B 71,
319–392.
Ruppert, D.R., Wand, M.P., and Carroll, R.J. (2003). Semiparametric Regression. New York: Cambridge University
Press.
Skene, A.M. and Wakefield, J.C. (1990). Hierarchical models for multi-centre binary response studies. Statistics in
Medicine 9, 919–929.
Spiegelhalter, D., Best, N., Carlin, B., and van der Linde, A. (1998). Bayesian measures of model complexity and fit(with discussion). Journal of the Royal Statistical Society, Series B 64, 583–639.
Spiegelhalter, DJ., Thomas, A., and Best, NG. (1998). WinBUGS User Manual, version 1.1.1. Cambridge, UK.
Thall, P.F. and Vail, S.C. (1990). Some covariance models for longitudinal count data with overdispersion. Biometrics 46, 657–671.
Verbeke, G. and Molenberghs, G. (2000). Linear Mixed Models for Longitudinal Data. New York: Springer-Verlag.
Verbeke, G. and Molenberghs, G. (2005). Models for Discrete Longitudinal Data. New York: Springer-Verlag.

Bayesian Generalized Linear Mixed Modeling27

Wakefield, J.C. (2007). Disease mapping and spatial regression with count data. Biostatistics 8, 158–183.
Wakefield, J.C. (2009). Multi-level modelling, the ecologic fallacy, and hybrid study designs. International Journalof Epidemiology 38, 330–336.
Wand, M.P. and Ormerod, J.T. (2008). On semiparametric regression with O’Sullivan penalised splines. Australianand New Zealand Journal of Statistics 50, 179–198.
Zeger, S.L. and Karim, M.R. (1991). Generalized linear models with random effects: A Gibbs sampling approach.
Journal of the American Statistical Association 86, 79–86.