Approximate Inference in Generalized Linear Mixed Models
N. E. Breslow; D. G. Clayton
Journal of the American Statistical Association, Vol. 88, No. 421. (Mar., 1993), pp. 9-25.
Stable URL:
http://links.jstor.org/sici?sici=0162-1459%28199303%2988%3A421%3C9%3AAIIGLM%3E2.0.CO%3B2-T
Journal of the American Statistical Association is currently published by American Statistical Association.

Your use of the JSTOR archive indicates your acceptance of JSTOR's Terms and Conditions of Use, available at
http://www.jstor.org/about/terms.html. JSTOR's Terms and Conditions of Use provides, in part, that unless you have obtained
prior permission, you may not download an entire issue of a journal or multiple copies of articles, and you may use content in
the JSTOR archive only for your personal, non-commercial use.
Please contact the publisher regarding any further use of this work. Publisher contact information may be obtained at
http://www.jstor.org/journals/astata.html.
Each copy of any part of a JSTOR transmission must contain the same copyright notice that appears on the screen or printed
page of such transmission.

The JSTOR Archive is a trusted digital repository providing for long-term preservation and access to leading academic
journals and scholarly literature from around the world. The Archive is supported by libraries, scholarly societies, publishers,
and foundations. It is an initiative of JSTOR, a not-for-profit organization with a mission to help the scholarly community takeadvantage of advances in technology. For more information regarding JSTOR, please contact support@jstor.org.

http://www.jstor.org
Mon Dec 10 03:44:51 2007

Approximate Inferencein Generalized Linear Mixed Models
N. E. BRESLOW and D. G. CLAYTON*
Statistical approaches to overdispersion, correlated errors, shrinkage estimation, and smoothing of regression relationships may be
encompassed within the framework of the generalized linear mixed model (GLMM). Given an unobserved vector of random effects,
observations are assumed to be conditionally independent with means that depend on the linear predictor through a specified linkfunction and conditional variances that are specified by a variance function, known prior weights and a scale factor. The random
effects are assumed to be normally distributed with mean zero and dispersion matrix depending on unknown variance components.
For problems involving time series, spatial aggregation and smoothing, the dispersion may be specified in terms of a rank deficient
inverse covariance matrix. Approximation of the marginal quasi-likelihood using Laplace's method leads eventually to estimating
equations based on penalized quasilikelihood or PQL for the mean parameters and pseudo-likelihood for the variances. Implementation
involves repeated calls to normal theory procedures for REML estimation in variance components problems. By means of informal
mathematical arguments, simulations and a series of worked examples, we conclude that PQL is of practical value for approximate
inference on parameters and realizations of random effects in the hierarchical model. The applications cover overdispersion inbinomial proportions of seed germination; longitudinal analysis of attack rates in epilepsy patients; smoothing of birth cohort effects
in an age-cohort model of breast cancer incidence; evaluation of curvature of birth cohort effects in a case-control study of childhood
cancer and obstetric radiation; spatial aggregation of lip cancer rates in Scottish counties; and the success of salamander matings in
a complicated experiment involving crossing of male and female effects. PQL tends to underestimate somewhat the variance componentsand (in absolute value) fixed effects when applied to clustered binary data, but the situation improves rapidly for binomial observations
having denominators greater than one.

KEY WORDS: Longitudinal data; Overdispersion; Penalized quasi-likelihood; Spatial aggregation; Variance components.

The generalized linear model (GLM) (McCullagh and
Nelder 1989) neatly synthesizes likelihood-based approaches
to regression analysis for a variety of outcome measures.
Several recent extensions of this useful theory involve models
with random terms in the linear predictor. Such generalized
linear mixed models (GLMM's) are useful for accommodating the overdispersion often observed among outcomesthat nominally have binomial (Williams 1982) or Poisson(Breslow 1984) distributions; for modeling the dependence
among outcome variables inherent in longitudinal or repeated measures designs (Stiratelli, Laird, and Ware 1984;
Zeger, Liang, and Albert 1988);and for producing shrinkage
estimates in multiparameter problems, such as the construction of maps of small area disease rates (Clayton and Kaldor
1987; Manton et al. 1989).
It is often a reasonable approximation and certainly traditional to assume that the random error terms have a multivariate normal distribution whose variance components areto be estimated from the data. When the outcomes come in
the form of proportions or counts, a full maximum likelihood
analysis based on their joint marginal distribution requires
numerical integration techniques for calculation of the loglikelihood, score equations, and information matrix. Thismethod has been implemented successfully in relatively
simple problems involving binomial (Brillinger and Preisler
1986; Crouch and Spiegelman 1990) and Poisson (Hinde1982) mixtures with a high degree of independence among
the observations. To date it has proved intractable for more
complicated problems involving irreducibly high-dimen-

sional integrals, however. Recent Bayesian procedures avoid
the need for numerical integration by taking repeated samples
from the posterior distributions using importance (Ii and
Raghunathan 1991) or Gibbs (Besag, York, and MolliC 1991;
Zeger and Karim 1991) sampling techniques. An attractive
feature of the Bayesian approach is its flexibility for full assessment of the uncertainty in the estimated random effects
and functions of model parameters. Potential drawbacks include the intensive computations and questions about when
the sampling process has achieved equilibrium (Ripley andKirkland 1990). There is still room for simple, approximate
methods both for exploratory analyses and to provide starting
values for use with other, more exact procedures.
This article considers two closely related approximate
methods of inference in GLMM's and investigates their suitability for practical work by means of Monte Carlo studies
and illustrative applications. Both have been considered previously, although not at the level of generality adopted here.
The penalized quasi-likelihood (PQL) method exploited by
Green (1987) for semiparametric regression analysis is available for inference in hierarchical models where the focus ison shrinkage estimation of the random effects (Robinson
1991). PQL was proposed as an approximate Bayes procedure for some commonly occurring GLMM's by Laird(1978) and by Stiratelli et al. (1984) and it has been used
more recently by Schall(199 1) and McGilchrist and Aisbett( 199 1). Marginal quasi-likelihood (MQL) is the name that
we give to the procedure proposed by Goldstein (1991) as
an extension to GLM's of his work on multilevel modeling* N. E. Breslow is Professor and Chairman, Department of Biostatistics,
University of Washington, Seattle, WA 98195. D. G. Clayton is Statistician,
MRC Biostatistics Unit, Cambridge CB2 2SR, United Kingdom. This work
was supported in part by United States Public Health Service Grant CA40644
and Senior International Fellowship TW01617. The authors thank the editors
and referees for helpful comments that substantiallyimproved the exposition.

O 1993 American Statistical Association
Journal of the American Statistical Association
March 1993, Vol. 88, NO. 421, Applications and Case Studies

Journal of the American Statistical Association. March 1993

10(Goldstein 1986, 1988). It is appropriate when interest is
focused more on the marginal relationship between covariables and outcome (Liang, Zeger, and Qaqish 1992; Zegeret al. 1988). A key feature of these methods is that they may
be implemented by repeated use of standard software forvariance components analysis of normally distributed observations, just as GLM's may be fitted by repeated calls toweighted least squares procedures. Both provide for shrinkage
estimates of random error terms and for analogs of restrictedmaximum likelihood (REML) estimates of the variancecomponents (Patterson and Thompson 197 1).
The organization of the article is as follows. After a brief
statement of the model in Section 1, the PQL criterion ismotivated in Section 2 by approximating the integrated
quasi-likelihood. In Section 3 an approximate GLM for themarginal distribution of the data is developed and related to
the generalized estimating equation approach of Zeger et al.
(1988). An outline of the computational procedure is presented in Section 4. A Monte Carlo investigation using binomial-normal observations is discussed in Section 5. Illustrative analyses presented in Section 6 indicate a wide range
of applications. In Section 7, PQL and MQL are related toother recent work and some suggestions for further research
are made. Readers who are less interested in theoretical andcomputational aspects may wish to concentrate their attention on Sections 1, 2.3, 3.1, 5, and 6.
I. THE HIERARCHICAL MODEL

Observations on the ith of n units consist of a univariate
response variable yi together with vectors xi and zi of explanatory variables associated with the fixed and random
effects. Units may be blocked in some way, for example when
they involve repeated measures on the same subject. We
suppose that, given a q-dimensional vector b of random effects, the yi are conditionally independent with means
E(yi ( b) = pP and variances var(yi ( b) = 4aiv(pP), wherev( ) is a specified variance function, ai is a known constant(e.g., the reciprocal of a binomial denominator) and 4 is adispersion parameter that may or may not be known. This
formulation encompasses situations where the random effects
are nested within subjects (e.g., Sections 6.1, 6.5) and when
they are not (Sections 6.2-6.4, 6.6). The conditional mean
is related to the linear predictor 11; = xfa zf b by the linkfunction g(pP) = ? f , with inverse h = g-', where a is a p
vector of fixed effects. Denoting the observation vector by y
= ( y l , . . . , y,)' and the design matrices with rows xf andzf by X and Z , the conditional mean satisfies.

+

The model is completed by the assumption that b has a multivariate normal distribution with mean 0 and covariancematrix D = D(0) depending on an unknown vector 8 ofvariance components. In all the examples we consider, which
involve binomial, Poisson, and hypergeometric specifications
for the conditional distribution of yi , the dispersion parameter 4 is fixed at unity. In other applications, however, it
may be estimated together with 0 as a parameter in the covariance matrix of the marginal distribution of y .

The integrated quasi-likelihood function used to estimate( a , 8) is defined bye4'(a'o)a~DI-'~'S 2"

1
exp[- - d i ( y i ;$ ) - - b'D-'b
24 , = ,
2

where

I

db,
(2)

denotes the deviance measure of fit. If, conditionally on b,
the observations are drawn from a linear exponential familywith variance function v( ), then the deviance is well known
to equal to the scaled difference 24{1(y; y, 4 ) - l(y; p, 4 ) ) ,where I(y; p, 4 ) denotes the conditional likelihood of y given
its mean p (see, for example, McCullagh and Nelder 1989).
In this case ql(a, 8) represents the true log-likelihood of thedata. As already mentioned, the primary difficulty in implementing full likelihood inference lies in the integrations
needed to evaluate ql and its partial derivatives.

.

2. PENALIZED QUASI-LIKELIHOOD
2.1 Motivation of the PQL Criterion
Writing Equation (2) in the form c(DI-""s-"'~'
db,
we apply Laplace's method for integral approximation(Barndorff-Nielsen and Cox 1989, sec. 3.3; Tierney and Kadane 1986). Let K' and K" denote the q vector and q X qdimensional matrix of first- and second-order partial derivatives of K with respect to b. Ignoring the multiplicative constant c, the approximation yields

where 6 = & ( a ,8) denotes the solution tothat minimizes K(b ) . Differentiating again with respect to b ,we haverz: Z'WZ+ D-I,(4)

where W is the n X n diagonal matrix with diagonal terms
wi = { +aiv(pf ) [gr(p; ) I ) that are recognizable as the
GLM iterated weights (Firth 1991, p. 63; McCullagh and
Nelder 1989, sec. 2.5). The remainder term

' -'has expectation 0 and is thus, in probability as a function ofn, of lower order than the two leading terms in Equation(4). R equals 0 for the canonical link functions, for which

Breslow and Clayton: Approximate Inference for GLMM's

II

gl(y) = v-'(y) (McCullagh and Nelder 1989, p. 32). Combining (2)-(4) and ignoring R leads towhere b is chosen to maximize the sum of the last two terms.
Assuming that the GLM iterative weights vary slowly (ornot at all) as a function of the mean, we ignore the first termin this expression and choose a to maximize the second.
Thus (&, b) = (&(8),b(8)), where b(8) = b(&(8)),jointly
maximize Green's ( 1 987) PQL

Differentiation with respect to a and b leads to the scoreequations for the mean parameters:andStiratelli et al. (1984) derived these equations for logistic
regression of binary data by maximizing the Bayes posterior
distribution for ( a , b) under a diffuse prior for a (cf. Schall1991).

2.2 Fisher Scoring
Green (1987) developed the Fisher scoring algorithm forsolution of Equations (7) and (8) as an iterated weighted
least squares (IWLS) problem involving a working dependent
variable and a weight matrix that are updated at each iteration. His development is modified slightly here to exploit
the close correspondence with the normal theory calculationsof Harville (1977). Defining the working vector Y to havecomponents Yi = ? ( yi - y? )g1(y: ), the solution to (7)and (8) via Fisher scoring may be expressed as the iterative
solution to the system

+

X'WX
[Z'WX

X'WZD ] ( a ) = [X'WY]
I+ZrWZD v
Z'WY

(9)

where b = Dv . Harville (1977) derived (9) for the best linear
unbiased estimation (BLUE) of a and b in the associatednormal theory model Y = Xa + Zb e, where e -- N(0,W-I ) and b -- N ( 0 , D), e and b independent. Equivalently,
one may first solve for a in+

(X'V-' x ) a = X'V-'Y,
where V = W(10)-' + ZDZ', and then set
b = D; = DZ'V-'(Y - x&).

theory linear model, provided that 8 is known. Standard errors for b may be calculated from (1 I). But both sets ofstandard errors ignore the additional variability stemmingfrom the need to estimate 8. Normal theory approximations
that account for this additional variability have been proposed (Kackar and Harville 1984), and it would be of interest
to explore their suitability for the more general models considered here.

2.3 Singular Variance Matrices
As seen in the examples in Sections 6.3-6.5, which involve
autoregressive models for smoothing and spatial aggregation,
it is sometimes useful to specify the random effects modelin terms of the inverse dispersion matrix R rather than D.
Reciprocals of the diagonal terms in R equal conditionalvariances given the remaining variables, whereas the off-diagonal terms determine the conditional regression and correlation relationships (Dempster 1972; Whittaker 1990). In
stochastic smoothing models of this type, the term b'Rb in
Equation (6) represents a "roughness penalty" (Good andGaskins 197 1). Typically, however, no penalty is associatedwith the overall average or level of the components of b, norsometimes even with their linear trend. Then R is singularand the probability distribution of the random effects is not
fully specified.
Suppose that r linearly independent combinations of b,
the "aliased" components, have no associated probability
distribution. Typically one constrains these components to
0, say by Gb = 0 where G is a r X q matrix of rank r, andincludes them instead in the fixed part of a model. There
are then two ways to proceed. Either one can reduce theproblem to that of estimating a reduced set of q - r random
effects that have a specified, full rank probability distribution(e.g., McGilchrist and Aisbett 1991). Or equivalently, one
can use for the dispersion matrix D the Moore-Penrose generalized inverse, R- (Graybill 1983, sec. 6.2.1). This sets to
0 the variances for the aliased linear combinations of b, sothat they are effectively constrained to take on the value 0.
The latter approach offers the greater flexibility and is theone chosen here.
2.4 Variance Component Estimation

Substitution of the maximized value of (6) into (5) andevaluation of W at (&(8),b(8)) generates an approximate
profile quasi-likelihood function for inference on 8. We makesome further approximations to motivate standard estimating equations in terms of the working vector Y, the iterated
weights W , and the design matrices X and Z . Ignoring
throughout the dependence of W on 8 and replacing thedeviance C di ( yi , y? ) by the Pearson chi-squared statistic
C (yi - y ~ ) 2 / [ u l v ()],y ~we have up to the usual additive
constantql(&(e),e)(1 1)This suggests that one take as an approximate covariancefor & the matrix (x'V-'X)-'. This is exact for the normal-

1--iogp I

2

1--(Y-X&)'v-'(Y-X&).
2(12)

Journal of the American Statistical Association, March 1993

12

This quantity, whose derivation used Harville's (1977)
Equations 5.1 and 5.2, may be recognized as the profile likelihood based on the associated normal theory model for Y.
To make degrees-of-freedom adjustments that account for
the fact that & rather than a appears in the quadratic formin (12), we use in practice the REML version (Patterson andThompson 1971):

become more normally distributed. Such "small dispersion
asymptotics" occur, for example, as the denominators ofbinomial proportions or the means of Poisson observations
increase (Jorgensen 1987).

11
qll (&(8), 8) = - - log ( V I - - log 1 X'V-'X I
2
2

A key feature of the hierarchical model is that the regression structure in Equation (1) is conditional on the valuesof the random effects b. When separate random effects are
estimated for each person in a medical study, for example,
this means that a represents covariable effects at the level ofthe individual subject. Because one generally desires estimates of covariable effects on population averages (e.g., on
the survival rates of specific subgroups), it is often more appropriate to specify the GLM in terms of the marginal meanas

For the normal theory linear model this adjustment corresponds to the profile likelihood correction of Cox and Reid
(1987). Full justification requires (a) that a and 8 be orthogonal parameters and (b) that the information matrix for&(8)be X'V-'X. Neither requirement holds exactly for thegeneral GLMM developed in this section. Both do hold,however, for the marginal model considered in the next section.
Following Harville (1977), we define P = V-' - V -'X
X (X'V-'X)-'X'V-' and differentiate (13) with respect tothe components of 8 to obtain estimating equations for thevariance parameters:dV(Y - Xa)'V-' - V-'(Y - Xa) - traej

The corresponding Fisher information matrix J has components

Note that because the dependence of W on 8 is ignored in
calculating dV/dO,, in (14) and (15) qll cannot be used as an
objective function to help solve the equations.

2.5 Remarks on Asymptotic Theory
Our "derivation" of the penalized quasi-likelihood (6) andmodified profile quasi-likelihood (13) involved several ad
hoc adjustments and approximations for which no formal
justification was given. It is best viewed as providing heuristic
motivation for the estimating Equations (7), (8), and (14)that may be studied in their own right. Informal considerations, however, do suggest circumstances in which the approximations should perform well. First, the equations are
the REML equations under the normal theory linear model,
for which the working and observation vectors coincide and
W is the identity matrix. Unless y can be partitioned into K
independent components where K increases with n, ofcourse, even this does not guarantee that the standard
asymptotic theory applies (Harville 1977, sec. 4.2). Second,
key portions of the argument involved approximating thedeviance increments by the normed, squared residuals or
the penalized deviance by a quadratic function of b. Bothapproximations are likely to improve as the individual yi

3. MARGINAL QUASI-LIKELIHOOD
3.1 The Marginal Model

Unless the link function is the identity, however, the marginal
mean so defined does not generally coincide with the marginal mean calculated from Equation (1).
One may, nonetheless, think of (16) as derived from arather crude, first-order approximation to the hierarchicalmodel that is valid in the limit as the components of dispersion approach 0. Writing the model in the form yi = pf
ci with var(ei ) = 4aiv(pf ) and b -- N ( 0 , D), one has y;
B h(x: a )hr(x:a)z: b t i (Goldstein 199 1). Defining Vo
and A to be the diagonal matrices with diagonal elements
4aiv(pi) and gl(pi) , the corresponding first-order varianceapproximation is

+

+

+

Marginal models of this form were investigated for longitudinal designs by Zeger et al. (1988). They showed thatthe true marginal mean for the hierarchical model with normally distributed random effects often could be expressedin the form of (16), at least approximately, but with alteredvalues for the regression variables or regression coefficients.
With the log link, for example, one finds E ( yi ) = exp(x; azfDzi/2), and thus that the random effects add an ofset
to the equation for the marginal mean. (Nonnormal random
effects lead to other offsets.) With the logit link, the a coefficients are attenuated:

+

+

+

where ci = I c 2 ~ z i z : I I - ' I 2 = (I
~ ~ z : ~ zand~ )c - ' ~ ~= 1 6 f i / ( 15a).They considered the approximation (17) sufficiently accurate for use as a "working covariance" matrix
in their iterative estimation procedure but based their inferences on an empirical covariance matrix derived from theestimating equations.
3.2 Estimation of Fixed and Random Effects
For fixed 8, we estimate the regression coefficients a inthe marginal model using the quasi-likelihood equations appropriate for dependent outcomes (McCullagh and Nelder

13Breslow and Clayton: Approximate Inference for GLMM's

1989, sec. 9.3). Denoting the marginal mean vector by p
= (pl,. . . , Iln)f,the estimating equations

take the form
With q = ( q l , . . . ,qn)'denoting the vector of linear predictors
q, = x ; a , Fisher scoring leads to IWLS regression of theworking vector Y = q A(y - p) on X with weight matrix+

expected to yield consistent estimates for 8. Using P rather
than V-I is intended to alleviate the small sample bias thatarises when & is substituted for a , and it does not affect theasymptotic argument.
The most rigorous demonstrations to date of these properties are those of Liang and Zeger (1986) and Prentice (1988)
for the case of block diagonal V. Their arguments lead toan alternative, empirical estimate of cov(&) that requires
only that the regression model (16) for the marginal means
be correctly specified.
4. COMPUTATIONAL ASPECTS
4.1 Initial Estimates of the Variance Parameters

where W is once again the diagonal matrix with the GLM
iterated weights
as diagonal elements. At each step in the iteration, theproblem is formally equivalent to that of estimating a
in the associated normal theory model Y = Xa Zb E ,with E -- N ( 0 , W-') and b -- N ( 0 , D). Shrinkage
estimates for the random effects are again obtained as b
= DZ'V-' (Y - X&).
The essential difference between the MQL estimating
equations for the marginal model and the PQL equations
for the hierarchical model is that the latter incorporate therandom effect terms z: b in the linear predictor.

++

3.3 Estimation of Variance Components

Variance parameters may be estimated using Carroll and
Ruppert's (1982) method of pseudolikelihood. Assuming that
E(y) is known, we consider the normal theory likelihood
based on the variance approximation (17) and take logarithmic derivatives with respect to 8. This leads to precisely thesame REML equations for the variance parameters as were
derived previously for PQL. Instead of iterating back andforth between (7)-(8) and (14), obtaining new values of b
at each iteration for use in Y and V, we iterate between (19)for a and ( 14) for 8, delaying the estimation of b from ( I I)
until convergence. Goldstein (199 I) and colleagues have implemented this procedure in the context of multilevel models
involving nested random effects. Approximate estimation of
the random effects is still a by-product of this approach, although their interpretation is not as clear as it is for thehierarchical model.
3.4 Asymptotic Justification

If (16) in fact correctly specifies the marginal mean of y,then the estimating equations (19) are unbiased so that, undersuitable regularity conditions, their root & is consistent fora and asymptotically normally distributed (McCullagh and
Nelder 1989, sec. 9.3). Provided also that (17) correctly
specifies the variance, the asymptotic covariance matrix may
be estimated by cov(&) = (x'V-'X)-'. Because E[dU/d8]= 0, furthermore, these asymptotic properties will continue
to hold even if the parameters in V are replaced by consistentestimates. Substituting V-I for P in (14), we see that these
also are unbiased estimating equations and hence may be

The computations associated with the two procedures arenearly identical and may be described in broad outline asfollows. Standard GLM techniques (McCullagh and Nelder
1989, sec. 2.5) lead to an initial estimate of a under theassumption that the n observations are independent ( D= 0). Residuals from this initial fit may be used to computeinitial values for the parameters in V, namely, 8 and also 4
if it is to be estimated. The exact procedure will depend on
the particular problem. Any reasonable method is likely to
be satisfactory for simple overdispersion models involving
independent observations and a single variance component.
We used moment equations (Moore 1986). In more complicated problems, the method of generalized least squares
may be tried. Here the upper triangular elements of Y *= (Y - X&)(Y - X&)'
X1(X'WX )-'X are regressed on
the appropriate design matrix for the variance components
using ordinary least squares. Goldstein (1986, 1988) has
shown that iterated generalized least squares, whereby thecovariance matrix of Y * is updated at each cycle and used
in a weighted analysis, produces REML estimates undernormal theory. In our experience with sparse discrete data,
however, the initial estimates obtained with ordinary leastsquares often were not satisfactory. In such cases we set theparameters corresponding to covariances to 0 and started
the variance parameters from small positive values.

+

4.2 Alternate Scoring for Means and Variances

The initial estimate 8' permits evaluation of the covariancematrix D that occurs in V. Updated values for (A,b) under
PQL are calculated from (7) and (8), or for & alone from(19) under MQL, with the fitted values from the GLM fit
being used to initialize the iterative solution of these equations by IWLS, as outlined in Section 2.2. Once the equations
for the mean parameters have been solved, the variancescores (14) and expected Hessian (15) are used to take a
Newton step towards a new value 8 ' , after which one returns
to (7)-(8) or (19) to reestimate the mean parameters. The
Newton step is halved if the resulting D(O1)is not positive
definite. This simple procedure led to joint solutions of themean and variance equations at an interior point of the parameter space in most of the examples considered. In general,
however, one must anticipate all the problems that attend
REML variance component estimation under standard normal theory. In the simulations we found that an adaptationof Marquardt's ( 1963) compromise between Fisher scoring

14

Journal of the American Statistical Association, March 1993

and the method of steepest descent, as described by Harville(1977), was useful in the initial stages of the search.
Evaluation of V-I and of the products of the n X n matrices
occurring in the gradient and expected Hessian for variancesis the most cumbersome aspect of the computation. Harville's
Equation (3.6), corresponding to Green's (5.1) and (5.2), is
useful for calculating V-' . Special structure in the covariancematrices must be exploited, however, if these procedures areto be used to estimate large numbers of random effects.
4.3 Scoring for REML Estimation in the Two-Level
ProblemLongford (1987), using formulas of Lamotte (1972), developed a "fast scoring algorithm" for maximum likelihood(ML) estimation of variance components in normal theorymodels involving nested random effects. We extended thedevelopment in his Appendix for the two-level problem toobtain REML rather than ML estimates. Further extensions
to accommodate more levels of nesting would be desirable.
Prosser, Rasbash, and Goldstein (1990) have implemented
the normal theory procedures for two and three levels using
iterated generalized least squares. Because the procedures we
describe involve repeated application of normal theory calculations to the working vector Y with associated covariancematrix V, such algorithms may be employed for approximate
inference in GLMM's (Goldstein 1991). The results reported
in Sections 5 and 6 were obtained using the matrix programming language GAUSS (Aptech Systems 1988).

5. A SIMULATION STUDY
The simulation study followed the design specificationsof Zeger and Karim (1991) to compare results with theirs
obtained using a Bayesian approach and the Gibbs sampler.
Each data set involved K = 100 clusters of size nk = 7. Conditionally independent binary observations ykl were generated
within each cluster with conditional response probabilities
given for k = 1, . . . , 100 bylogit E(yklIbk) = a 0+ a l t l + a2xk + a3xktl+ b t + bitl,

where xk = 1 for half the sample and 0 for the other half
and tl = I - 4 for I = 1, . . . , 7. The regression coefficientswere fixed at a' = (-2.5, 1, - 1, -.5) while the random effects(b:, bh) were generated as a series of 100 iid normal variables
with mean 0 and covariance structure(Zeger and Karim in fact used .49 rather than .50 for cooin
D2 and set a3 = .5 rather than -.5, but this has negligible
effect on the comparative results.) Two hundred data sets of
700 observations each were generated with Dl and 100 eachwere generated with D2.The entire experiment was replicatedwith binomial observations ykl whose denominators were m= 1 , 2 , 4 , 8.
Average values of the PQL regression coefficients were, inabsolute value, less than the true values but approached thelatter as the denominators of the individual.binomia1 observations increased (see Table 1). Positive estimates of coo= var(by) were obtained for all data sets sampled with D= D l . This single variance component was seriously underestimated when m = 1 but was only moderately underestimated when m = 8. With binary data generated under D2,
PQL frequently converged toward a non-positive definite
covariance matrix, in which case we set the smaller varianceand the covariance term, as well as the corresponding elements of their dispersion matrix, to 0. We estimated coo(resp.
a l l )to be 0 in 40% (resp. 12%) of samples with m = 1, 5%(8%) with m = 2, 1% (0%) with m = 4, and 0% (0%) with
m = 8. Moreover, PQL tended to underestimate both variance components. Rather large discrepancies were observed
between the simulated and estimated standard errors of thevariance components for small m (see Table 2). The standard
errors estimated for the regression coefficients, on other hand,agreed reasonably well with the simulated standard errors.
Even though the simulated data were generated under thehierarchical model, we also analyzed them under the marginal model using MQL to study the accuracy of the ap-

Table 1. Mean Values of Parameter Estimates in the Simulation Study
Parameterm

Method

uoo

501

01I

ffo

D = Dl (200replications)True value:
PQL

Gibbs samplera
D = D, (1 00 replications)

PQL

Gibbs samplera

True value:
1
2
4
8
1

From Zeger and Karim (1991).
See the text.

ffl

ff2

013

Breslow a n d Clayton: Approximate Inference for GLMM's

15

tained under PQL. With D2, in contrast, MQL tended tooverestimate aoo and to give a negative value (average of-.15) to sol. Such behavior would be anticipated from thediscussion in Section 3. According to (18), the "true" linear

Table 2. Comparison of Simulated and Estimated Standard Errors
Parameter
Methodm

2
4
8
Gibbs samplerC

1

goo

oo1

Estb
Sim
st
Sim
ESt
Sim
Est

.35 .24 .25 .20 .19 .18 .17 -

Sim
st

.63
.60

-

017

a
.

a,

a2

a3

- .28 .12 .46 .20- .22 .08 .31 .13- .22 .08 .34 .i4 where c:-.O6.18 .06,05.17 .05= ckl- Fand Fdenotes the average of the Ckl. This+- ,40 ,16 ,61 .25 into a negative covariance when modeled via the random- .36 .15 .56 .24 t e r m b ~ + b ~ t l .
6. ILLUSTRATIVE EXAMPLES

D = DP(100 replications)

PQL

1

8
Gibbs samplerc

1

+

.28
equals cG(-2.5
1.Otl) when xk = 0 and czl(- 1.5 1.Otl).28 .10
when
xk
=
1.
Becausein both cases the constant and t coef,24 ,07.25 .07 ficient are of opposite signs, the extra variability translates.lo .lo

Sima

.30

Est
Sim
Est
Sim
Est

.30 .I2.06.18 .08.14 .05.13 .06

Sim
Est.12 ,14 ,38.59 .21 .22 .37

.08
.06
.06
,04
.05

6.1 Overdispersion.28 .12 .4i .i9
.20
.15
,14
.13

.lo
.09
.09
.07
.08
,17
.19

.31
.24
.23
.22
.20
,54
.57

of seeds that germinated on each of 21 plates arranged according to a 2 X 2 factorial layout by seed variety and typeof root extract. He noted that the within-group variation
exceeded that predicted by binomial sampling theory, andhe was concerned that his logistic regression analysis oftreatment and interaction effects should account appropriately
for such overdispersion. One rather natural way of accounting for the extraneous plate-to-plate variability in this situation is by means of a GLMM that has linear predictor

.16
.13

.lo

.12
,31
.28

Simulated standard error.
ROO^ mean estimated variance.
From Zeger and Karim (1991).

logit Pr[y, = 1 / x i , b,] = 17f = x f a

proximation (18) and to evaluate the distortions in the estimated covariance matrix caused by the failure to model
the mean correctly. For this balanced design the MQL estimates of the fixed effects were identical to the standard
logistic regression estimates there were used as starting values(McCullagh and Nelder 1989, ex. 14.8 and 14.9). Thus, regardless of the binomial denominator m , & was estimating
the same quantity. With m = 8 and D l , the average valueof & was (-2.19, .87, -.92, -.40)', which may be comparedwith the approximate value of .862 X (-2.5, 1, - 1, -.5)
= (-2.16, .86, -.86, -.43) based on (18). The estimates ofaooand their standard errors were quite similar to those ob-

i = 1, . . . , 21, where a represents the fixed effects associatedwith seed and extract and the bi, assumed iid N(0, a2),
represent random effects associated with each plate. For the
agriculturalist interested in the effects of seed variety androot extract treatment on germination rates, it is more appropriate to model the marginal probabilities of germination(averaged over plates). In contrast, the hierarchical model isof interest in selecting plates containing subgroups of seeds
that may have particularly high germination rates.
Table 3 presents the regression coefficients in linear logistic
models fitted to the 21 binomial proportions of seed ger-

Table 3. Model Fits to Crowder's Seed Data
Method of analysis

Variable

LR a

PQL

MQL

ML

B+SE

B+SE

B+SE

B+SE

Main effects model
Constant
Seed (2)Extract (2)u-.430+ ,114-.270 k ,155
1.065+ ,144-.375 k .I82-.363 + .228
1.012 k .224
.352+ .I18

Constant
Seed (2)Extract (2)Interactionu-.558 + .I26
.I46+ ,223
1.318+ ,177-.778 + ,306-.542 k ,190
.077+ .308
1.339+ .270-.825 + ,430,313k ,121--.369 + ,180-.357+ ,227.998k ,222
,349+ ,117-.389 + ,166-.347 k .215
1.029 + ,205,295k .I12

lnteraction model-

LR = ordinary logistic regression.

+ bi,-.536 + ,190
.074+ .308
1.326 + ,269-.816 + ,429
,313k .I20-.548 k ,167
.097+ ,278
1.337+ ,237-.811 k .385
,236+ ,110

16

Journal of the American Statistical Association, March 1993

mination. The results for ML analysis under the hierarchicalmodel were obtained using the program EGRET (SERC
1989), which evaluates the integrated likelihood (2) and its
logarithmic derivatives using Gaussian quadrature. This gavesmaller estimates of the overdispersion variance component
than did PQL and MQL, and the interaction between seedtype and root extract appeared slightly more significant as aconsequence. The parameter estimates using MQL were noticeably attenuated in comparison with those for PQL, aswould be anticipated from Sections 3 and 5. From a practical
viewpoint, however, there is little to choose between the twoanalyses. Figure 1 presents the observed proportions andcorresponding fitted values from PQL under the main effects
model. There was substantial shrinkage toward the estimatedgroup means when the random effects were incorporated in
the linear predictor, especially for proportions with small
denominators.
Rotnitzky and Jewel1 (1990, sec. 4.4) considered robust
Wald and score tests for these data based on the empiricalvariance. Their Wald tests using a "working" exchangeable
correlation structure agreed quite well with those calculated
from the ML estimates and standard errors shown in Table
3. The corresponding tests for PQL and MQL were slightly
conservative in comparison.

analysis indicated that the counts were substantially lower
during the fourth visit and a binary variable (V4 = 1 for
fourth visit, 0 otherwise) was constructed to model such effects.
Thall and Vail's analysis focused on estimating the regression coefficients in a log-linear model for the marginal event
rates and the covariance parameters in various patterned
covariance matrices. Most of the latter were derived from ahierarchical, conditionally Poisson model that included one
set of independent random effects associated with each subject and another set associated with each visit. By carefully
examining residuals comparing the observed and fitted
counts of each subject at each visit, they identified as "outliers" a number of patients who either had particularly large
counts relative to the fitted model or had marked changes
over time in their counts, or both.
Our reanalysis of these data using a GLMM was oriented
primarily toward more systematic identification of patients
who had extreme levels, or extreme degrees of change overtime, in their event rates. We assumed that y,k, the seizure
count for patient j on the kth visit, was (conditionally) Poisson-distributed with mean d k . The most general model considered was

6.2 Longitudinal Data

where xjk denotes the vector of treatment, visit, and covariable main effects and interactions; ( bf , b j ) are bivariate normal random effects that represent the residual level and rateof change in the event rate for the jth subject; and the b$
are additional random error terms that represent nonspecificoverdispersion beyond that introduced by the subject-tosubject variation. In view of the discussion in Section 3.1,
the regression coefficients a for all but the constant term
may be interpreted also as log (relative) event rates in thecorresponding marginal model.
Table 4 presents the results obtained with PQL. Model I
involved a standard Poisson regression analysis, without accounting for the intrasubject correlation and overdispersion;
as a result, the standard errors of the subject level variables
were seriously in error. There was a marked attenuation in
the regression coefficients, especially the constant, and a sharp
increase in the standard errors when random subject effects
were introduced (Model 11) and a further reduction in the
fixed effect associated with the fourth visit with the additionof unit level random variation (Model 111). (Models I1 andI11 correspond to Models 52 and 42 of Thall and Vail; ournumerical results were similar to those shown in their table
4 for the closely related Model 22.)
Model IV was the most interesting. There was substantial
heterogeneity among subjects, even after accounting for thetreatment and baseline variables, both in the overall level
and in the trend in mean seizure counts. We estimated thecorrelation between these two random effects to be effectively
0. Figure 2 graphs the random effects estimated for the 59
subjects and identifies by number those with particularly
extreme levels or changes in attack rates, even after covariableadjustment. We identified patient 135 as having marked improvement over time after an initially high seizure rate, andidentified patients 227, 225, and 112 as having the highest

Thall and Vail(1990, table 2) presented data from a clinical trial of 59 epileptics who were randomized to a new drug
(Trt = 1) or a placebo (Trt = 0) as an adjuvant to the standard
chemotherapy. Baseline data available at entry into the trial
included the number of epileptic seizures recorded in thepreceding 8-week period and age in years. The logarithm ofI
7 the number of baseline seizures (Base) and the logarithmof age (Age) were treated as covariables in the analysis. A
multivariate response variable consisted of the counts of seizures during the 2-weeks before each of four clinic visits(Visit, coded Visit] = -3, - 1, 1, Visit4 = 3). Preliminary

.0

10

~

'

4

"

"

8

"

~

"

12

'

16

~

20

24

Observation Number

Figure 7. Proportion of Seeds That Germinated. The observed proportions are plotted as circles (a),followed by their denominators. Fitted
proportions under the fixed and random effects models are plotted as
squares ( E l ) and plusses (+) .

log /.lfk = xjka

+ bf + b? visitk/ 10 + byk,

Breslow a n d Clayton: Approximate Inference for GLMM's
Table 4. PQL Model Fits to Thall and Vail's Epilepsy Data
Model

Variable

I1

111

IV

pi
SE

B i SE

piSE

6iSE

+ .41

-1.25 k 1.2

-1.27

+ 1.2

-1.27 i1.2

.95 k .04
- 1 . 3 4 i .16.56 i.06.90 + .12-.I6 + .05-.87 k .14-.91 i .41.33 i .21
.47 + .36
- . I 6 k .05-.86 i .13-.93 i .40.34 k .21
.47 i .35-.lo + .09-.87 + .14-.91 i .41.33 i .21
.46 i .36
- . 2 6 i .16-.53 i .06-

.48 i .06-

.52 i .06
.74 i .16-.01 k .03--.36 k .04

-

-2.76

Fixed effects
Constant
Base
Trt
Base x Trt
Age
V4Visit110
Subject level random effects
Constant
Visit110
Covariance (u,,)

(G)
(G)

Unit level random effects

ConstantI(G)overall count levels relative to expectation based on the covariables. Because we expressed the random effects on thesame log scale as the fixed effects, moreover, we also identified
patients with especially low or zero counts (e.g., patient 232)
that were not so apparent in Thall and Vail's analysis.
6.3

Smoothing of Birth Cohort Effects in an
Age-Cohort Model

Breslow and Day (1975, table 2) analyzed breast cancer
rates in Iceland according to year of birth in K = 1 1 cohorts
from 1840-1 849 to 1940-1 949 and age in J = 13 groups
from 20-24 years to 80-84 years. They fitted a log-linear
model with fixed effects for age and cohort, both treated asfactors, to the two-way table of rates using Poisson regression
with the logarithm of the person-years denominators as anoffset in the regression equation. This yielded an excellent

I -1.0

-0.6

-0.2

0.2

0.6

1 .O

1.4

Intercept ( b l )Figure 2. Random Intercepts and Slopes for Epilepsy Patients. Individual
patients are identified according to the ID number in Table 2 of Thall andvail(1990).

fit, with a deviance of49.7 on 54 degrees of freedom. Because
case ascertainment was limited to the period 19 10- 1971, nodata were recorded for the younger age groups in the older
cohorts nor for the older age groups in the recent cohorts;
63 of the 11 X 13 = 143 cells in the two-way layout were
empty. Moreover, the 1840 (1 1 cases) and 1940 (7 cases)
cohorts were so small that the fixed effects estimated for themwere particularly unstable. To try to reduce the random error
in the relative risks estimated for these extreme cohorts, we
fitted a log-linear model with a single term for the logarithmof birth cohort number in place of the cohort factor. This
increased the deviance by only 10.8, with an increase of 9
degrees of freedom. Figure 3 plots the two sets of fixed cohort
effects, with the fifth (1880) cohort selected as baseline.
Despite the excellent fit of the linear model, however, there
was concern that the strong parametric assumption was unwarranted. To explore this possibility, we used a nonparametric smoother based on a GLMM with an autoregressive
error component. The logarithm of the mean number ofbreast cancer cases in the jth age group and kth birth cohort
was assumed to satisfy

where n,k denotes the person-years denominator and aj denotes the fixed effect of age. The fixed effect trend 0and therandom effect vectors u and v modeled three aspects of thevariation of rates with date of birth.
The effects vk were assumed to be independent N(O, I),
modeling unstructured heterogeneity of risk over birth cohorts. The model for u, intended to represent smooth variation over time, was specified in the forward direction as a
Gaussian autoregressive model. We initially considered asimple random walk model, as is used commonly in Bayesian
forecasting (Harrison and Stevens 1976). But to improve
performance at the endpoints, we decided instead on themodel in which each point is predicted by linear extrapolation from its two immediate predecessors rather than from

18

Journal of the American Statistical Association, March 1993

so that the conditional expectation is obtained by cubic interpolation from the two points on either side. The first (andlast) point in the series may be seen to have conditional
expectations given by linear extrapolation from the adjacent
two points and conditional variance 1. The second (and penultimate) point has conditional expectation intermediate between quadratic and linear interpolation, which would correspond torowsofRof(-1, 2, - l , O , . . .)and(-1, 3, -3,
1, . . .). R has rank K - 2 in this model, reflecting the factthat both level and trend are aliased. The former is taken up
by the period parameters aj. The latter is included as a fixed
effect, pk, for technical reasons described in Section 2.3.
Not surprisingly, in view of the goodness of fit of the modelwith a strong linear trend in cohort effects, the estimate forul converged to 0 when both u and v terms were includedin the model. When the autoregressive component alone
Log birth cohort number
was included, we estimated Go = .12 f .06 using PQL. With
Figure 3. Relative Rates of Breast Cancer Incidence in Iceland, by Coded independent random effects alone, excluding also pk, weYear of Birth. The logarithms of fitted rate ratios, using birth cohort number
found GI = .69 _t .17. Figure 3 shows the cohort effects es5 as standard, are plotted for the fixed-effects model as squares (a),for
the independent random effects model as plusses (+), and for the au- timated under these two mixed models as well as those fortoregressive random effects model as diamonds (6).The straight line the fixed effects models. The greatest differences between the(-)
represents the fitted rate ratios for the fixed-effects model with a estimated random effects occur at the endpoints, where thelinear term for the logarithm of the coded birth cohort.
autoregressive effects are fitted by linear extrapolation but
the independent effects are pulled in towards the common
mean. Differences also are evident at interior points, where
just the latest one. Thus for k > 2 we assumed
the log-relative risk is pulled towards the local cubic fit byE(ukIUj,j < k) = 2~k-I- uk-2
the autoregressive model and towards the overall mean byandthe independence model. Both mixed models yield morevar(uk)uj,j < k) = 1,(21) reasonable estimates than does the highly variable fixed effectat the upper endpoint.
which treats the second dzferences of the series as independent N ( 0 , 1) variates. The inverse dispersion matrix for u is 6.4 A Mixed Model for the Log Odds Ratio
given byKneale (197 1) classified deaths from childhood cancer andmatched controls in the Oxford region by age at death, yearof birth, and whether or not the mother reported having
-2
1
0 ...
received pelvic radiation during pregnancy. These data werearranged by Breslow (1976) into a series of 2 X 2 tables ofcases versus controls and x-ray versus no x-ray, one table foreach of 120 combinations of age 0-9 and birth year 1944. . . . . . . . . etc.
1964. Conditioning on the marginal totals in each table, he
fitted log-odds-ratio regression models using maximum likelihood techniques based on the noncentral hypergeometric
distribution (Zelen 1971). Although the relative risk for radiation was reasonably constant across age groups, there wasa marked decline with year of birth, which was interpreted
as an effect of temporal decline in radiation dose. The fixedeffects analysis seemed to suggest some curvature, as measured by a quadratic term in the regression of the log-relativerisk on year of birth (Fig. 4). But doubts about the significanceof such curvature were raised by the excess scatter of theEquation (23) gives an alternative specification of the individually estimated relative risks (shown as "fixed effects"
model as an undirected conditional regression model. For 2 in Fig. 4) about the regression line. This source of variation< k < K - 1,
had not been formally considered in the earlier analysis.
Mixed models for log-odds-ratio regression are easily accommodated within this article's general framework. Although the exact conditional means and variances of theandnoncentral hypergeometric required for (quasi)likelihood
analysis are prohibitively complicated when the marginal
totals are large, simple approximations suggested by Mc-

Breslow and Clayton: Approximate Inference for GLMM's

Year of birthFigure 4. Relative Risks of Childhood Cancer in the Oxford Region forChildren Exposed to In Utero Irradiation Versus the Nonexposed, by Yearof Birth. The logarithms of fitted relative risks are plotted for the fixedeffects model as squares (m) and for the autoregressive random effects
The curved line (-)
represents the fitted
model as diamonds (0).
relative risks in a fixed-effects model with linear and quadratic terms foryear of birth.

Cullagh (1984) are accurate enough for most practical purposes (Breslow and Cologne 1986). Denote the odds ratioof the expected values in the ith table by $i and denote thetwo sets of marginal totals by (ni 1, ni2) and (mi 1, mi2),
where ni ni = mi mi = Ni . Rather than parameterize
the model in terms of the mean pi = E(yi), where yi is theobserved frequency in the upper left cell, it is more convenient here to parameterize it directly in terms of the odds
ratio. Thus we suppose that I); = exp(x:a zf b) under
PQL or qi = exp(xfa ) under MQL. Using McCullagh's approximation, the mean pi and variance vi jointly satisfy theequations

+

+

+

19

where lCjk represents the log-relative risk of radiation in the2 X 2 table formed for the jth age group and the kth yearof birth. Yeark is coded -10 for 1944, -9 for 1945,. . . , 10
for 1964, and pk is an iid N(0, 1) error term representing
extraneous year-to-year variation. For a = 0 the results wereidentical to those reported by Breslow and Cologne (1986,
table 3) for the McCullagh approximation. Note the marked
decrease in the estimated variance component as the linear
birth cohort effect is added to the fixed part of the model.
There is no clear evidence for extraneous variation of the
log-relative risk about the regression line. Even if such variation is accounted for in the model, however, the coefficientof the quadratic term remains statistically significant.
We also fitted an autoregressive model, without the quadratic term, to estimate "nonparametrically" the evolutionof risk with time. Using the specifications of Section 6.3 for
pk, this model clearly showed the flattening of risk as radiation dosage was controlled during the mid-l 950s (Fig. 4).
The fixed-effects quadratic fit was less satisfactory, because
it suggested a sharper increase in risk in the early 1960s,
which seemed more a consequence of the parametric formulation than of any clear feature of the data.
6.5 Spatial Aggregation in Scottish Lip Cancer
Rates
Clayton and Kaldor (1987) analyzed observed and expected numbers of lip cancer cases in the 56 counties of
Scotland with a view toward producing a map that would
display regional variations in cancer incidence yet avoid thepresentation of unstable rates for the smaller counties. The
expected numbers had been calculated accounting for thedifferent age distributions in the counties using a fixed-effects
multiplicative model but were regarded for purposes of analysis as constants based on an external set of standard rates.
Thus, conditional on a set of values bi representing countyspecific log-relative risks (i.e., standardized morbidity ratios,;r SMR's), the observed numbers of cases yi, i = 1, . . . , 56
were assumed to have independent Poisson distributions withmeans pf = niexp(a bi). Here the ni denote the expected
numbers and a, the grand mean, plays the role of the logarithm of the overall SMR, which one would expect to benear unity in view of the way the expected numbers were
derived.
Clayton and Kaldor (1987) proposed empirical Bayes estimation of the county-specific SMR's using several alternative assumptions about the distribution of the random
effects. Specification of these as a random sample from a+and+ (mil - Pi) + (ni2 - Pi)
Because log rC/ is the canonical parameter for the noncentral
hypergeometric distribution (McCullagh and Nelder 1989,
sec. 7.3.2), the link derivative gl(p) equals the inverse variance v-'(p), so that the denominator terms drop out of the
PQL score equations (7) and (8). The link derivative so defined is needed, however, to calculate the GLM iterative
weights that enter into the PQL and MQL estimation procedures.
Table 5 shows the results of several PQL model fits to the
Oxford childhood cancer data using the previously described
approximation. The general form of the fitted models was
log $,k= a+ PIYeark+ P2(Year2 - 22) + apk,

Table 5. PQL Model Fits to the Oxford Childhood Cancer Data
Regression coefficients* standard error

Constant

YEAR

YEAR^ - 22

u

*
*
*
*
*

-

-

+ .07
.16 + .09,505 .056,531 .076.516 .056,536 +_ ,070
,565 ,061.566 ,070

+
+

,0144
-.0385-.0406 +_ .0162-.0445.0149-.0469,0167*+,0067 ,0030.0071 + ,0033

-

.23

-

-

.15

* .10

Journal of the American Statistical Association, March 199320

log-gamma distribution led to an exact analysis based on the
resulting independent negative binomial distributions for theyi . To account for spatial aggregation, they also considereda model with normally distributed random effects. Approximation of the conditionally Poisson log-likelihood by aquadratic expansion centered at the empirical log SMR, 6i= log[(d,;)/nil, led to a solution based on the EM algorithm (Dempster, Laird, and Rubin 1977). There was,however, an inconsistency in their specification of the spatial
process (Clayton and Bernadelli 1991).

+

Presumably the spatial aggregation is due in large part tothe effects of environmental risk factors. Data were available
on the percentage of the work force in each county employedin agriculture, fishing, or forestry. The authors who compiled
the data (Kemp, Boyle, Smans, and Muir 1985) noted that
this covariable, xi, exhibited spatial aggregation parallelingthat for lip cancer itself. Because all three occupations involve
outdoor work, the authors suggested that exposure to sunlight, the principal known risk factor for lip cancer, might
be the explanation. We analyzed the augmented data (Table

Table 6. Observed and Fitted SMR's for Lip Cancer in 56 Scottish Counties
Fitted SMR's
Co

Obs

EXP

Cov

Obs
SMR

NOTE: I, Clayton and Kaldor; II, POL, independence; and Ill, PQL, spatial correlation

I

11

111

Adjacent Counties

Breslow and Clayton: Approximate Inference for GLMM's

6) using a conditionally independent Poisson model with theconditional means for county i specified bylog p! = log ni+ a. + a l x i / 1 0 + bi.

Two separate assumptions were made regarding the distribution of the random effects bi: (a) iid N(0, c 2 ) and (b)Gaussian intrinsic autoregression. The latter specified
the inverse variance matrix R according to Besag et al.
(1991) and gave b an improper density proportional to-- jdenotesadjacente ~ p { - C ~ - ~-( bJ)2/(2u2)),whereib~counties. Here the conditional expectation of bi given b,, j# i equals the mean of the bj in contiguous counties, whereas
the conditional variance equals u 2 divided by their number.
Table 7 presents regression coefficients estimated by PQL
for models with and without the covariable x, and with both
independence and autoregressive structures for the random
effects. Note the slight attenuation in the coefficient of thecovariable when independent random effects are incorporated and the much more marked attenuation under the
spatial structure. Similarly, there is a reduction in the variance components when the covariable is included. This confirms that much of the spatial aggregation is explained by
the percentage of the work force engaged in outdoor occupations. Caution should be exercised in interpreting this result, however, because the model assumed that the random
effects were unrelated to the covariables. If the other factorsthat contributed to the spatial correlation were positively
correlated with the work force covariable, which seems likely,
then the regression coefficient estimated for this variable
would be attenuated, with some of the work force effect being
attributed instead to location.
Table 6 presents the SMRs estimated for each county under Clayton and Kaldor's (1987) approximation to the normal theory independence model (Model I), under PQL withindependent random effects (Model 11), and under PQL withspatial autocorrelation (Model 111). Because Clayton and
Kaldor had not considered the covariable, it was omitted
from the PQL fits to achieve comparability. All models realized the fundamental goal of pulling in the extreme SMRs
based on small numbers of observed cases. Estimated SMRs
for Models I1 and 111, however, were markedly different for
those counties for which the mean SMR of the neighbors
differed notably from the overall mean. For example, theSMR for county 19 (note that the counties are ordered according to the observed SMR) was pulled down from 162.7
to 153.6 by the independence model, but upwards to 204.2
by the spatial model. Its SMR was higher than average, but
those of the neighboring counties were higher still. For coun-

21

ties with 20 or more cases, the PQL estimates under independence (Model 11) agreed closely with those of Clayton
and Kaldor (Model I). For smaller counties, the PQL estimates were generally lower. The difference can be attributed
largely to their use of the .5 correction in the empiricallog SMR used to center the approximation, which acts toincrease this quantity above the level actually observed.
6.6 Crossed Random Effects: The Salamander
Data

McCullagh and Nelder (1989, sec. 14.5) published an interesting set of data on the success of matings between maleand female salamanders drawn from two populations, theroughbutts (RB) and the whitesides (WS), that had been geographically isolated from each other. In the first of three experiments, conducted during the summer of 1986, 10 RB
females and 10 WS females were mated with three RB males
and three WS males, for a total of six matings each over 24
days. Each of 10 RB males and 10 WS males likewise served
as mates for three females of each type. These same 40 salamanders were used in a repeat experiment conducted in
the fall that involved no repetitions of the earlier male-female
pairs. A third experiment, also conducted in the fall, used a
new set of 40 animals. Each experiment involved 30 matingsof each of the four gender-population combinations. Simple
inspection of the data revealed that three of the crosses had
success rates of about 7070,whereas the mating of WS femaleswith RB males was successful only 25% of the time. Evaluating the statistical significance of these differences was complicated by the fact that the 360 binary responses were notindependent.
McCullagh and Nelder considered a linear logistic modelfor the marginal probabilities of success, using a linearization
as in Section 3.1 to derive an approximate covariance matrix
for the vector of 360 binary outcomes. Because of the balanced design, their quasi-likelihood estimates of the regression parameters were identical to those obtained from standard logistic regression under independence (McCullagh and
Nelder 1989, ex. 14.8 and 14.9). They used a method of
moments procedure based on the observed covariance matrixof the residuals (observed binary minus fitted values) to estimate the variance components.
Karim and Zeger (1992) reanalyzed these data under ahierarchical model, using the Gibbs sampler with a noninformative prior on the dispersion matrix to approximate theposterior distributions of the parameters and random effects
of interest. Denoting by y f i the response for female i andmale j in experiment k, their model B may be written

Table 7. PQL Model Fits to the Scottish Lip Cancer Data
Regression coefficients +_ standard error
Constantx/10

u (spatial)u (independence)

where bf and bj" each have two components representing
the random effects associated with the indicated animal insummer and fall. The fixed effects consisted of a constant,an indicator WSF of whiteside females, an indicator WSMof
whiteside males, their interaction, and an indicator of fallseason. They also considered Model A, in which the genderspecific random effects from each experiment were assumed
independent with equal variances, and Model C, in which

22

Journal of t h e American Statistical Association, March 1993

separate fixed and random effects were estimated for eachexperiment. Models A and C were also considered by
McCullagh and Nelder. Neither contained the fixed effectfor fall season.
We fitted these same models with PQL. For Model B thevariance matrix of Equation (20) took the form
V = W-'+ z f D f ( z f ) '+ Z m D m ( Z m ) t ,where zf and Z m identified the female and male animals
involved in each mating. The covariance matrix D~of female
random effects was

Table 9. Regression Coefficients for the Salamander Data
Model A

Model BGibbs
sampler *PQLparameterbkSE(~)+ S E ( ~ ) 6 + SE(b)6B S E ( ~ )Constant
Fall season

1.03 k .43.79 k .32+

-2.29 + .43-.54+.39
2.82k.50

1.18 k .49-.50 k .41-2.43 + .44-.62k.46
3.01k.52

ws,

W ~ M
WSFXSWM

-

-3.01
.60
-.69k.50
3.74k.68-

Gibbs
sampler*

1.48 k .64-.62 + .54-3.13 + .62-.76+.62
3.90k.72

PQL' Zeger and Karim (1991), table 3. SE computed as range of 90% CI t 3.3.

,
where I is the identity matrix of dimension 20, a l l is thevariance of the summer effects, az2is the variance of the fall
effects, and a12is the covariance. D mhad a parallel structure.
Using Harville's (1977) Equation (3.6) and standard formulas
for the inversion of partitioned matrices, we calculated V-I
by inversion of matrices no larger than 60 X 60. But theneed to consider 360 X 360 matrices in the variance equations (14) and (15) meant that this problem was about aslarge as could be handled reasonably by a microcomputer.
Table 8 presents variance components found with the three
procedures. Those estimated by PQL and by McCullagh and
Nelder's moments method were in reasonable agreement,with the exception of the summer male variance. Variances
estimated by the Gibbs sampler were substantially larger, in
accordance with our simulation results. The covariance matrix for males under Model B converged towards singularity
with PQL; we thus constrained the summer and fall effects
for individual males to be equal, reducing the number of
estimated variance parameters from six to four. A high correlation between the summer and fall male effects was also
evident with the Gibbs sampler. Absolute values of theregression coefficients for Models A and B (Table 9) wereTable 8. Variance Components for the Salamander Data
Method of estimation
Moments8d

Seasona,

Gibbs
sampler

d

PQL

a,

d

a,

Model A: Pooled over independent experiments
Total

.91

.88

1.50

1.36

.72

.63

1.09
.62
-.I2

.9Oc
.9Oc
.9Oc

Model B: Correlated random effects
SummerFall
Covariance

-

-

1.92
1.37
-.25

1.25
2.02
1.52

Model C: Separate effects each experiment
Summer '86
Fall '86 rerun
Fall '86

1.37
.98
.40

.70
.60
1.34

McCullagh and Nelder (1989), table 14.10.
Karim and Zeger (1992), table 3, medians.
' constrained to be equal.

2.35
2.99
.33

.14
1.42
2.89

1.41
1.26
.26

.09
.62
1.50

25-30% larger under the Gibbs sampler than under PQL,
which also accords with our simulations (see Sec. 5). The
key interaction effect was highly significant regardless of themethod of analysis. MQL actually would have been more
appropriate in this example, because one is primarily concerned with estimating the marginal success rates associatedwith different types of matings, not in identifying individual
salamanders with particularly high success rates.

7. DISCUSSION AND CONCLUSIONSThe limited simulation results and illustrative analyses
presented in this article suggest that reasonably simple approximate methods are available for inference on random
effects in the context of GLM's. These methods may be implemented by repeated application of normal theory mixed
model procedures, and they are subject to the same computational limitations as regards sample size and variancestructure. When the random effects are nested, even at threeor more levels, algorithms developed by Goldstein (1986,
1988) and Longford (1987) suffice to treat large complex
problems. The procedures discussed herein are still rather
limited in their ability to handle crossed designs, however.
The salamander problem, involving two sets of 60 random
effects each, probably represents an approximate break evenpoint, after which the Bayesian treatment based on Gibbs
sampling becomes progressively easier to implement incomparison with PQL or MQL, due to their requirement
for manipulation of large matrices. Of course, such comparisons are highly dependent on the available computing
resources.
Our simulation results were encouraging as regards theability of PQL to render approximately correct inferences
on regression coefficients in hierarchical models. As anticipated by the development of Section 2, accuracy improved
as the binomial denominators increased. Nonetheless, evenwith binary data the results were sufficiently accurate formany practical purposes. Inference on variance parameters
was less satisfactory under PQL, due largely to the tendencyof the procedure to converge to a nonpositive definite variance matrix when the binomial denominator was 1 or 2.
When the response probabilities are small and the data arehighly discrete, only limited information is present for estimating random effects and their associated variances andcovariances. The Bayesian formulation enjoys an advantage
here because of the information on variance components

Breslow and Clayton: Approximate Inference for GLMM'scontributed by the prior distribution. It makes full use of thenormal theory distributional assumptions, however, whereas
PQL requires specification only of the first and second moments.
Another limitation of PQL, in common with empirical
Bayes methods more generally, is the failure to account forthe contribution of the estimated variance components whenassessing the uncertainty in both random and fixed effects.
One major distinction between PQL and MQL is the factthat the regression coefficients of the former, but not of thelatter, depend strongly on the estimated variance components
when the link function is not the identity, even in large samples. Development of an approximate covariance estimatefor ( 2 , 8) and refinements in the expressions given for thevariance estimates of & and 8 individually remain important
problems for further research. (In the meantime, some encouragement can perhaps be taken from the results forregression coefficients in Table 2.) Bootstrap methods (Lairdand Louis 1987) are used increasingly as a means of constructing confidence intervals for the estimated random effects, but these are controversial as regards their ability to
fully reflect the relevant uncertainties. Such problems arehandled more easily within the Bayesian context by examining posterior distributions.
Our simulation results for MQL, summarized briefly inthe text, confirm the anticipated attenuation in the estimatedregression coefficients when the binary data are generated
under the hierarchical model yet analyzed under the marginalone. A notable feature is the spurious correlation inducedin the estimated random effects by the failure to correctlymodel the heterogeneity in the fixed effects; see Equation(18). Much of the bias in the estimation of both mean andvariance parameters might well be alleviated for the logistic
model by treating the terms ci of (18) as a multiplicative
oflset whose values depended on the current 8. We plan further research to investigate this question. Results of anothersimulation study of Liang and Zeger's (1986) generalized
estimating equation (GEE) approach to clustered binary data
(Sharples and Breslow 1992), for which the marginal meansof the simulated data actually satisfied the assumed linear
logistic equation, suggested that the regression parameters
and their standard errors were well estimated, but the correlation parameters less so. Because the estimating equationsfor MQL and GEE are identical as regards the mean parameters, one might anticipate similar results for MQL. MQL(or GEE) is the method of choice when interest is focusedon the marginal relationship between covariables and response, and the random effects model serves mainly to suggest a plausible covariance structure, as expressed in V, thatenables one to get reasonably efficient estimating equations
for the mean value parameters. By contrast, PQL is themethod of choice for estimating parameters in the hierarchical model, especially when attention is focused on the
random effects.
Our work is closely related to that of several other researchteams. The equivalence of MQL and Goldstein's (199 1) procedure for GLMM's with nested random effects has been
mentioned already. Earlier work on this topic was reported
by Morton (1987, 1988). When applied to normally distrib-

23

uted data with nonlinear link functions, PQL is equivalent
to the nonlinear mixed model procedure of Lindstrom andBates (1990). Note in particular the correspondence between
their Equation (3.3) and our Equation (6) and between their(4.3)and our (12). Their development extends ours to general
nonlinear structures for the mean, whereas ours extends
theirs to nonnormal observations having defined meanvariance relationships. A synthesis is clearly possible. Since
this article was first submitted for publication, papers by
Schall(1991) and McGilchrist and Aisbett (1991) have appeared, both of which discussed PQL for GLMM's that involve sums of independent random effects. Schall (1991)
obtained estimates of regression coefficients and variancecomponents for Model C of the salamander data that wereidentical to ours; he did not consider the correlated random
effects Model B.
Liang and Waclawiw (1990) and Waclawiw and Liang(199 1) proposed an interesting estimating equation approachto the random effects in GLMM's, in which the fixed effects
were estimated from the approximating marginal model as
in MQL and were then corrected for attenuation using
Equation (1 8). Optimal estimating equations in the sense of
Godambe (1960) were suggested for the random effects.
When applied to the simple overdispersion problem involving a log-linear, conditionally Poisson model, their smoothed
estimate incorporating the random effect is a weighted average of the actual observation and the fitted mean. WithPQL and MQL, in contrast, such averaging is canied out
with the working vector Y on the scale of the linear predictor.
This approach seems more natural. Detailed comparisonsof the properties of the random effects predicted by these
two approaches would be a worthy subject for further research.
Other proposals have been made to introduce autoregressive variance structures into GLMM's. Zeger (1988) considered a time series model for the random effects in a simple
overdispersion model for count data. The conditionally loglinear mean structure was preserved in the marginal model.
He approximated the marginal variance matrix by a band
diagonal matrix having correlation structure appropriate forthe assumed autoregressive model, which was not precisely
preserved. This enabled him to avoid the large-scale matrix
inversion required by our approach to autoregressive structures. Development of similar approximations for other
GLMM's is a high priority. Goldstein, Healy, and Rasbash(1991) also considered an autoregressive structure for theunit level random effects in the context of multilevel mixed
models for continuous responses and estimated the nonlinear
correlation parameters by the method of constructed variables. There are close connections here with the Bayesian
approach to dynamic generalized linear models (West, Harrison, and Migon 1985).
The preceding citations give some indication of the current
high level of interest in GLMM's. Our goal has been to provide a unified framework in which much of this work may
be discussed and compared and to demonstrate by meansof simulations and worked examples its potential for applications.
[Received March 1991. Revised April 1992. ]24

Journal of the American Statistical Association, March 1993

REFERENCES
Aptech Systems (1988), The GAUSS SystemLVersion2.0, Kent, WA: Author.
Barndorff-Nielsen, 0 . E., and Cox, D. R. (1989), Asymptotic Techniquesfor Use in Statistics, London: Chapman and Hall.
Besag, J., York, J., and MolliC, A. (1991), "Bayesian Image Restoration,with Two Applications in Spatial Statistics" (with discussion), Annals of
the Institute of Statistical Mathematics, 43, 1-59.
Breslow, N. E. (1976), "Regression Analysis of the Log Odds Ratio: A Method
for Retrospective Studies," Biornetrics, 32, 409-4 11.
(1984), "Extra-Poisson Variation in Log-Linear Models," Applied
Statistics, 33, 38-44.
Breslow, N. E., and Cologne, J. (1986), "Methods of Estimation in Log
Odds Ratio Regression Models," Biornetrics, 42, 949-954.
Breslow, N. E., and Day, N. E. (1975), "Indirect Standardization and Multiplicative Models for Rates, With Reference to the Age Adjustment ofCancer Incidence and Relative Frequency Data," Journal ofchronic Diseases, 28, 289-303.
Brillinger, D., and Preisler, H. K. (1986), "Two Examples of Quanta1 Data
Analysis," in Proceedings ofthe 13th International Biornetrics Conference,
Seattle: The Biometric Society. pp. 94-1 13.
Carroll, R. J., and Ruppert, D. (1982), "Robust Estimation in Heteroscedastic
Linear Models," The Annals of Statistics, 10, 429-441.
Clayton, D., and Bernardinelli, L. (1991), "Bayesian Methods for Mapping
Disease Risk," in Small Area Studies in Geographical and Environmental
Epidemiology, eds. J. Cuzick and P. Elliot, Oxford: Oxford University
Press.
Clayton, D., and Kaldor, J. (1987), "Empirical Bayes Estimates of AgeStandardized Relative Risks for Use in Disease Mapping," Biornetrics,
43, 671-681.
Cox, D. R., and Reid, N. (1987), "Parameter Orthogonality and Approximate
Conditional Inference," Journal of the Royal Statistical Society, Ser. B,
49, 1-39.
Crouch, E. A. C., and Spiegelman, D. (1990), "The Evaluation of Integralsof the Form f(t)exp(-t2)dt: Application to Logistic Normal Models,"
Journal of the American Statistical Association, 85, 464-469.
Crowder, M. J. (1978), "Beta-Binomial ANOVA for Proportions," Applied
Statistics, 27, 34-37.
Dempster, A. P. (1972), "Covariance Selection," Biornetrics, 28, 157-175.
Dempster, A. P., Laird, N. M., and Rubin, D. B. (1972), "Maximum Likelihood From Incomplete Observations," Journal of the Royal StatisticalSociety, Ser. B, 39, 1-38.
Firth, D. (1991), "Generalized Linear Models," in Statistical Theory andModelling, eds. D. V. Hinkley, N. Reid, and E. J. Snell, London: Chapman
and Hall, pp. 55-82.
Godambe, V. P. (1960), "An Optimum Property of Regular Maximum
Likelihood Estimation," The Annals of Mathematical Statistics, 3 1, 12081211.
Goldstein, H. (1986), "Multilevel Mixed Linear Model Analysis Using Iterative Generalized Least Squares," Biometrika, 73, 43-56.
(1988), "Restricted Unbiased Iterative Generalized Least-Squares
Estimation," Biometrika. 76, 622-623.
(199 I), "Nonlinear Multilevel Models, With an Application to Discrete Response Data," Biometrika, 78, 45-5 1.
Goldstein, H., Healy, M., and Rasbash, J. (1992), "Time Series Models for
Repeated Measures Data," unpublished manuscript, submitted to Biometrics.
Good, I. J., and Gasluns, R. A. (1971), "Nonparametric Roughness Penalties
for Probability Densities," Biometrika, 58, 255-277.
Graybill, F. A. (1983), Matrices With Applications in Statistics (2nd ed.),Belmont, CA: Wadsworth.
Green, P. J. (1987), "Penalized Likelihood for General Semi-parametric
Regression Models," International Statistical Review, 55, 245-259.
Hamson, P. J., and Stevens, C. F. (1976), "Bayesian Forecasting" (withdiscussion), Journal ofthe Royal Statistical Society, Ser. B, 38, 205-247.
Harville, D. A. (1977), "Maximum Likelihood Approaches to Variance
Component Estimation and to Related Problems," Journal ofthe American Statistical Association, 72, 320-340.
Hinde, J. (1982), "Compound Poisson Regression Models," in GLIM 82:
Proceedings ofthe International Conference on Generalized Linear Models.
ed. R. Gilch&t, Berlin: Springer, pp. 109-121.
Ii, Y., and Raghunathan, T. E. (1991), "Bayesian Analysis of Series of2 x 2 ~ ~ bunpublishedl ~ ~manuscript, ~submitted to statistics in ~ ~ d iJorgensen, B. (1987), "Exponential Dispersion Models" (with discussion),
Journal ofthe Royal Statistical Society, Ser. B, 49, 127-162.
Kackar, R. N., and Harville, D. A. (1984), "Approximations for Standard
Errors of Estimators of Fixed and Random Effects in Mixed Linear Models," Journal ofthe American Statistical Association, 79, 853-862.

S

Karim, M. R., and Zeger, S. L. (1992), "Generalized Linear Models With
Random Effects: Salamander Mating Revisited," Biornetrics, 48, 631644.
Kemp, I., Boyle, P., Smans, M., and Muir, C. (1985), Atlas of Cancer in
Scotland, 1975-1980. Incidence and Epidemiologic Perspective, IARC
Scientific Publication 72, Lyon, France: International Agency for Research
on Cancer.
Kneale, G. W. (197 I), "Problems Arising in Estimating From Retrospective
Survey Data the Latent Period of Juvenile Cancers Initiated by ObstetricRadiography," Biornetrics, 27, 563-590.
Laird, N. M. (1978), "Empirical Bayes Methods for Two-Way Contingency
Tables," Biometrika, 65, 581-590.
Laird, N. M., and Louis, T. A. (1987), "Empirical Bayes Confidence Intervals
Based on Bootstrap Samples," Journal ofthe American Statistical Association, 82, 739-757.
Lamotte, L. R. (1972), "Notes on the Covariance Matrix of a Random,
Nested ANOVA Model," The Annals ofMathematica1 Statistics, 43, 659662.
Liang, K. Y., and Waclawiw, M. A. (1990), "Extension of the Stein Estimating Procedure Through the Use of Estimating Functions," Journalof the American Statistical Association, 85, 435-440.
Liang, K. Y., and Zeger, S. L. (1986), "Longitudinal Data Analysis Using
Generalized Linear Models," Biometrika, 73, 13-22.
Liang, K. Y., Zeger, S. L., and Qaqish, B. (1992), "Multivariate Regression
Analyses for Categorical Data" (with discussion), Journal ofthe Royal
Statistical Society, Ser. B, 54, 3-40.
Lindstrom, M. J., and Bates, D. M. (1990), "Nonlinear Mixed Effects Models
for Repeated Measures Data," Biornetrics, 46, 673-687.
Longford, N. T. (1987), "A Fast Scoring Algorithm for Maximum Likelihood
Estimation in Unbalanced Mixed Models With Nested Random Effects,"
Biometrika, 74, 8 17-827.
Manton, K. G., Woodbury, M. A,, Stallard, E., Riggan, W. B., Creason,
J. P., and Pellom, A. C. (1989), "Empirical Bayes Procedures for Stabilizing
Maps of U.S. Cancer Mortality Rates," Journal of the American Statistical
Association, 84, 637-650.
Marquardt, D. W. (1963), "An Algorithm for Least Squares Estimation of
Nonlinear Parameters," SIAM Journal, 11, 43 1-44 1.
McCullagh, P. (1984), "On the Elimination of Nuisance Parameters in the
Proportional Odds Model," Journal ofthe Royal Statistical Society, Ser.
B, 46, 250-256.
McCullagh, P., and Nelder, J. A. (1989), Generalized Linear Models (2nded.), London: Chapman and Hall.
McGilchrist, C. A,, and Aisbett, C. W. (1991), "Restricted BLUP for Mixed
Linear Models," Biometrical Journal, 33, 131-14 1.
Moore, D. F. (1986), "Asymptotic Properties of Moment Estimators for
Overdispersed Counts and Proportions," Biometrika, 73, 583-588.
Morton, R. (1987), "A Generalized Linear Model With Nested Strata ofExtra-Poisson Variation," Biometrika. 74, 247-257.
(1988), "Analysis of Generalized Linear Models With Nested Strata
of Variation," Australian Journal of Statistics, 30A, 2 15-224.
Patterson, H. D., and Thompson, R. (1974), "Recovery of Interblock Information When Block Sizes Are Unequal," Biometrika, 58, 545-554.
Prentice, R. L. (1988), "Correlated Binary Regression With Covariates Specific to Each Binary Observation," Biometrics, 44, 1033-1048.
Prosser, R., Rasbash, J., and Goldstein, H. (199 I), ML3 Sofwarefor ThreeLevel Analysis, Users' Guide for V. 2, London: Institute of Education.
Ripley, B. D., and Kirkland, M. D. (1990), "Iterative Simulation Methods,"
Journal of Computational and Applied Mathematics, 3 1, 165-172.
Robinson, G. K. (1991), "That BLUP Is a Good Thing: The Estimation of
Random Effects." Statistical Science., 6., 15-5 1.
Rotnitzky, A,, and Jewell, N. P. (1990), "Hypothesis Testing of Regression
Parameters in Semiparametric Generalized Linear Models for Cluster
Correlated Data," Biometrika, 77, 484-497.
Schall, R. (1991), "Estimation in Generalized Linear Models With Random
Effects," Biometrika, 40, 7 19-727.
SERC (1989), EGRET Users' Manual, Seattle: Author.
Sharpies, K., and Breslow, N. (1992), "Regression Analysis of Correlated
Binary Data: Some Small Sample Results for the Estimating Equation
Approach," Journal of Statistical Computation and Simulation, 42, 1A A

LU.

Stiratelli, R., Laird, N., and Ware, J. (1984), "Random Effects Models forWith Binary Responses," Biometrics. 40, 96 1-97 1.
~ Serial
i ~ Observations~ ,Thall, P. F., and Vail, S. C. (1990), "Some Covariance Models for Longitudinal Count Data With Overdispersion," Biometrics. 46, 657-67 1.
Tierney, L., and Kadane, J. B. (1986), "Accurate Approximations for Posterior Moments and Marginal Densities," Journal ofthe American Statistical Association, 8 1, 82-86.

reslow and Clayton: Approximate Inference for GLMM's

Vaclawiw, M. A., and Liang, K. Y. (199 I), "Prediction of Random Effects
in the Generalized Linear Model," Technical Report 704, The Johns
Hopkins University, Dept. of Biostatistics.
Sest, M., Hamson, P. J., and Migon, H. S. (1985), "Dynamic Generalized
Linear Models and Bayesian Forecasting" (with discussion), Journal of
the American Statistical Association, 80, 73-97.
Whittaker, J. (1990), Graphical Models in Applied Multivariate Statistics,
New York: John Wiley.
Williams, D. A. (1982), "Extra-Binomial Variation in Logistic Linear Models," Applied Statistics, 3 1, 144-148.

25
Zeger, S. L. (1988), "A Regression Model for Time Series of Counts," Biometrika, 75, 62 1-629.
Zeger, S. L., and Karim, M. R. (1991), "Generalized Linear Models with
Random Effects; A Gibbs Sampling Approach," Journal ofthe American
Statistical Association, 86, 79-86.
Zeger, S. L., Liang, K. Y., and Albert, P. S. (1988), "Models for Longitudinal
Data: A Generalized Estimating Equation Approach," Biometrics, 44,1049- 1060.
Zelen, M. (197I), "The Analysis of Several Contingency Tables," Biometrika,
58, 129-137.