Bayesian and likelihood methods for fitting multilevel models with complex level-1 variation.
In multilevel modelling it is common practice to assume constant variance at level 1 across individuals. In this paper we consider situations where the level-1 variance depends on predictor variables. We examine two cases using a dataset from educational research; in the first case the variance at level 1 of a test score depends on a continuous "intake score" predictor, and in the second case the variance is assumed to differ according to gender. We contrast two maximum-likelihood methods based on iterative generalised least squares with two Markov chain Monte Carlo (MCMC) methods based on adaptive hybrid versions of the Metropolis-Hastings (MH) algorithm, and we use two simulation experiments to compare these four methods. We find that all four approaches have good repeated-sampling behaviour in the classes of models we simulate. We conclude by contrasting raw- and log-scale formulations of the level-1 variance function, and we find that adaptive MH sampling is considerably more efficient than adaptive rejection sampling when the heteroscedasticity is modelled polynomially on the log scale.
Introduction
Over the past 15 years or so, tting multilevel models to data with a hierarchical or
nested structure has become increasingly common for statisticians in many application
areas (e.g., Goldstein 1986, 1995; Bryk and Raudenbush 1992; Draper 2000). The main
purpose of tting such models is to partition the variation in a response variable as a
function of levels in the hierarchy and relate this variability to descriptions of the data
structure. In education, for example, multilevel modelling can be used to calculate the
proportion of variation in an observation that is explained by the variability between
students, classes, and schools in a 3{level nested structure. Random-eects modelling
of this kind is generally combined with xed-eects modelling, in which predictors are
additionally related to the response variable as covariates.
Generally these models assume a constant level{1 variance for the error or residual
term for all observations (in our notation students are at level 1 in the 3-level structure
above), but there is no reason why this should be true in all applications. An alternative
is to allow heteroscedasticity|in other words, to t models that relate the amount of
level{1 variability to predictor variables. We will refer to this here as complex level{1
variation. Heteroscedasticity is a common modelling concern in the standard tting of
Institute of Education, University of London, 20 Bedford Way, London WC1H 0AL, UK (email
bwjsmsr@ioe.ac.uk [WJB], h.goldstein@ioe.ac.uk [HG], and teuejra@ioe.ac.uk [JR]).
y Department of Mathematical Sciences, University of Bath, Claverton Down, Bath BA2 7AY, UK
(email d.draper@maths.bath.ac.uk, web http://www.bath.ac.uk/masdd).
W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash

Table

1: A comparison of means and variances of normalised exam scores for various
partitions of the GCSE dataset.
Partition Size Mean Variance
Whole dataset 4,059 0.000 1.000
Boys 1,623 0:140 1.052
Girls 2,436 0.093 0.940
Standardised LRT < 1 612 0:887 0.731
0:5 < Standardised LRT < 0:1 619 0:191 0.650
0:1 < Standardised LRT < 0:3 710 0.044 0.658
0:3 < Standardised LRT < 0:7 547 0.279 0.659
0:7 < Standardised LRT < 1:1 428 0.571 0.678
1:1 < Standardised LRT 549 0.963 0.703
linear models to data lacking a hierarchical or multilevel structure (e.g., Weisberg 1985),
but far less attention has been paid to this topic with multilevel data.
As our main motivating example we consider a dataset studied in Rasbash et al. (2000),
which was originally analysed in Goldstein et al. (1993). This dataset contains exam results
for 4,059 pupils from schools sampled from six inner London Education Author-
ities. The response variable of interest is the total score achieved in GCSE examinations
(a standardised test taken at age 16 by these pupils). This variable has already been
normalised (transformed by replacing each value by its standard normal score) in the
dataset we consider.

Table

contains mean and variance estimates for the response variable for various
partitions of the dataset. One of the main predictors of interest is a score on a reading
test (LRT) that all pupils took at age 11. For purposes of partitioning we have divided
the pupils into 7 groups of roughly equal sample size based on a standardised version of
the LRT score. From the mean column of the table it is clear that girls generally do a bit
better than boys and that the LRT score is positively correlated with the exam score. It
can also be seen that boys' exam scores are slightly more variable than girls' scores and
that the variance of the exam score bears a roughly quadratic relationship to LRT score.
Both of these conclusions mean that in tting a multilevel model to this dataset it will
be worth considering the need for complex variation at level 1.
The plan of the paper is as follows. In Section 2 we describe two variations on a
maximum-likelihood approach to the tting of multilevel models with complex level{1
variation and examine several examples of complex variance structures. Sections 3 and 4
present a Markov chain Monte Carlo (MCMC) method for Bayesian tting of such models
based on adaptive Metropolis-Hastings sampling, using two dierent proposal distribu-
tions. In Section 5 we give results from two simulation studies investigating the bias and
interval coverage properties, in repeated sampling, of the four tting methods described
in the previous three sections. Section 6 examines alternatives (a) to our MCMC methods
and (b) to our formulation of complex variance structures, and Section 7 discusses our
conclusions and suggests extensions of the work presented here.
Fitting multilevel models with complex level{1 variation 3
Maximum-likelihood-based methods and complex variance
structures
We begin by describing a general 2{level model with complex variation (later sections
will examine methods to t alternatives to this general model with additional constraints
added). The basic structure for a general Gaussian multilevel model is
y  N n (X; V
Here y is an (n  1) vector of responses, not necessarily independently distributed, with
a (p f  1) vector of xed-eect coe-cients of the predictors in the (n
n is the total number of level{1 observations in the data set (4,059 students, in the
example in Section 1), and p f is the number of xed eects in the model. The (n  n)
covariance matrix V for the responses contains all the random structure in the model; for
the two-level case we can write the variance term V
level{2 unit j. In this expression the variance has been partitioned into separate terms
for the two levels, with e and u denoting random eects at levels 1 and 2, respectively.
The covariances between the responses have the form V ij;i
if the two
observations are in the same level{2 unit, and V ij;i This means that
if the y vector is ordered so that all the observations in each level{2 unit are grouped
together, V has a block diagonal form.
In this general formulation the level{1 and level{2 variances and covariances are potentially
dierent for each pair of observations, but important special cases exist with
simpler structure, e.g., variance-components models where both the level{1 and level{
2 variances are constant across observations. The covariates X may make an appearance
in the random structure of the model, leading to a further partition of  u;ij . An
example is a random-slopes regression model with a single predictor X 1;ij , in which
u;ij
u;11 .
Here
u consists of the variance and co-variance
terms at level 2 expressed as a matrix (with structural zeroes where neces-
sary). Using this notation the (general) within-block covariance term can be written
a vector of predictors.
In the language of this section, what was referred to earlier as complex variation at
level 1 simply means partitioning the level{1 variance so that it depends in a natural way
on predictor variables. Figure 1 presents several potential variance structures that can
be tted to the GCSE dataset described earlier. The corresponding models are
e;ij
u;ij
e;ij
u;ij
e;ij
4 W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash
Standardised LRT score
Levelvariance
Model 2
Standardised LRT score
Variance
Level 1 variance
Level 2 variance
Standardised LRT score
Variance
Level 1 variance
Level 2 variance
Standardised LRT score
Variance
Model 5
Level 1 variance Boys
Level 1 variance Girls
Level 2 variance

Figure

1: Four dierent variance structures tted to the GCSE dataset.
u;ij
e;ij
In all these models X 1 refers to the standardised LRT score and X 2 refers to gender
(coded 0 for boys and 1 for girls). In equation (2) we have a simple one-level regression
model with a quadratic variance relationship with LRT; the other models involve tting
increasingly complex variance structures to the data in a two-level framework.
One approach to tting models such as (2{5) via maximum likelihood (ML) is based on
iterative generalised least squares (IGLS), and its restricted variant (RIGLS, also known
as REML) which corrects for bias. The basic idea is similar to that of the EM algorithm
in that (a) an estimate ^
of  in (1) is obtained using a current estimate of V and (b) an
estimate of V is then obtained using the ^
from (a), but in IGLS/RIGLS the estimation
Fitting multilevel models with complex level{1 variation 5

Table

2: IGLS estimates for models (2{5) tted to the GCSE dataset (standard errors
(SEs) in parentheses).
Model
Parameter (2) (3) (4) (5)
u;00 | 0.094 (0.018) 0.091 (0.018) 0.086 (0.017)
u;01 | | 0.019 (0.007) 0.020 (0.007)
u;11 | | 0.014 (0.004) 0.015 (0.004)
e;12 | | | 0.032 (0.013)
e;22 | | | 0:058 (0.026)
of the covariance matrix V is recast as a regression problem and weighted least squares
is used in both steps (see Goldstein 1986, 1989 for details).

Table

estimates obtained for models (2{5) applied to the GCSE data.
Both gender and LRT score are evidently useful in predicting GCSE score. Model (2),
which naively ignores the hierarchical nature of the data, hints at heteroscedasticity (the
ML estimate
of
e;11 is about as big as its standard error (SE)), but (from the estimates
of
u;00 in equations (3{5)) there is a clear need for two-level modelling, and the full
complexity of what is required to describe the data only comes into focus with model (5)
(in which every estimate is at least 2.2 times as large as its SE).
3 An MCMC method for a general 2{level Gaussian model
with complex level{1 variation
Browne and Draper (2000a and 2000b) gave Gibbs-sampling algorithms for Bayesian
tting of 2{level variance-components and random-slopes-regression models, respectively.
In this section we consider a general 2{level model with complex variation at level 1; this
can easily be generalised to an N{level model via an approach similar to the method
detailed in Browne (1998). For MCMC tting of model (1) it is useful to rewrite it as
follows, with y ij denoting the (scalar) outcome for (level{1) observation i in level{2 unit
(0;
(0;
e
are (p f  1); (p 2  1), and (p 1  1) vectors of xed-eects parameters
and level{2 and level{1 residuals,
are vectors of predictor
values; and p 1 and p 2 are the numbers of parameters specifying the random eects at
levels 1 and 2, respectively. The IGLS/RIGLS methods do not directly estimate the
but they can be estimated after tting the model using a method given in
6 W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash
Goldstein (1995). In equation
(6),
e
and
u are the variance terms at level 1 and level
Gibbs sampling procedures for tting multilevel models such as (6) proceed most
smoothly by treating the level{2 residuals as latent variables when forming the full conditional
posterior distributions. In a multilevel model with simple (homoscedastic) variation
at level 1, the level{1 residuals may be calculated at each iteration by subtraction. In
the above model we cannot explicitly compute the individual level{1 residuals; instead we
deal with the \composite" residuals X C
these can be calculated by subtraction. The
important part of the algorithm that follows is to store the composite level{1 variance
function for each individual,
All the other parameters then depend on the level{1 covariance
matrix
e through these
individual variances. This means that the algorithm that follows, apart from the updating
step
for
e , is almost identical to the algorithm for the same model without complex
variation (Browne 1998).
3.1 Inverse-Wishart proposals for the level{1 covariance matrix
In the rst MCMC method examined in this paper, we collect together the terms in
the variance equation at level 1,  e;ij , into the covariance
matrix
e .
Updating
e
using a Metropolis-Hastings (MH) algorithm therefore requires a proposal distribution
that generates positive-denite matrices (later we will relax this restriction). We use
an inverse-Wishart proposal distribution with expectation the current
estimate
e at
iteration t to
generate
e . In the parameterisation used, for example, by Gelman,
Carlin, et al. (1995), the inverse-Wishart distribution W 1
expectation
e , where w is a positive
integer degrees of freedom parameter, this will produce a distribution with expectation
e . The parameter w is a tuning constant which may be set to an (integer) value that
gives the desired MH acceptance rate.
For prior distributions on the parameters in model (6) we make the following choices
in this algorithm: a generic prior
e ) (to be specied in Section 4.2) for the level{1
covariance matrix, an inverse-Wishart
prior
for the level{2 covariance
matrix, and a multivariate normal prior   N p f
for the xed eects parameter
vector. The algorithm, which is detailed in the Appendix, is a hybrid of Gibbs and MH
steps; it divides the parameters and latent variables in (6) into four blocks and uses
multivariate normal Gibbs updates for  and the u j , inverse-Wishart Gibbs updates for
u , and inverse-Wishart MH proposals
for
e .
3.2 An adaptive method for choosing the tuning constant w
Browne and Draper (2000b) describe an adaptive hybrid Metropolis-Gibbs sampler for
tting random-eects logistic regression models. Gibbs sampling may be used in such
models for variance parameters, but Metropolis updates are needed for xed eects and
latent residuals. Browne and Draper employ a series of univariate normal proposal distributions
(PDs) for these quantities, and give a procedure for adaptive choice of appropriate
Fitting multilevel models with complex level{1 variation 7

Table

3: An illustration of the adaptive MH procedure with model (4) applied to the GCSE
data.
Acceptance Within
Iterations Rate w Tolerance?
100 20% 138 0
200 19% 195 0
300 30% 208 1
500 30% 229 3
values for the variances of these PDs to achieve e-cient MH acceptance rates. Here we
provide a modication of this procedure for the case of inverse-Wishart proposals.
We set the tuning parameter w described above to an arbitrary starting value (in
the example that follows 100) and run the algorithm in batches of 100 iterations. The
goal is to achieve an acceptance rate for the level{1 covariance matrix that lies within
a specied tolerance interval (r We compare the empirical acceptance rate
r  for the current batch of 100 iterations with the tolerance interval, and modify the
proposal distribution appropriately before proceeding with the next batch of 100. The
modication performed at the end of each batch is as follows:
If r
r
where only the integer part of w is used in (8). The amount by which w is altered in each
iteration of this procedure is an increasing function of the distance between r and r  . The
adaptive procedure ends when three successive r  values lie within the tolerance interval;
the value of w is then xed and we proceed with the usual burn-in and monitoring periods.
3.3 An example
We consider the model in Section 2 which has a quadratic relationship between the variance
and the LRT predictor (model (4)). The adaptive procedure was run for this model
with a target acceptance rate of (based on a recommendation in Gelman, Roberts
and Gilks 1995) and a tolerance of summarises the progress of the adaptive
method in this example; here only 500 iterations are required to adjust the proposal
distribution to give the desired acceptance rate (500{2000 iterations are typically needed
in the applications we have examined).

Table

4 compares the estimates produced by this MCMC method for model (4) to
those (a) from the IGLS and RIGLS procedures and (b) from another MCMC method to
be described in the next section (here and throughout the paper, MCMC point estimates
are posterior means). We used a slightly informative inverse-Wishart prior for the level{2
covariance matrix for the MCMC methods based on the RIGLS estimate, and a uniform
prior for the level{1 covariance matrix. In this case the results for all the methods are
fairly similar, with one exception: the
parameter
e;11 is noticeably larger using MCMC
8 W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash

Table

4: Parameter estimates for four methods of tting model (4) to the London schools
dataset (SEs/posterior standard deviations in parentheses). The MCMC methods were
monitored for 50,000 iterations after the adaptive procedure and a burn-in of 500 iterations

MCMC Method
Parameter IGLS RIGLS 1 2
method 1. This dierence highlights the fact that the rst MCMC approach actually ts
a model with an extra positive-denite constraint: we are
forcing
e;11 to be positive,
which in
ates the point estimate. The second MCMC method, which we consider below,
is based on dierent constraints; when we examined the chain of values it produced for
e;11 we found that nearly 40% of the values were negative. There is no inconsistency in
this result: in the model to be examined in the next
section,
e;11 is not a variance.
4 Truncated normal proposals for the level{1 variance func-
tion
The inverse-Wishart updating method assumes that the variance function at level 1 arises
from a positive-denite covariance matrix. We now consider an alternative method that,
in a manner similar to IGLS and RIGLS, only requires the variance at level 1 to be a
linear function of the parameters. This MCMC solution will still have more constraints
than the IGLS solution, because we are still considering the level{1 and level{2 variances
separately and both of these quantities must be positive.
The constraint used in MCMC method 1 that the covariance matrix at level 1 is
positive-denite is actually stronger than necessary. Positive-denite matrices will guarantee
that any vector X C
ij will produce a positive variance in equation (6); a milder but still
scientically reasonable constraint is to allow all values
of
e such that
for all i and j. This restriction appears complicated to work with, but if we consider each
of the parameters
in
e separately and assume the other variables are xed the constraint
becomes manageable. It is once again useful to rewrite model (1), this time as follows:
e
(0;
where e
e;ij is given by equation (7). Here the composite level{1 residuals
e
are normally distributed with variances that depend on the predictors; consequently
Fitting multilevel models with complex level{1 variation 9
the constraint that the level{1 variance is always positive is still satised
but
e need not
be positive-denite.
4.1 MH updating: method 2
Our second method is identical to the rst for , the u j ,
and
in the
Appendix) but involves a Hastings update with a dierent proposal distribution
for
e .
We update each parameter in the level{1 variance equation in turn, always requiring for
all i and j at every iteration t in the Markov chain that
0: (10)
Considering rst the diagonal
terms,
e;kk , for each can be
d C
here X C
is the kth element of the vector X C
ij . This is equivalent to requiring that
d C
We use a normal proposal distribution with variance s 2
kk but reject generated values
that fail to satisfy (12). This amounts to using a truncated normal proposal, as shown in

Figure

2(i). The Hastings ratio R can then be calculated as the ratio of the two truncated
normal distributions shown in Figure 2(i) and (ii). Letting the value
for
e;kk at time t
be A and the proposed value for time
s kk
s kk
The update step is then as
follows:
e;kk with probability min4
e;kk otherwise
and the corresponding density in the denominator of (14) are
given by (28).
The diagonal terms are a special case as they are always multiplied by a positive quantity
in the variance equation, so that the proposal distribution needs only one truncation
point. More generally for the non-diagonal
terms
e;kl we get the following. As before,
at time t for all i and j constraint (10) must be satised; for each 1  k < l  p 1 this can
be rewritten
e;kl d C
d C
e;kl
W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash
(iv)

Figure

2: Plots of truncated univariate normal proposal distributions for a parameter .
A is the current value  c and B is the proposed new value   . M is max  and m is
min  , the truncation points. The distributions in (i) and (iii) have mean  c , while the
distributions in (ii) and (iv) have mean   .
This is equivalent to the two constraints
over all (i; j) such that X C
e;kl < min e;kl  min ij
over all (i; j) such that X C
We again use a normal proposal distribution, this time with variance s 2
kl , and again values
failing to satisfy (16) are rejected. This leads to the truncated normal proposal shown in

Figure

2(iii). The Hastings ratio R is then simply the ratio of the two truncated normal
distributions shown in Figure 2(iii) and (iv). Letting the value
for
e;kl at time t be A
Fitting multilevel models with complex level{1 variation 11
and the proposed value for time
s kl
s kl
min e;kl A
s kl
s kl
The update step is then similar to (14) with subscripts kl in place of kk in
the
e terms.
4.2 Proposal distribution variances and prior distributions
In the method outlined above we consider each parameter
in
e separately. This means
that we use a separate truncated univariate normal proposal distribution for each parameter
subject to the constraints that the value generated will produce a positive level{1
variance  e;ij for all i and j. We therefore need to choose a proposal distribution variance
for each parameter. Two possible solutions are to use the variance of the parameter
estimate from the RIGLS procedure multiplied by a suitable positive scale factor, or to
use an adaptive approach before the burn-in and monitoring run of the simulation. See
Browne and Draper (2000a) for a description of both of these methods in the case of
random eects logistic regression models.
Prior distributions using this method must take account of the constraints imposed on
the parameters. In all the analyses we perform in this paper with this method, we use a
series of marginal uniform priors for the level{1 variance terms subject to the constraints;
in other words, all valid combinations of parameter estimates
for
e are a priori equally
likely. Other prior distributions may be problematic.
4.3 Examples
Model (4) was tted to the GCSE data in Section 3.3, and the estimates produced by
both MCMC methods are shown in Table 4. For the truncated normal method we used
the adaptive MH procedure, in this case with a desired acceptance rate of 50% as the
parameters are updated separately. The advantage of the truncated normal method is
that it can handle variance functions that would not necessarily have a positive-denite
matrix form. For illustration we consider a simple case which the inverse-Wishart method
cannot t. Our model is as follows:
This model includes a variance for boys and a term that represents the dierence in
variance between boys and girls. The results from tting this model are given in Table 5;
all methods give roughly the same estimates for the level{1 variance terms. The total
variances produced by the model for
are similar to the values given in the part of Table 1 where the variance in the response
is calculated for boys and girls separately.
12 W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash

Table

5: Parameter estimates for three methods tted to model (18) for the GCSE dataset.
method 2, using truncated normal proposals, was monitored for 50,000 iterations
following the adapting period and a burn-in of 500 iterations.
MCMC
Parameter IGLS RIGLS Method 2
5 Simulation studies
In this section we examine the bias and interval-coverage properties, in repeated sampling,
of the four methods described above, in two sets of simulated models with complex level{
variation based on the GCSE example. We rst consider model (4), which features
a quadratic variance relationship with the input reading test (LRT) predictor. As true
(population) parameters for our simulation we used values close to the estimates obtained
in the actual data, with one exception: we
increased
e;11 so that the correlation of the
random eects at level 1 was reduced. This is because sample datasets drawn from
multilevel models with high correlation cause convergence problems with the IGLS and
RIGLS methods (Browne and Draper 2000b).
One thousand datasets were generated randomly according to model (4)|with the
same numbers of level{1 and level{2 units (4,069 and 65, respectively) as in the original
GCSE data set, and the same distribution of level{1 observations within level{2
units|and tted using the four methods, with the results presented in Table 6. For the
MCMC methods (in both of the simulation studies) the posterior distribution with each
dataset was monitored for 10,000 iterations after the adapting period and a burn-in of
500 from IGLS starting values. Uniform priors were used for the level{1 variances and
xed eects. A (slightly) informative inverse-Wishart prior was used for the level{2 co-variance
matrix in line with the results in Browne and Draper (2000b). Interval estimates
at nominal level 100(1 )% with the IGLS and RIGLS approaches were of the form
based on the large-sample normal approximation; this is what
users of most multilevel packages such as MLwiN (Rasbash et al. 2000) and HLM (Bryk
et al. 1988) would report, if they report interval estimates at all (since such packages
routinely only report estimates and standard errors with maximum-likelihood methods).
With the Bayesian MCMC methods we give results based on posterior means as point
estimates and 90%/95% central posterior intervals.
Our second simulation study (Table 7) was based on model (18) from Section 4.3, in
which the male and female subsamples had dierent level{1 variances (MCMC method 1 is
not available for this model). We again created 1,000 simulation datasets with population
values similar to the estimates obtained with the GCSE dataset. With the Bayesian
approach to tting, uniform priors were used for the level{1 variances and xed eects,
Fitting multilevel models with complex level{1 variation 13

Table

Summary of results for the rst simulation study, with LRT score random at
levels 1 and 2. Bias results in (a) are relative except those in brackets, which are absolute
(the true value in those cases is zero). Monte Carlo standard errors (SEs) in (a) are
given in parentheses. The Monte Carlo SEs for the estimated interval coverages in (b)
range from 0.7% to 1.0%.
(a) Relative bias of point estimates (%)
Parameter MCMC MCMC
fTrue Valueg IGLS RIGLS Method 1 Method 2
u;00 f0:1g 1:81 (0.59) 0:08 (0.60) 2.79 (0.62) 2.79 (0.62)
(b) Interval coverage probabilities at nominal levels 90%/95%
MCMC MCMC
Parameter IGLS RIGLS Method 1 Method 2
u;00 89.4/93.1 90.7/93.6 91.1/96.0 91.1/96.0
u;01 90.0/94.4 90.3/94.6 88.7/94.1 88.8/94.1
e;00 90.7/94.1 90.7/94.1 90.2/94.1 90.9/94.8
e;11 90.6/95.1 90.7/95.1 90.0/95.0 90.9/95.4
(c) Mean interval widths at nominal levels 90%/95%
MCMC MCMC
Parameter IGLS RIGLS Method 1 Method 2
14 W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash

Table

7: Summary of results for the second simulation study, with separate variances
at level 1 for boys and girls. Monte Carlo standard errors (SEs) in (a) are given in
parentheses. The Monte Carlo SEs for the estimated interval coverages in (b) range from
0.7% to 1.0%.
(a) Relative bias of point estimates (%)
Parameter MCMC
fTrue Valueg IGLS RIGLS Method 2
(b) Interval coverage probabilities at nominal levels 90%/95%
MCMC
Parameter IGLS RIGLS Method 2
u00 88.3/92.5 89.2/93.4 90.0/95.3
(c) Mean interval widths at nominal levels 90%/95%
MCMC
Parameter IGLS RIGLS Method 2
Fitting multilevel models with complex level{1 variation 15
and a 1 (; ) prior (with used for the level{2 variance
parameter
in line with the results in Browne and Draper (2000a).
It is evident from Tables 6 and 7 that all four methods performed reasonably well in
both models. RIGLS succeeded in reducing the (already small) biases arising from IGLS
estimation in most cases, and the relative biases of the MCMC methods are also small
(ranging from 0% to 4.8%, with a median absolute value of 1.5%). Interval coverages
for all four methods were all close to nominal, with actual coverages ranging from 86{
91% and 91{96% at nominal 90% and 95%, respectively; all four methods achieved this
level of coverage with intervals of comparable length; and the ratios of 95% and 90%
interval lengths for each method were all close to the value (  1 (0:975)
to be expected
under normality. The ML methods have the clear advantage of speed (on the original
GCSE data set IGLS/RIGLS and MCMC methods 1 and 2 took 2, 168, and 248 seconds
on a 500{MHz Pentium PC, respectively, with the MCMC methods based on 10,000
monitoring iterations), but the ML approach has two potential disadvantages: on data sets
with small numbers of level{1 and level{2 units, it requires more sophisticated methods
for constructing interval estimates for variance parameters (to achieve good coverage
properties) than the large-sample normal approximation used here (Browne and Draper
2000ab), and it may fail to converge when
the
e
and/or
matrices exhibit a high degree
of correlation between the parameters quantifying the random eects. The Bayesian
methods are considerably slower but have the additional advantage that inferences about
arbitrary functions of the model parameters are automatic once the model parameters
themselves have been monitored.
6 Other MCMC methods
6.1 Gibbs Sampling
There are special cases of the problem of complex level{1 variation that can be tted
using a standard Gibbs sampler. The model (equation 18) used in the second simulation,
where we use a dierent level{1 variance term for each gender, is one such example. Here
we could reparameterise the model with two variances, one for boys ( 2
b ) and one for girls
rather than a boys' variance plus a dierence. Scaled-inverse- 2 priors (see, e.g.,
Gelman et al. 1995) can be used for these two variances, with parameters
respectively. If we divide the children into boys' and girls' subgroups B and G,
of size n b and n g , then step 3 of the algorithm given in the Appendix can be rewritten as
two Gibbs sampling steps as follows: the full conditional for  2
b is
a
and the full conditional for  2
g is exactly analogous. Now the level{1 variance is
and the other steps of the algorithm are as before.
6.2 Modelling the variance on the log scale
The developers of the software package BUGS (Spiegelhalter et al. 1997) use a dierent
approach to tting complex level{1 variation in one of their examples, the Schools data
W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash
set (example 9 in Volume 2 of Spiegelhalter et al. 1996). They model the logarithm of
the level{1 precision as a function of predictors and other parameters:
(0;
This results in a multiplicative, rather than an additive, variance function:
exp
The advantages of this approach are that the parameters are now unconstrained, the level{
variance will never be negative, and it is easier to specify a prior with this method. The
disadvantages are that the interpretation of the individual coe-cients is not as easy and
computation for these models is slower. The interpretation di-culty will be apparent
mainly when the X variables are categorical.
Model (20) can be tted in BUGS using adaptive rejection (AR) sampling (Gilks and
Wild 1992). Alternatively the adaptive MH method used in the truncated normal algorithm
in Section 4.1 can be used, this time with no parameter constraints and hence no
truncation in the normal proposal distributions. Goldstein (1995, Appendix 5.1) shows
how to obtain ML estimates for this model; see Yang et al. (2000) for a set of MLwiN
macros to do this.
To explore the dierences between log-variance modelling and our earlier approach,
we tted four dierent level{1 variance functions to the GCSE dataset to model the eect
of LRT score (X 1 ) on the level{1 variance. We considered the quadratic relationship
examined earlier (model (4)), and the simpler linear relationship
e;ij
we also considered two exponential relationships:
In each of the four models the level{2 variance structure and xed eects were as in
equation (4). Figure 3 plots the resulting level{1 estimated variances as a function of
LRT score. In this case for the majority of the data the variance estimates produced by
all four models are fairly similar, with any discrepancies between the models occurring at
the extremes of the LRT range where there are relatively few observations.

Table

presents estimates of the four models (using MH method 2), together with
Raftery-Lewis (1992) default diagnostics (for both MH and AR sampling in the exponential
models (23)) and comparative timings. From part (b) of the table it is evident
that the parameter with the worst MCMC mixing is the intercept  0 . This means that,
although the MH method requires longer monitoring runs than the AR approach for the
Fitting multilevel models with complex level{1 variation 17
Standardised LRT score
LevelVariance
Quadratic
Exp. Linear
Exp. Quadratic

Figure

3: Four ways to model the eect of standardised LRT score on the level{1 variance
in the GCSE dataset.
level{1 variance parameters, the run lengths required to ensure that all parameter estimates
have a specied accuracy (with respect to 95% interval estimation) will be roughly
equal. From part (c) of the table it can be seen that the MH approach is 4{9 times faster
in real-time execution speed in this example. Results in Table 8 are based on a single
data set but are typical of ndings we have obtained with other similar models.
7 Conclusions and extensions
In this paper we have presented several methods for modelling non-constant level{1 variance
functions with multilevel data. We have introduced two new adaptive Metropolis-Hastings
sampling methods for tting such functions subject to dierent constraints. The
two methods give similar estimates for models where the true parameter values are not
aected by the constraints, but if the true values do not satisfy the additional positive-
denite matrix constraint of the inverse-Wishart proposal method then the estimates from
the two methods will dier.
The main advantage of the inverse-Wishart method is that it models the level{1
variance function as a matrix in a manner analogous to the usual treatment of the level{2
variance function, meaning (among other things) that informative inverse-Wishart priors
at level 1 can be used with this approach. The main advantage of the truncated normal
proposal method is that it is more general and can deal with any variance function at level
W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash

Table

8: Parameter estimates for four dierent level{1 variance functions applied to
the GCSE dataset and t by MCMC. All methods used a monitoring period of 50,000
iterations after a burn-in of 500 iterations. Methods with an adaptive MH step at level
were run using a development version of MLwiN; those with an adaptive rejection (AR)
step at level 1 were run in WinBUGS. Posterior standard deviations are given in parentheses
in (a).
(a) Parameter estimates
Exponential Exponential
Parameter Linear Quadratic Linear Quadratic
e;11 | 0.003 (0.009) | 0:0005 (0.016)
(b) Raftery-Lewis values (in thousands of iterations); main entries
apply to MH method 2, with the corresponding values for AR in parentheses
Exponential Exponential
N Linear Quadratic Linear Quadratic
u;00 4.3 4.3 4.4 (4.3) 4.3 (4.1)
e;11 | 16.9 | 15.6 (4.6)
(c) Timings (in minutes at 500 Pentium MHz)
Exponential Exponential
Method Linear Quadratic Linear Quadratic
Metropolis-Hastings 19 20 23 25
Adaptive Rejection | | 95 227
Fitting multilevel models with complex level{1 variation 19
1. Both methods have bias and interval coverage properties that are similar to those from
the maximum-likelihood (IGLS and RIGLS) approaches, and all four methods perform
satisfactorily in repeated sampling in this regard.
In Section 6.2 we considered an alternative formulation of the level{1 variance function
in terms of the log of the precision at level 1. This method has two advantages: there is
no need to impose constraints on the terms in the resulting variance function, and it is
therefore easier to contemplate a variety of prior distributions for the resulting variance
parameters. The main disadvantage of this approach is that the individual terms in
the variance function may not be as easily interpreted, making it potentially di-cult
to construct sensible informative priors. Table 8 shows clearly, however, that adaptive-
rejection sampling is much less e-cient than adaptive Metropolis-Hastings sampling to
achieve default MCMC accuracy standards with variance (or precision) functions that are
exponential in the parameters.
There are two obvious extensions of this work, to arbitrary variance structures at
higher levels and to multivariate normal responses. Two approaches to tting random
eects at level 2 and above appear common in current applied work: modelling all random
eects independently, or tting fully dependent random eects with a complete covariance
matrix at each level (see the Birats example in Spiegelhalter et al. (1996) for an illustration
of both formulations). It is fairly easy to t any block-diagonal covariance structure at a
higher level using Gibbs sampling, in a straightforward extension of the approach given
in Section 6.1. The adaptive MH sampler with a truncated normal proposal (method 2,
Section 4.1) can be used to t any dependence structure among the random eects at the
higher levels, including non block-diagonal covariance matrices.
With multivariate normal response models the variance function at the lowest level
includes variances for each response plus covariances between responses. This variance
function could also be extended to include predictors that may in
uence the variance
of individual responses in an analogous way to the univariate model. We intend to report
on MCMC sampling algorithms for general multivariate-response multilevel models
elsewhere.

Acknowledgements

We are grateful to the EPSRC, ESRC, and European Commission for nancial support,
and to David Spiegelhalter, Nicky Best, and other participants in the BUGS project for
references and comments on the set of multilevel modelling papers based on the rst
author's PhD dissertation. Membership on this list does not imply agreement with the
ideas expressed here, nor are any of these people responsible for any errors that may be
present.


Appendix

: Details of MCMC method 1
In step 1 of the algorithm described in Section 3.1, the full conditional distribution in the
Gibbs update for the xed eects parameter vector  is multivariate normal: with p f the
W. J. Browne, D. Draper, H. Goldstein, and J. Rasbash
number of xed eects,
e;ij
and
e;ij
involves a Gibbs update of the level{2 residuals, u j , also with a multivariate normal
full conditional distribution: with p 2 the number of parameters describing the random
eects at level{2 and n
j the number of level{1 units in level{2 unit j,
Pn
e;ij
and
Pn
e;ij
Step 3 employs a Hastings update using an inverse-Wishart proposal distribution for the
level{1 covariance
matrix
e . Specically, the Markov chain moves
from
e at time
to
e as
follows:
e with probability min
e
e
e otherwise
Here
(a)
e  W 1
e
, where w is chosen as in Section 3.2 and p 1 is the
number of rows or columns
in
(b) the Hastings ratio R in (26) is

exp
tr
e
e
tr
e
(c) the full conditional distribution
for
e in (26) is
Y
1e;ij exp
where we have expressed the right-hand side of (28) for convenience in terms of  e;ij as in
equation (7).
Finally, step 4 involves a Gibbs update of the level{2 covariance
matrix
expressed as a
ofu the full conditional is
e
is the number of rows or columns
in
u and J is the number of level{2 units in the
data set. An improper uniform prior
on
u corresponds to the choice
Fitting multilevel models with complex level{1 variation 21



--R

Applying MCMC Methods to Multilevel Models.


Bryk AS
Bryk AS
Bayesian Hierarchical Modeling.
Bayesian Data Analysis.

Adaptive rejection sampling for Gibbs sampling.
Multilevel mixed linear model analysis using iterative generalised least squares.
Restricted unbiased iterative generalised least squares estimation.
Multilevel Statistical Models (second edition).
A multilevel analysis of school examination results.
How many iterations in the Gibbs sampler?
A User's Guide to MLwiN (Version 2.1).
BUGS 0.5 Examples (Version ii).
Cambridge: Medical Research Council Biostatistics Unit.
BUGS: Bayesian Inference Using Gibbs Sampling (Version 0.60).
Applied Linear Regression
MLwiN Macros for Advanced Multilevel Modelling (Version 2.0).
--TR
maximum-likelihood methods;adaptive Metropolis-Hastings sampling;hierarchical modelling;educational data;IGLS;MCMC efficiency;Markov chain Monte Carlo MCMC;multilevel modelling;heteroscedasticity;RIGLS