=Paper=
{{Paper
|id=Vol-1218/bmaw2014_paper_1
|storemode=property
|title=A Comparison of Two MCMC Algorithms for Hierarchical Mixture Models
|pdfUrl=https://ceur-ws.org/Vol-1218/bmaw2014_paper_1.pdf
|volume=Vol-1218
|dblpUrl=https://dblp.org/rec/conf/uai/Almond14
}}
==A Comparison of Two MCMC Algorithms for Hierarchical Mixture Models==
https://ceur-ws.org/Vol-1218/bmaw2014_paper_1.pdf
A Comparison of Two MCMC Algorithms for Hierarchical Mixture Models
Russell G. Almond⇤
Florida State University
Abstract 1 Introduction
Mixture models (McLachlan & Peel, 2000) are a fre-
Mixture models form an important class of
quently used method of unsupervised learning. They
models for unsupervised learning, allowing
sort data points into clusters based on just their val-
data points to be assigned labels based on
ues. One of the most frequently used mixture models
their values. However, standard mixture
is a mixture of normal distributions. Often the mean
models procedures do not deal well with rare
and variance of each cluster is learned along with the
components. For example, pause times in
classification of each data point.
student essays have di↵erent lengths depend-
ing on what cognitive processes a student As an example, Almond, Deane, Quinlan, Wagner, and
engages in during the pause. However, in- Sydorenko (2012) fit a mixture of lognormal distribu-
stances of student planning (and hence very tions to the pause time of students typing essays as
long pauses) are rare, and thus it is dif- part of a pilot writing assessment. (Alternatively, this
ficult to estimate those parameters from a model can be described as a mixture of normals fit to
single student’s essays. A hierarchical mix- the log pause times.) Almond et al. found that mix-
ture model eliminates some of those prob- ture models seems to fit the data fairly well. The mix-
lems, by pooling data across several of the ture components could correspond to di↵erent cogni-
higher level units (in the example students) tive process used in writing (Deane, 2012) where each
to estimate parameters of the mixture com- cognitive process takes di↵erent amounts of time (i.e.,
ponents. One way to estimate the parame- students pause longer when planning, than when sim-
ters of a hierarchical mixture model is to use ply typing).
MCMC. But these models have several issues
Mixture models are difficult to fit because they dis-
such as non-identifiability under label switch-
play a number of pathologies. One problem is com-
ing that make them difficult to estimate just
ponent identification. Simply swapping the labels of
using o↵-the-shelf MCMC tools. This paper
Components 1 and 2 produces a model with identical
looks at the steps necessary to estimate these
likelihoods. Even if a prior distribution is placed on
models using two popular MCMC packages:
the mixture component parameters, the posterior is
JAGS (random walk Metropolis algorithm)
multimodal. Second, it is easy to get a pathological
and Stan (Hamiltonian Monte Carlo). JAGS,
solution in which a mixture component consists of a
Stan and R code to estimate the models and
single point. These solutions are not desirable, and
model fit statistics are published along with
some estimation tools constrain the minimum size of
the paper.
a mixture component (Gruen & Leisch, 2008). Fur-
thermore, if a separate variance is to be estimated for
each mixture component, several data points must be
assigned to each component in order for the variance
estimate to have a reasonable standard error. There-
Key words: Mixture Models, Markov Chain Monte
fore, fitting a model with rare components requires a
Carlo, JAGS, Stan, WAIC
large data size.
⇤
Paper presented at the Bayesian Application Workshop at Un-
certainty in Artificial Intelligence Conference 2014, Quebec City,
Almond et al. (2012) noted these limitations in
Canada. their conclusions. First, events corresponding to the
1
highest-level planning components in Deane (2012)’s
cognitive model would be relatively rare, and hence
would lumped in with other mixture components due
to size restrictions. Second, some linguistic contexts
(e.g., between Sentence pauses) were rare enough that
fitting a separate mixture model to each student would
not be feasible.
One work-around is a hierarchical mixture model. As
with all hierarchical models, it requires units at two
di↵erent levels (in this example, students or essays
are Level 2 and individual pauses are Level 1). The
assumption behind the hierarchical mixture model is
that the mixture components will look similar across
the second level units. Thus, the mean and variance of
Mixture Component 1 for Student 1 will look similar
Figure 1: Non-hierarchical Mixture Model
to those for Student 2. Li (2013) tried this on some of
the writing data.
tion 2.4 describes prior distributions and constraints
One problem that frequently arises in estimating mix-
on parameters which keep the Markov chain away from
ture models is determining how many mixture compo-
those points.
nents to use. What is commonly done is to estimate
models for K = 2, 3, 4, . . . up to some small maximum
number of components (depending on the size of the 2.1 Mixture of Normals
data). Then a measure of model–data fit, such as AIC,
Consider a collection observations, Yi =
DIC or WAIC (see Gelman et al., 2013, Chapter 7), is
(Yi,1 , . . . , Yi,Ji ) for a single student, i. Assume
calculated for each model and the model with the best
that the process that generated these data is a mix-
fit index is chosen. These methods look at the deviance
ture of K normals. Let Zi,j ⇠ cat(⇡ i ) be a categorical
(minus twice the log likelihood of the data) and adjust
latent index variable indicating which component Ob-
it with a penalty for model complexity. Both DIC and ⇤
servation j comes from and let Yi,j,k ⇠ N (µi,k , i,k )
WAIC require Markov chain Monte Carlo (MCMC) to
be the value of Yi,j which would be realized when
compute, and require some custom coding for mixture
Zi,j = k.
models because of the component identification issue.
Figure 1 shows this model graphically. The plates in-
This paper is a tutorial for replicating the method used
dicate replication over categories (k), Level 1 (pauses,
by Li (2013). The paper walks through a script writ-
j) and Level 2 (students, i) units. Note that there is
ten in the R language (R Core Team, 2014) which per-
no connection across plates, so the model fit to each
forms most of the steps. The actual estimation is done
Level 2 unit is independent of all the other Level 2
using MCMC using either Stan (Stan Development
units. This is what Gelman et al. (2013) call the no
Team, 2013) or JAGS (Plummer, 2012). The R scripts
pooling case.
along with the Stan and JAGS models and some sam-
ple data are available at http://pluto.coe.fsu.edu/ The latent variables Z and Y ⇤ can be removed from
mcmc-hierMM/. the likelihood for Yi,j by summing over the possible
values of Z. The likelihood for one student’s data, Yi ,
2 Mixture Models is
Ji X
Y K
Yi,j µi,k
Let i 2 {1, . . . , I} be a set of indexes over the sec- Li (Yi |⇡ i , µi , i) = ⇡i,k ( ) (1)
i,k
ond level units (students in the example) and let j=1 k=1
j 2 {1, . . . , Ji } be the first level units (pause events
where (·) is the unit normal density.
in the example). A hierarchical mixture model is by
adding a Level 2 (across student) distribution over the Although conceptually simple, there are a number of
parameters of the Level 1 (within student) mixture issues that mixture models can have. The first is-
model. Section 2.1 describes the base Level 1 mixture sue is that the component labels cannot be identified
model, and Section 2.2 describes the Level 2 model. from data. Consider the case with two components.
Often MCMC requires reparameterization to achieve The model created by swapping the labels for Compo-
better mixing (Section 2.3). Also, there are certain pa- nents 1 and 2 with new parameters ⇡ 0i = (⇡i,2 , ⇡i,1 ),
rameter values which result in infinite likelihoods. Sec- µ0i = (µi,2 , µi,1 ), and 0i = ( i,2 , i,1 ) has an identical
2
likelihood. For the K component model, any permu-
tation of the component labels produces a model with
identical likelihood. Consequently, the likelihood sur-
face is multimodal.
A common solution to the problem is to identify the
components by placing an ordering constraint on one
of the three parameter vectors: ⇡, µ or . Section 4.1
returns to this issue in practice.
A second issue involves degenerate solutions which
contain only a single data point. If a mixture com-
ponent has only a single data point, its standard devi-
ation, i,k will go to zero, and ⇡i,k will approach 1/Ji ,
which is close to zero for large Level 1 data sets. Note
that if either ⇡i,k ! 0 or i,k > 0 for any k, then the
likelihood will become singular. Figure 2: Hierarchical Mixture Model
Estimating i,k requires a minimum number of data
points from Component k (⇡i,k Ji > 5 is a rough mini- ing mixture model. To make the hierarchical model,
mum). If ⇡i,k is believed to be small for some k, then a add across-student prior distributions for the student-
large (Level 1) sample is needed. As K increases, the specific parameters parameters, ⇡ i , µi and i . Be-
smallest value of ⇡i,k becomes smaller so the minimum cause JAGS parameterizes the normal distribution
sample size increases. with precisions (reciprocal variances) rather than stan-
dard deviations, ⌧ i are substituted for the standard
2.2 Hierarchical mixtures deviations i . Figure 2 shows the hierarchical model.
Gelman et al. (2013) explains the concept of a hier- Completing the model requires specifying distributions
archical model using a SAT coaching experiment that for the three Level 1 (student-specific) parameters. In
took place at 8 di↵erent schools. Let Xi ⇠ N (µi , i ) particular, let
be the observed e↵ect at each school, where the school
⇡ i = (⇡i,1 , . . . , ⇡i,K ) ⇠ Dirichlet(↵1 , . . . , ↵k ) (2)
specific standard error i is known (based mainly on
the sample size used at that school). There are three µi,k ⇠ N (µ0,k , 0,k ) (3)
ways to approach this problem: (1) No pooling. Es- log(⌧i,k ) ⇠ N (log(⌧0,k ), 0,k ) (4)
timate µi separately with no assumption about about
the similarity of µ across schools. (2) Complete pool- This introduces new Level 2 (across student) param-
ing. Set µi = µ0 for all i and estimate µ0 (3) Par- eters: ↵, µ0 , 0 , ⌧ 0 , and 0 . The likelihood for a
tial pooling. Let µi ⇠ N (µ0 , ⌫) and now jointly es- single student, Li (Yi |⇡ i , µi , ⌧ i ) is given (with a suit-
timate µ0 , µ1 , . . . , µ8 . The no pooling approach pro- able chance of variable) by Equation 1. To get the
duces unbiased estimates for each school, but it has complete data likelihood, multiply across the students
the largest standard errors, especially for the smallest units (or to avoid numeric overflow, sum the log likeli-
schools. The complete pooling approach ignores the hoods). If we let ⌦ be the complete parameter (⇡ i ,µi ,
school level variability, but has much smaller standard ⌧ i for each student, plus ↵, µ0 , 0 , ⌧ 0 , and 0 ), then
errors. In the partial pooling approach, the individual I
school estimates are shrunk towards the grand mean, X
L(Y|⌦) = log Li (Yi |⇡ i , µi , ⌧ i ) . (5)
µ0 , with the amount of shrinkage related to the size of
i=1
the ratio ⌫ 2 /(⌫ 2 + i 2 ); in particular, there is more
shrinkage for the schools which were less precisely mea-
Hierarchical models have their own pathologies which
sured. Note that the higher level standard deviation,
require care during estimation. If either of the stan-
⌫ controls the amount of shrinkage: the smaller ⌫ is
dard deviation parameters, 0,k or 0,k , gets too close
the more the individual school estimates are pulled to-
to zero or infinity, then this could cause the log poste-
wards the grand mean. At the limits, if ⌫ = 0, the par-
rior to go to infinity. These cases correspond to the no
tial pooling model is the same as the complete pooling
pooling and complete pooling extremes of the hierar-
model and if ⌫ = 1 then the partial pooling model is
chical model. Similarly, the varianceP of the Dirichlet
the same as the no pooling model. K
distribution is determined by ↵N = k=1 ↵k . If ↵N
Li (2013) builds a hierarchical mixture model for is too close to zero, this produces no pooling in the
the essay pause data. Figure 1 shows the no pool- estimates of ⇡ i and if ↵N is too large, then there is
3
nearly complete pooling in those estimates. Again,
those values of ↵N can cause the log posterior to be-
come infinite. ↵k = ↵0,k ⇤ ↵N (6)
µi,k = µ0,k + ✓i,k 0,k (7)
✓i,k ⇠ N (0, 1) (8)
log(⌧i,k ) = log(⌧0,k ) + ⌘i,k 0,k (9)
⌘i,k ⇠ N (0, 1) (10)
2.3 Reparameterization
2.4 Prior distributions and parameter
constraints
It is often worthwhile to reparameterize a model in or- Rarely does an analyst have enough prior informa-
der to make MCMC estimation more efficient. Both tion to completely specify a prior distribution. Usu-
random walk Metropolis (Section 3.2) and Hamilto- ally, the analyst chooses a convenient functional form
nian Monte Carlo (Section 3.3) work by moving a point and chooses the hyperparameters of that form based
around the parameter space of the model. The geom- on available prior information (if any). One popular
etry of that space will influence how fast the Markov choice is the conjugate distribution of the likelihood.
chain mixes, that is, moves around the space. If the For example, if the prior for a multinomial probability,
geometry is unfavorable, the point will move slowly ⇡ i follows a Dirichlet distribution with hyperparame-
through the space and the autocorrelation will be high ter, ↵, then the posterior distribution will also be a
(Neal, 2011). In this case a very large Monte Carlo Dirichlet with a hyperparameter equal to the sum of
sample will be required to obtain reasonable Monte the prior hyperparameter and the data. This gives the
Carlo error for the parameter estimates. hyperparameter of a Dirichlet prior a convenient in-
terpretation as pseudo-data. A convenient functional
Consider once more the Eight Schools problem where
form is one based on a normal distribution (sometimes
µi ⇠ N (µ0 , ⌫). Assume that we have a sampler that
on a transformed parameter, such as the log of the
works by updating the value of the µi ’s one at a time
precision) whose mean can be set to a likely value of
and then updating the values of µ0 and ⌫. When up-
the parameter and whose standard deviation can be
dating µi , if ⌫ is small then values of µi close to µ0 will
set so that all of the likely values for the parameter
have the highest conditional probability. When updat-
are within two standard deviations of the mean. Note
ing ⌫ if all of the values of µi are close to µ0 , then small
that for location parameters, the normal distribution
values of ⌫ will have the highest conditional probabil-
is often a conjugate prior distribution.
ity. The net result is a chain in which the movement
from states with small ⌫ and the µi ’s close together Proper prior distributions are also useful for keeping
to states with large ⌫ and µi ’s far apart takes many parameter values away from degenerate or impossible
steps. solutions. For example, priors for standard deviations,
0,k and 0,k must be strictly positive. Also, the group
A simple trick produces a chain which moves much
precisions for each student, ⌧ i , must be strictly pos-
more quickly. Introduce a series of auxiliary variables
itive. The natural conjugate distribution for a preci-
✓i ⇠ N (0, 1) and set µi = µ0 + ⌫✓i . Note that the
sion is a gamma distribution, which is strictly positive.
marginal distribution of µi has not changed, but the
The lognormal distribution, used in Equations 4 and 9,
geometry of the ✓, µ0 , ⌫ space is di↵erent from the ge-
has a similar shape, but its parameters can be inter-
ometry of the µ, µ0 , ⌫ space, resulting in a chain that
preted as a mean and standard deviation on the log
moves much more quickly.
scale. The mixing probabilities ⇡ i must be defined
A similar trick works for modeling the relationships over the unit simplex (i.e., they must all be between 0
between ↵ and ⇡ i . Setting ↵k = ↵0,k ⇤ ↵N , where and 1 and they must sum to 1), as must the expected
↵0 has a Dirichlet distribution and ↵N has a gamma mixing probabilities ↵0 ; the Dirichlet distribution sat-
distribution seems to work well for the parameters ↵ of isfies this constraint and is also a natural conjugate.
the Dirichlet distribution. In this case the parameters Finally, ↵N must be strictly positive; the choice of the
have a particularly easy to interpret meaning: ↵0 is chi-squared distribution as a prior ensures this.
the expected value of the Dirichlet distribution and ↵N
There are other softer restraints on the parameters. If
is the e↵ective sample size of the Dirichlet distribution.
0,k or 0,k gets too high or too low for any value of
Applying these two tricks to the parameters of the k, the result is a no pooling or complete pooling so-
hierarchical model yields the following augmented pa- lution on that mixture component. For both of these
rameterization. parameters (.01, 100) seems like a plausible range. A
4
lognormal distribution which puts the bulk of its prob- comes up with a single parameter estimate but does
ability mass in that range should keep the model away not explore the entire space of the distribution. In con-
from those extremes. Values of ⌧i,k that are too small trast, MCMC algorithms explore the entire posterior
represent the collapse of a mixture component onto distribution by taking samples from the space of possi-
a single point. Again, if ⌧0,k is mostly in the range ble values. Although the MCMC samples are not inde-
(.01, 100) the chance of an extreme value of ⌧i,k should pendent, if the sample is large enough it converges in
be small. This yields the following priors: distribution to the desired posterior. Two approaches
to MCMC estimation are the random walk Metropolis
log( 0k ) ⇠ N (0, 1) (11) algorithm (RWM; used by JAGS, Section 3.2) and the
log( 0k ) ⇠ N (0, 1) (12) Hamiltonian Monte Carlo algorithm (HMC; used by
log(⌧0k ) ⇠ N (0, 1) (13) Stan, Section 3.3).
3.1 EM Algorithm
High values of ↵N also correspond to a degenerate so-
lution. In general, the gamma distribution has about McLachlan and Krishnan (2008) provides a review of
the right shape for the prior distribution of ↵N , and we the EM algorithm with emphasis on mixture models.
expect it to be about the same size as I, the number The form of the EM algorithm is particularly simple
of Level-2 units. The choice of prior is a chi-squared for the special case of a non-hierarchical mixture of
distribution with 2 ⇤ I degrees of freedom. normals. It alternates between an E-step where the
2 p(Zi,j = k) = pi,j,k is calculated for every observation,
↵N ⇠ (I ⇤ 2) (14)
j, and every component, k and an M-step where the
maximum likelihood values for ⇡i,k , µi,k and i,k are
The two remaining parameters we don’t need to con- found by taking moments of the data set Yi weighted
strain too much. For µ0,k we use a di↵use normal by the component probabilities, pi,j,k .
prior (one with a high standard deviation), and for ↵0
we use a Je↵rey’s prior (uniform on the logistic scale) A number of di↵erent problems can crop up with using
which is the Dirichlet distribution with all values set EM to fit mixture models. In particular, if ⇡i,k goes to
to 1/2. zero for any k, that component essentially disappears
from the mixture. Also, if i,k goes to zero the mix-
µ0,k ⇠ N (0, 1000) (15) ture component concentrates on a single point. Fur-
↵0 ⇠ Dirichlet(0.5, . . . , 0.5) (16) thermore, if µi,k = µi,k0 and i,k = i,k0 for any pair
of components the result is a degenerate solution with
K 1 components.
The informative priors above are not sufficient to al-
ways keep the Markov chain out of trouble. In partic- As the posterior distribution for the mixture model is
ular, it can still reach places where the log posterior multimodal, the EM algorithm only finds a local max-
distribution is infinite. There are two di↵erent places imum. Running it from multiple starting points may
where these seem to occur. One is associated with high help find the global maximum; however, in practice it
values of ↵N . Putting a hard limit of ↵N < 500 seems is typically run once. If the order of the components
to avoid this problem (when I = 100). Another possi- is important, the components are typically relabeled
ble degenerate spot is when ↵k ⇡ 0 for some k. This after fitting the model according to a predefined rule
means that ⇡i,k will be essentially zero for all students. (e.g., increasing values of µi,k ).
Adding .01 to all of the ↵k values in Equation 2 seems Two packages are available for fitting non-hierarchical
to fix this problem. mixture models using the EM algorithm in R (R Core
Team, 2014): FlexMix (Gruen & Leisch, 2008) and
mixtools (Benaglia, Chauveau, Hunter, & Young,
⇡ i = (⇡i,1 , . . . , ⇡i,K ) ⇠ Dirichlet(↵1 +.01, . . . , ↵k +.01) 2009). These two packages take di↵erent approaches
to how they deal with degenerate solutions. FlexMix
will combine two mixture components if they get too
3 Estimation Algorithms
close together or the probability of one component gets
too small (by default, if ⇡i,k < .05). Mixtools, on the
There are generally two classes of algorithms used
other hand, retries from a di↵erent starting point when
with both hierarchical and mixture models. The ex-
the the EM algorithm converges on a degenerate so-
pectation maximization (EM) algorithm (Section 3.1)
lution. If it exceeds the allowed number of retries, it
searches for a set of parameter values that maximizes
gives up.
the log posterior distribution of the data (or if the prior
distribution is flat, it maximizes the log likelihood). It Neither mixtools nor FlexMix provides standard er-
5
rors for the parameter estimates. The mixtools pack- ter for the correctness of the algorithm, the most com-
age recommends using the bootstrap (resampling from mon method is to go one parameter at a time and add
the data distribution) to calculate standard errors, and a random o↵set (a step) to its value, accepting or re-
provides a function to facilitate this. jecting it according to the Metropolis rule. As this
distribution is essentially a random walk over the pa-
rameter space, this implementation of MCMC is called
3.2 Random-walk Metropolis Algorithm
random walk Metropolis (RWM). The step size is a
(RWM; used by JAGS)
critical tuning parameter. If the average step size is
Geyer (2011) gives a tutorial summary of MCMC. The too large, the value will be rejected nearly every cycle
basic idea is that a mechanism is found for construct- and the autocorrelation will be high. If the step size
ing a Markov chain whose stationary distribution is the is too small, the chain will move very slowly through
desired posterior distribution. The chain is run until the space and the autocorrelation will be high. Gibbs
the analyst is reasonably sure it has reach the station- sampling, where the step is chosen using a conjugate
ary distribution (these early draws are discarded as distribution so the Metropolis-Hastings ratio always
burn-in). Then the the chain is run some more until it accepts, is not necessarily better. Often the e↵ective
is believed to have mixed throughout the entire poste- step size of the Gibbs sampler is small resulting in high
rior distribution. At this point it has reached pseudo- autocorrelation.
convergence (Geyer calls this pseudo-convergence, be- Most packages that use RWM do some adaptation on
cause without running the chain for an infinite length the step size, trying to get an optimal rejection rate.
of time, there is no way of telling if some part of the During this adaptation phase, the Markov chain does
parameter space was never reached.) At this point not have the correct stationary distribution, so those
the mean and standard error of the parameters are observations must be discarded, and a certain amount
estimated from the the observed mean and standard of burn-in is needed after the adaptation finishes.
deviation of the parameter values in the MCMC sam-
ple. MCMC and the RWM algorithm were made popular
by their convenient implementation in the BUGS soft-
There are two sources of error in estimates made from ware package (Thomas, Spiegelhalter, & Gilks, 1992).
the MCMC sample. The first arises because the ob- With BUGS, the analyst can write down the model
served data are a sample from the universe of potential in a form very similar to the series of equations used
observations. This sampling error would be present to describe the model in Section 2, with a syntax de-
even if the posterior distribution could be computed rived from the R language. BUGS then compiles the
exactly. The second is the Monte Carlo error that model into pseudo-code which produces the Markov
comes from the estimation of the posterior distribu- chain, choosing to do Gibbs sampling or random walk
tion with the Monte Carlo sample. Because the draws Metropolis for each parameter depending on whether
from the Markov chain are not statistically indepen- or not a convenient conjugate proposal distribution
dent,pthis Monte Carlo error does not fall at the rate was available. The output could be exported in a
of 1/ R (where R is the number of Monte Carlo sam- form that could be read by R, and the R package coda
ples). It is also related to the autocorrelation (corre- (Plummer, Best, Cowles, & Vines, 2006) could be used
lation between successive draws) of the Markov chain. to process the output. (Later, WinBUGS would build
The higher the autocorrelation, the lower the e↵ective some of that output processing into BUGS.)
sample size of the Monte Carlo sample, and the higher
the Monte Carlo error. Although BUGS played an important role in encourag-
ing data analysts to use MCMC, it is no longer actively
Most methods for building the Markov chain are based supported. This means that the latest developments
on the Metropolis algorithm (see Geyer, 2011, for de- and improvements in MCMC do not get incorporated
tails). A new value for one or more parameter is pro- into its code. Rather than use BUGS, analysts are
posed and the new value is accepted or rejected ran- advised to use one of the two successor software pack-
domly according to the ratio of the posterior distribu- ages: OpenBUGS (Thomas, O’Hara, Ligges, & Sturtz,
tion and the old and new points, with a correction fac- 2006) or JAGS (just another Gibbs sampler; Plummer,
tor for the mechanism used to generate the new sam- 2012). The R package rjags allows JAGS to be called
ple. This is called a Metropolis or Metropolis-Hastings from R, and hence allows R to be used as a scripting
update (the latter contains a correction for asymmet- language for JAGS, which is important for serious an-
ric proposal distributions). Gibbs sampling is a special alytic e↵orts. (Similar packages exist for BUGS and
case in which the proposal is chosen in such a way that OpenBUGS.)
it will always be accepted.
As the form of the proposal distribution does not mat-
6
3.3 Hamiltonian Monte Carlo (HMC; used 3.4 Parallel Computing and Memory Issues
by Stan)
As most computers have multiple processors, paral-
lel computing can be used to speed up MCMC runs.
Hamiltonian Monte Carlo (HMC) (Neal, 2011) is a
Often multiple Markov chains are run and the results
variant on the Metropolis Algorithm which uses a dif-
are compared to assess pseudo-convergence and then
ferent proposal distribution than RWM. In HMC, the
combined for inference (Gelman & Shirley, 2011). This
current draw from the posterior is imagined to be a
is straightforward using the output processing package
small particle on a hilly surface (the posterior distri-
coda. It is a little bit trickier using the rstan package,
bution). The particle is given a random velocity and is
because many of the graphics require a full stanfit
allowed to move for several discrete steps in that direc-
object. However, the conversion from Stan to coda
tion. The movement follows the laws of physics, so the
format for the MCMC samples is straightforward.
particle gains speed when it falls down hills and loses
speed when it climbs back up the hills. In this manner In the case of hierarchical mixture models, there is
a proposal is generated that can be a great distance an even easier way to take advantage of multiple pro-
from the original starting point. The proposed point is cesses. If the number of components, K, is unknown,
then accepted or rejected according to the Metropolis the usual procedure is to take several runs with dif-
rule. ferent values of K and compare the fit. Therefore, if
4 processors were available, one could run all of the
The software package Stan (Stan Development Team,
chains for K = 2, one for K = 3, and one for K = 4,
2013) provides support for HMC. As with BUGS and
leaving one free to handle operating system tasks.
JAGS, the model of Section 2 is written in pseudo-
code, although this time the syntax looks more like In most modern computers, the bottleneck is usually
C++ than R. Rather than translate the model into not available CPU cycles, but available memory. For
interpreted code, Stan translates it into C++ then running 3 chains with I = 100 students and R = 5000
compiles and links it with existing Stan code to run MCMC samples in each chain, the MCMC sample can
the sampler. This has an initial overhead for the com- take up to 0.5GB of memory! Consequently, it is crit-
pilation, but afterwards, each cycle of the sampler runs ically important to monitor memory usage when run-
faster. Also, as HMC generally has lower autocorrela- ning multiple MCMC runs. If the computer starts re-
tion than random walk Metropolis, smaller run lengths quiring swap (disk-based memory) in addition to phys-
can be used, making Stan considerably faster than ical memory, then running fewer processes will proba-
JAGS in some applications. A package rstan is avail- bly speed up the computations.
able to link Stan to R, but because of the compilation
Another potential problem occurs when storing the re-
step, it requires that the user have the proper R de-
sult of each run between R sessions in the .RData file.
velopment environment set up.
Note that R attempts to load the entire contents of
HMC has more tuning parameters than random walk that data file into memory when starting R. If there
Metropolis: the mass of the particle, the distribution are the results of several MCMC runs in there, the
of velocities and the number of steps to take in each .RData file can grow to several GB in size, and R can
direction must be selected to set up the algorithm. take several minutes to load. (In case this problem
Stan uses a warm-up phase to do this adaptation. The arises, it is recommended that you take a make a copy
recommended procedure is to use approximately 1/2 of the .RData after the data have been cleaned and all
the samples for warm-up as a longer warm-up produces the auxiliary functions loaded but before starting the
lower autocorrelations when actively sampling. MCMC runs, and put it someplace safe. If the .RData
file gets to be too large it can be simply replaced with
Stan has some interesting features that are not present
the backup.)
in BUGS or JAGS. For one, it does not require every
parameter to have a proper prior distribution (as long In order to prevent the .RData file from growing
as the posterior is proper). It will simply put a uni- unmanageably large, it recommended that the work
form prior over the space of possible values for any space not be saved at the end of each run. Instead,
parameter not given a prior distribution. However, run a block of code like this
using explicit priors has some advantages for the ap-
plication to student pause data. In particular, when assign(runName,result1)
data for a new student become available, the poste- outfile <-
rior parameters for the previous run can be input into gzfile(paste(runName,"R","gz",sep="."),
Stan (or JAGS) and the original calibration model can open="wt")
be reused to estimate the parameters for new student dump(runName,file=outfile)
(Mislevy, Almond, Yan, & Steinberg, 1999). close(outfile)
7
after all computations are completed. Here result1 5. Identify the mixture components. Ideally, the
is a variable that gathers together the portion of the Markov chains have visited all of the modes of the
results to keep, and runName is a character string that posterior distribution, including the ones which
provides a unique name for the run. di↵er only by a permutation of the component la-
bels. Section 4.1 describes how to permute the
Assume that all of the commands necessary to perform
component labels are permuted so that the com-
the MCMC analysis are in a file script.R. To run this
ponent labels are the same in each MCMC cycle.
from a command shell use the command:
6. Check pseudo-convergence. Several statistics and
R CMD BATCH --slave script.R plots are calculated to see if the chains have
reached pseudo-convergence and the sample size
This will run the R code in the script, using the .RData
is adequate. If the sample is inadequate, then ad-
file in the current working directory, and put the out-
ditional samples are collected (Section 4.3).
put into script.Rout. The --slave switch performs
two useful functions. First, it does not save the .RData 7. Draw Inferences. Summaries of the posterior dis-
file at the end of the run (avoiding potential memory tribution for the the parameters of interest are
problems). Second, it suppresses echoing the script to computed and reported. Note that JAGS o↵ers
the output file, making the output file easier to read. some possibilities here that Stan does not. In par-
Under Linux and Mac OS X, the command ticular, JAGS can monitor just the cross-student
parameters (↵0 , ↵N , µ0 , 0 , log(⌧ 0 ), and 0 ) for
nohup R CMD BATCH --slave script.R & a much longer time to check pseudo-convergence,
and then a short additional run can be used to
runs R in the background. The user can log out of the
draw inferences about the student specific param-
computer, but as long as the computer remains on, the
eters, ⇡ i , µi and ⌧ i (for a considerable memory
process will continue to run.
savings).
4 Model Estimation 8. Data point labeling. In mixture models, it is some-
times of interest to identify which mixture com-
A MCMC analysis always follows a number of similar ponent each observation Yi,j comes from (Sec-
steps, although the hierarhical mixture model requires tion 4.5).
a couple of additional steps. Usually, it is best to write 9. Calculate model fit index. If the goal is to compare
these as a script because the analysis will often need models for several di↵erent values of K, then a
to be repeated several times. The steps are as follows: measure of model fit such as DIC or WAIC should
be calculated (Section 5).
1. Set up parameters for the run. This includes
which data to analyze, how many mixture compo- 4.1 Component Identification
nents are required (i.e., K), how long to run the
Markov chain, how many chains to run, what the The multiple modes in the posterior distribution for
prior hyperparameters are. mixture models present a special problem for MCMC.
In particular, it is possible for a Markov chain to get
2. Clean the data. In the case of hierarchical mixture stuck in a mode corresponding to a single component
models, student data vectors which are too short labeling and never mix to the other component la-
should be eliminated. Stan does not accept miss- belings. (This especially problematic when coupled
ing values, so NAs in the data need to be replaced with the common practice of starting several paral-
with a finite value. For both JAGS and Stan the lel Markov chains.) Früthwirth-Schnatter (2001) de-
data are bundled with the prior hyperparameters scribes this problem in detail.
to be passed to the MCMC software.
One solution is to constrain the parameter space to fol-
3. Set initial values. Initial values need to be chosen low some canonical ordering (e.g., µi,1 µi,2 · · ·
for each Markov chain (see Section 4.2). µi, K). Stan allows parameter vectors to be specified
as ordered, that is restricted to be in increasing order.
4. Run the Markov Chain. For JAGS, this consists
This seems tailor-made for the identification issue. If
of four substeps: (a) run the chain in adaptation
an order based on µ0 is desired, the declaration:
mode, (b) run the chain in normal mode for the
burn-in, (c) set up monitors on the desired param- ordered[K] mu0;
eters, and (d) run the chain to collect the MCMC
sample. For Stan, the compilation, warm-up and enforces the ordering constraint in the MCMC sam-
sampling are all done in a single step. pler. JAGS contains a sort function which can achieve
8
a similar e↵ect. Früthwirth-Schnatter (2001) rec- for each Markov chain c: do
ommends against this solution because the resulting for each MCMC sample r in Chain c: do
Markov chains often mix slowly. Some experimenta- Find a permutation of indexes k10 , . . . , kK 0
so
tion with the ordered restriction in Stan confirmed this that !c,r,k10 · · · !c,r,k10 .
finding; the MCMC runs took longer and did not reach for ⇠ in the Level 2 parameters
pseudo-convergence as quickly when the ordered re- {↵0 , µ0 , 0 , ⌧ 0 , 0 }: do
striction was not used. Replace ⇠ c,r with (⇠c,r,k10 , . . . , ⇠c,r,kK
0 ).
end for{Level 2 parameter}
Früthwirth-Schnatter (2001) instead recommends let- if inferences about Level 1 parameters are de-
ting the Markov chains mix across all of the possi- sired then
ble modes, and then sorting the values according to for ⇠ in the Level 1 parameters {⇡, µ, ⌧ , }
the desired identification constraints post hoc. JAGS do
provides special support for this approach, o↵ering a for each Level 2 unit (student), i do
special distribution dnormmix for mixtures of normals. Replace ⇠ c,r,i with
This uses a specially chosen proposal distribution to (⇠c,r,i,k10 , . . . , ⇠c,r,i,kK
0 ).
encourage jumping between the multiple modes in the end for{Level 1 unit}
posterior. The current version of Stan does not pro- end for{Level 1 parameter}
vide a similar feature. end if
One result of this procedure is that the MCMC sample end for{Cycle}
will have a variety of orderings of the components; each end for{Markov Chain}
draw could potentially have a di↵erent labeling. For
example, the MCMC sample for the parameter µi,1 Figure 3: Component Identification Algorithm
will in fact be a mix of samples from µi,1 , . . . , µi,K .
The average of this sample is not likely to be a good
reach convergence more quickly. Gelman et al. (2013)
estimate of µi,1 . Similar problems arise when looking
recommend starting at the maximum likelihood esti-
at pseudo-convergence statistics which are related to
mate when that is available. Fortunately, o↵-the shelf
the variances of the parameters.
software is available for finding the maximum likeli-
To fix this problem, the component labels need to be hood estimate of the individual student mixture mod-
permuted separately for each cycle of the MCMC sam- els. These estimates can be combined to find starting
ple. With a one-level mixture model, it is possible values for the cross-student parameters. The initial
to identify the components by sorting on one of the values can be set by the following steps:
student-level parameters, ⇡i,k , µi,k or ⌧i,k . For the
hierarchical mixture model, one of the cross-student 1. Fit a separate mixture model for students using
parameters, ↵0,k , µ0,k , or ⌧0,k , should be used. Note maximum likelihood. The package mixtools is
that in the hierarchical models some of the student- slightly better for this purpose than FlexMix as
level parameters might not obey the constraints. In it will try multiple times to get a solution with
the case of the pause time analysis, this is acceptable the requested number of components. If the EM
as some students may behave di↵erently from most of algorithm does not converge in several tries, set
their peers (the ultimate goal of the mixture modeling the values of the parameters for that student to
is to identify students who may benefit from di↵erent NA and move on.
approaches to instruction). Choosing an identification
rule involves picking which parameter should be used 2. Identify Components. Sort the student-level pa-
for identification at Level 2 (cross-student) and using rameters, ⇡ i , µi , i and ⌧ i according to the de-
the corresponding parameter is used for student-level sired criteria (see Section 4.1.
parameter identification.
3. Calculate the Level 2 initial values. Most of the
Component identification is straightforward. Let ! be cross-student parameters are either means or vari-
the parameter chosen for model identification. Fig- ances of the student-level parameters. The initial
ure 3 describes the algorithm. value for ↵0 is the mean of ⇡ i across the students
i (ignoring NAs). The initial values for µ0 and 0
4.2 Starting Values are the mean and standard deviation of µi . The
initial values for log(⌧ 0 ) and 0 are the mean and
Although the Markov chain should eventually reach standard deviation of log(⌧ i ).
its stationary distribution no matter where it starts,
starting places that are closer to the center of the dis- 4. Impute starting values for maximum likelihood es-
tribution are better in the sense that the chain should timates that did not converge. These are the NAs
9
from Step 1. For each student i that did not the size of the Monte Carlo error could still be an is-
converge in Step 1, set ⇡ i = ↵0 , µi = µ0 and sue. Both rstan and coda compute an e↵ective sample
⌧ i = ⌧ 0. size for the Monte Carlo sample. This is the size of
a simple random sample from the posterior distribu-
5. Compute initial values for ✓ i and ⌘ i . These can tion that would have the same Monte Carlo error as
be computed by solving Equations 7 and 9 for ✓ i the obtained dependent sample. This is di↵erent for
and ⌘ i . di↵erent parameters. If the e↵ective sample size is too
small, then additional samples are needed.
6. Set the initial value of ↵N = I. Only one param-
eter was not given an initial value in the previous Note that there are K(3I + 5) parameters whose con-
steps, that is the e↵ective sample size for ↵. Set vergence needs to be monitored. It is difficult to
this equal to the number of Level 2 units, I. The achieve pseudo-convergence for all of these parame-
initial values of ↵ can now be calculated as well. ters, and exhausting to check them all. A reasonable
compromise seems to be to monitor only the 5K cross-
This produces a set of initial values for a single chain. student parameters, ↵, µ0 , 0 , log(⌧ 0 ) and 0 . JAGS
Common practice for checking for pseudo-convergence makes this process easier by allowing the analyst to
involves running multiple chains and seeing if they ar- pick which parameters to monitor. Using JAGS, only
rive the same stationary distribution. Di↵erent start- the cross-student parameters can be monitored dur-
ing values for the second and subsequent chains can ing the first MCMC run, and then a second shorter
be found by sampling some fraction of the Level 1 sample of student-level parameters can be be obtained
units from each Level 2 unit. When Ji is large, = .5 through an additional run. The rstan package always
seems to work well. When Ji is small, = .8 seems monitors all parameters, including the ✓ i and ⌘i pa-
to work better. An alternative would be to take a rameters which are not of particular interest.
bootstrap sample (a new sample of size Ji drawn with When I and Ji are large, a run can take several hours
replacement). to several days. As the end of a run might occur dur-
ing the middle of the night, it is useful to have a au-
4.3 Automated Convergence Criteria tomated test of convergence. The rule I have been
b are less than a certain
using is to check to see if all R
Gelman and Shirley (2011) describe recommended maximum (by default 1.1) and if all e↵ective sample
practice for assessing pseudo-convergence of a Markov sizes are greater than a certain minimum (by default
chain. The most common technique is to run several 100) for all cross-student parameters. If these criteria
Markov chains and look at the ratio of within-chain are not met, new chains of twice the length of the old
variance to between-chain variance. This ratio is called chain can be run. The traceplots are saved to a file for
b and it comes in both a univariate (one parameter
R, later examination, as are some other statistics such
at a time) and a multivariate version. It is natural to as the mean value of each chain. (These traceplots
look at parameters with the same name together (e.g., and outputs are available at the web site mentioned
µ0 , 0 , ⌧ 0 , and 0 ). Using the multivariate version of above.) It is generally not worthwhile to restart the
Rb with ↵0 and ⇡ i requires some care because calcu- chain for a third time. If pseudo-convergence has not
lating the multivariate R b involves inverting a matrix been achieved after the second longer run, usually a
that does not have full rank when the parameters are better solution is to reparameterize the model.
restricted to a simplex. The work-around is to only
look at the first K 1 values of these parameters. It is easier to extend the MCMC run using JAGS than
using Stan. In JAGS, calling update again samples
The other commonly used test for pseudo-convergence additional points from the Markov chain, starting at
is to look at a trace plot for each parameter: a time the previous values, extending the chains. The current
series plot of the sampled parameter value against the version of Stan (2.2.0)1 saves the compiled C++ code,
MCMC cycle number. Multiple chains can be plot- but not the warm-up parameter values. So the sec-
ted on top of each other using di↵erent colors. If the ond run is a new run, starting over from the warm-up
chain is converged and the sample size is adequate, the phase. In this run, both the warm-up phase and the
traceplot should look like white noise. If this is not sampling phase should be extended, because a longer
the case, a bigger MCMC sample is needed. While
it is not particularly useful for an automated test of 1
The Stan development team have stated that the abil-
pseudo-convergence, the traceplot is very useful for di- ity to save and reuse the warm-up parameters are a high
agnosing what is happening when the chains are not priority for future version of Stan. Version 2.3 of Stan was
converging. Some sample traceplots are given below. released between the time of first writing and publication
of the paper, but the release notes do not indicate that the
Even if the chains have reached pseudo-convergence, restart issue was addressed.
10
warm-up may produce a lower autocorrelation. In ei- timated number of components) was odd or vice versa.
ther case, the data from both runs can be pooled for The values of Rb for the constrained Stan model were
inferences. substantially worse, basically confirming the findings
of Früthwirth-Schnatter (2001).
4.4 A simple test The e↵ective sample sizes were much better for Stan
and HMC. For the K 0 = 2 case, the smallest e↵ective
To test the scripts and the viability of the model, I sample size ranged from 17–142, while for the uncon-
created some artificial data consistent with the model strained Stan model it ranged from 805–3084. Thus,
in Section 2. To do this, I took one of the data sets roughly 3 times the CPU time is producing a 5-fold
use by Li (2013) and ran the initial parameter algo- decrease in the Monte Carlo error.
rithm described in Section 4.2 for K = 2, 3 and 4.
I took the cross-student parameters from this exer- The following graphs compare the JAGS and Stan
cise, and generated random parameters and data sets (Unconstrained model) outputs for four selected cross-
for 10 students. This produced three data sets (for student parameters. The output is from coda and pro-
K = 2, 3 and 4) which are consistent with the model vides a traceplot on the left and a density plot for the
and reasonably similar to the observed data. All corresponding distribution on the right. Recall that
three data sets and the sample parameters are avail- the principle di↵erence between JAGS and Stan is the
able on the web site, http://pluto.coe.fsu.edu/ proposal distribution: JAGS uses a random walk pro-
mcmc-hierMM/. posal with a special distribution for the mixing param-
eters and Stan uses the Hamiltonian proposal. Also,
For each data set, I fit three di↵erent hierarchical mix- note that for Stan the complete run of 15,000 samples
ture models (K 0 = 2, 3 and 4, where K 0 is the value is actually made up of a sort run of length 5,000 and
used for estimation) using three di↵erent variations a longer run of length 10,000; hence there is often a
of the algorithm: JAGS (RWM), Stan (HMC) with discontinuity at iteration 5,000.
the cross-student means constrained to be increasing,
and Stan (HMC) with the means unconstrained. For Although not strictly speaking a parameter, the de-
the JAGS and Stan unconstrained run the chains were viance (twice the negative log likelihood) can easily
sorted (Section 4.1) on the basis of the cross-students be monitored in JAGS. Stan does not o↵er a deviance
means (µ0 ) before the convergence tests were run. In monitor, but instead monitors the log posterior, which
each case, the initial run of 3 chains with 5,000 ob- is similar. Both give an impression of how well the
servations was followed by a second run of 3 chains model fits the data. Figure 4 shows the monitors for
with 10,000 observations when the first one did not the deviance or log posterior. This traceplot is nearly
converge. All of the results are tabulated at the web ideal white noise, indicating good convergence for this
value. The value of R b is less than the heuristic thresh-
site, including links to the complete output files and
all traceplots. old of 1.1 for both chains, and the e↵ective sample size
is about 1/6 of the 45,000 total Monte Carlo observa-
Unsurprisingly, the higher the value of K 0 (number of tions.
components in the estimated model), the longer the
run took. The runs for K 0 = 2 to between 10–15 min- Figure 5 shows the grand mean of the log pause times
utes in JAGS and from 30–90 minutes in Stan, while for the first component. The JAGS output (upper
the runs for K 0 = 4 too just less than an hour in JAGS row) shows a classic slow mixing pattern: the chains
and between 2–3 hours in Stan. The time di↵erence are crossing but moving slowly across the support of
between JAGS and Stan is due to the di↵erence in the parameter space. Thus, for JAGS even though
the chains have nearly reached pseudo-convergence (R b
the algorithms. HMC uses a more complex proposal
distribution than RWM, so naturally each cycle takes is just slightly greater than 1.1), the e↵ective sample
longer. The goal of the HMC is to trade the longer size is only 142 observations. A longer chain might
run times for a more efficient (lower autocorrelation) be needed for a good estimate of this parameter. The
sample, so that the total run time is minimized. The Stan output (bottom row) looks much better, and the
constrained version of Stan seemed to take longer than e↵ective sample size is a respectable 3,680.
the unconstrained. The Markov chains for ↵0,1 (the proportion of pause
In all 27 runs, chain lengths of 15,000 cycles were not times in the first component) have not yet reached
sufficient to reach pseudo-convergence. The maximum pseudo convergence, but they are close (a longer run
b varied
(across all cross-student parameters) value for R might get them there). Note the black chain often
between 1.3 and 2.2, depending on the run. It seemed ranges above the values in the red and green chains
to be slightly worse for cases in which K (number of in the JAGS run (upper row). This is an indication
components for data generation) was even and K 0 (es- that the chains may be stuck in di↵erent modes; the
11
Figure 4: Deviance (JAGS) and Log Posterior (Stan) plots
Note: To reduce the file size of this paper, this is a bitmap picture of the traceplot. The original pdf version is
available at http://pluto.coe.fsu.edu/mcmc-hierMM/DeviancePlots.pdf.
12
Figure 5: Traceplots and Density plots for µ0,1
Note: To reduce the file size of this paper, this is a bitmap picture of the traceplot. The original pdf version is
available at http://pluto.coe.fsu.edu/mcmc-hierMM/mu0Plots.pdf.
13
Figure 6: Traceplots and Density plots for ↵0,1
Note: To reduce the file size of this paper, this is a bitmap picture of the traceplot. The original pdf version is
available at http://pluto.coe.fsu.edu/mcmc-hierMM/alpha0Plots.pdf.
14
shoulder on the density plot also points towards mul- 4.6 Additional Level 2 Units (students)
tiple modes. The plot for Stan looks even worse, after
iteration 5,000 (i.e., after the chain was restarted), the In the pause time application, the goal is to be able
red chain has moved into the lower end of the support to estimate fairly quickly the student-level parameters
of the distribution. This gives a large shoulder in the for a new student working on the same essay. As the
corresponding density plot. MCMC run is time consuming, a fast method for an-
alyzing data from new student would be useful.
The chains for the variance of the log precisions for the
first component, 0,1 (Figure 7), are far from pseudo- One approach would be to simply treat Equations 2,
convergence with R b = 1.73. In the JAGS chains (top 3 and 4 as prior distributions, plugging in the point
row) the black chain seems to prefer much smaller val- estimates for the cross-student parameters as hyper-
ues for this parameter. In the Stan output, the green parameters, and find the posterior mode of the ob-
chain seems to be restricted to the smaller values in servations from the new data. There are two problems
both the first and second run. However, between the with this approach. First, the this method ignores any
first and second run the black and red chains swap remaining uncertainty about the cross-student param-
places. Clearly this is the weakest spot in the model, eters. Second, the existing mixture fitting software
and perhaps a reparameterization here would help. packages (FlexMix and mixtools) do not provide any
way to input the prior information.
Looking at Figures 6 and 7 together, a pattern
emerges. There seems to be (at least) two distinct A better approach is to do an additional MCMC run,
modes. One mode (black chain JAGS, green chain in using the posterior distributions for the cross-student
Stan) has higher value for ↵0,1 and a lower value for parameters from the previous run as priors from the
previous run (Mislevy et al., 1999). If the hyperpa-
0,1 , and the other one has the reverse pattern. This
is an advantage of the MCMC approach over an EM rameters for the cross-student priors are passed to
approach. In particular, if the EM algorithm was only the MCMC sampler as data, then the same JAGS or
run once from a single starting point, the existence of Stan model can be reused for additional student pause
the other mode might never be discovered. records. Note that in addition to the MCMC code
Stan provides a function that will find the maximum
of the posterior. Even if MCMC is used at this step,
4.5 Data Point Labeling
the number of Level 2 units, I, will be smaller and the
For some applications, the analysis can stop after es- chains will reach pseudo-convergence more quickly.
timating the parameters for the students units, ⇡ i , µi If the number of students, I, is large in the original
and ⌧ i (or i as desired). In other applications it is sample (as is true in the original data set), then a
interesting to assign the individual pauses to compo- sample of students can be used in the initial model
nents, that is, to assign values to Zi,j . calibration, and student-level parameters for the oth-
When ⇡ i , µi and i are known, it is simple to calculate ers can be estimated later using this method. In par-
(c,r) (c,r) (c,r) ticular, there seems to be a paradox estimating mod-
pi,j,k = Pr(Zi,j = k). Let ⇡ i , µi and i be
els using large data sets with MCMC. The larger the
the draws of the Level 1 parameters from Cycle r of
data, the slower the run. This is not just because
Chain c (after the data have been sorted to identify
each data point needs to be visited within each cycle
the components). Now define,
of each Markov chain. It appears, especially with hi-
(c,r) erarchical models, that more Level 2 units cause the
(c,r) (c,r) Yi,j µi,k
qi,j,k = ⇡i,k ( ), (17) chain to have higher autocorrelation, meaning larger
(c,r)
i,k runs are needed to achieve acceptable Monte Carlo er-
ror. For the student keystroke logs, the “initial data”
(c,r) PK
(c,r) used for calibration was a sample of 100 essays from
and set Qi,j =
k=1 qi,j,k . The distribution for Zi,j
(c,r) (c,r) (c,r) the complete collection. There is a obvious trade-o↵
for MCMC draw (c, r) is then pi,j,k = qi,j,k /Qi,j .
between working with smaller data sets and being able
The MCMC estimate for pi,j,k is just the aver-
(c,r)
to estimate rare events: I = 100 seemed like a good
age over all the MCMC draws of pi,j,k , that is compromise.
1
PC PR (c,r)
CR c=1 r=1 pi,j,k .
If desired, Zi,j can be estimated as the maximum over 5 Model Selection
k of pi,j,k , but for most purposes simply looking at the
probability distribution pi,j provides an adequate, if A big question in hierarchical mixture modeling is
not better summary of the available information about “What is the optimal number of components, K?”
the component from which Yi,j was drawn. Although there is a way to add a prior distribution
15
Figure 7: Traceplots and Density plots for 0,1
Note: To reduce the file size of this paper, this is a bitmap picture of the traceplot. The original pdf version is
available at http://pluto.coe.fsu.edu/mcmc-hierMM/gamma0Plots.pdf.
16
over K and learn this as part of the MCMC pro- Gelman et al. (2013) advocate using a new measure of
cess, it is quite tricky and the number of parameters model selection called WAIC. WAIC makes two sub-
varies for di↵erent values of K. (Geyer, 2011, describes stitutions in the definition of AIC above. First, it uses
the Metropolis-Hastings-Green algorithm which is re- D(⌦) in place of D(⌦). Second, it uses a new way of
quired for variable dimension parameter spaces). In calculating the dimensionality of the model. There are
practice, analysts often believe that the optimal value two possibilities:
of K is small (less than 5), so the easiest way to de- Ji
I X
X
termine a value of K is to fit separate models for
pWAIC1 = 2 (log E[Qi,j (⌦)] E[log Qi, j(⌦)]) ,
K = 2, 3, 4, . . . and compare the fit of the models.
i=1 j=1
Note that increasing K should always improve the fit (21)
of the model, even if the extra component is just used Ji
I X
X
to trivially fit one additional data point. Consequently, pWAIC2 = 2 Var(Qi,j (⌦)) . (22)
analysts use penalized measures of model fit, such as i=1 j=1
AIC to chose among the models.
The expectation and variance are theoretically defined
Chapter 7 of Gelman et al. (2013) gives a good over the posterior distribution and are approximated
overview of the various information criteria, so the using the MCMC sample. The final expression for
current section will provide only brief definitions. Re- WAIC is:
call Equation 5 that described the log likelihood of
the complete data Y as a function of the complete WAIC = D(⌦) + 2 ⇤ mWAIC . (23)
parameter ⌦. By convention, the deviance is twice
the negative log likelihood: D(⌦) = 2L(Y|⌦). The For all of the model fit statistics discussed here: AIC,
lower the deviance, the better the fit of the data to the WAIC1 and WAIC2 , the model with the lowest value
model, evaluated at that parameter. As students are is the best. One way to chose the best value of K is to
independent given the cross-student parameters, and calculate all of these statistics for multiple values of K
pause times are independent given the student-level and choose the value of K which has the lowest value
e for some value of the parameters
parameters, D(⌦), for all 5 statistics. Hopefully, all of the statistics will
e
⌦ can be written as a sum: point at the same model. If there is not true, that is
often evidence that the fit of two di↵erent models is
Ji
I X
X
e = e , too close to call.
D(⌦) log Qi,j (⌦) (18)
i=1 j=1 In the 27 runs using the data from known values of K,
simply using the lowest WAIC value was not sufficient
e is the likelihood of data point Yi,j , that
where Qi,j (⌦) to recover the model. As the Markov chains have not
is: yet reached pseudo-convergence, the WAIC values may
Xk not be correct, but there was fairly close agreement on
e = Yi,j µ gi,k
Qi,j (⌦) g
⇡ i,k ( ). (19) both WAIC1 and WAIC2 for both the JAGS and Stan
g
i,k
k=1 (unconstrained) models. However, the values of both
Note that the cross-student parameters drop out of WAIC1 and WAIC2 were also fairly close for all values
these calculations. of K 0 (number of estimated components). For example
when K = 2 the WAIC1 value was 1132, 1132 and
The Akaike information criterion (AIC) takes the de- 1135 for K 0 = 2, 3 and 4 respectively when estimated
viance as a measure of model fit and penalizes it for with JAGS and 1132, 1131, and 1132 when estimated
b be the
the number of parameters in the model. Let ⌦ with Stan. The results for WAIC2 were similar. It
maximum likelihood estimate of the parameter. Then appears as if the common practice of just picking the
AIC is defined as: best value of WAIC to determine K is not sufficient.
b +2⇤m In particular, it may be possible to approximate the
AIC = D(⌦) (20)
data quite well with a model which has an incorrect
where m = K(2IK + 4) + (K 1)(I + 1) + 1 is the value of K.
number of free parameters in the model. Note that
the MCMC run does not provide a maximum likeli- 6 Conclusions
hood estimate. A pseudo-AIC can be computed by
substituting ⌦, the posterior mean of the parameters The procedure described above is realized as code in
⌦. Note that ⌦ must be calculated after the parame- R, Stan and JAGS and available for free at http://
ters in each MCMC cycle are permuted to identify the pluto.coe.fsu.edu/mcmc-hierMM/. Also available,
components. are the simulated data for trying out the models.
17
These scripts contain the functions necessary to auto- for large data sets indicates that there may be addi-
matically set starting parameters for these problems, tional transformation of this model that are necessary
identify the components, and to calculate WAIC model to make it work properly. Hopefully, the publication of
fit statistics. the source code will allow other to explore this model
more fully.
The scripts seem to work fairly well for a small number
of students units (5 I 10). With this sample size,
although JAGS has the faster run time, Stan has a References
larger e↵ective sample size. Hamiltonian Monte Carlo.
With large sample sizes I = 100, Ji ⇡ 100 even the Almond, R. G., Deane, P., Quinlan, T., Wagner,
HMC has high autocorrelation and both JAGS and M., & Sydorenko, T. (2012). A prelimi-
Stan are extremely slow. Here the ability of JAGS to nary analysis of keystroke log data from a
allow the user to selectively monitor parameters and timed writing task (Research Report No.
to extend the chains without starting over allows for RR-12-23). Educational Testing Service. Re-
better memory management. trieved from http://www.ets.org/research/
policy research reports/publications/
The code supplied at the web site solves two key tech-
report/2012/jgdg
nical problems: the post-run sorting of the mixture
Benaglia, T., Chauveau, D., Hunter, D. R., & Young,
components (Früthwirth-Schnatter, 2001) and the cal-
D. (2009). mixtools: An R package for ana-
culation of the WAIC model-fit indexes. Hopefully,
lyzing finite mixture models. Journal of Sta-
this code will be useful for somebody working on a
tistical Software, 32 (6), 1–29. Retrieved from
similar application.
http://www.jstatsoft.org/v32/i06/
The model proposed in this paper has two problems, Brooks, S., Gelman, A., Jones, G. L., & Meng, X.
even with data simulated according to the model. The (Eds.). (2011). Handbook of markov chain monte
first is the slow mixing. This appears to be a prob- carlo. CRC Press.
lem with multiple modes, and indicates some possible Deane, P. (2012). Rethinking K-12 writing assessment.
non-identifiability in the model specifications. Thus, In N. Elliot & L. Perelman (Eds.), Writing as-
further work is needed on the form of the model. The sessment in the 21st century. essays in honor of
second problem is that the WAIC test does not ap- Edward M. white (pp. 87–100). Hampton Press.
pear to be sensitive enough to recover the number of Früthwirth-Schnatter, S. (2001). Markov chain
components in the model. As the original purpose of Monte Carlo estimation of classical and dynamic
moving to the hierarchical mixture model was to dis- switching and mixture models. Journal of the
cover if there were rare components that could not be American Statistical Association, 96 (453), 194–
detected in the single-level hierarchical model, the cur- 209.
rent model does not appear to be a good approach for Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B.,
that purpose. Vehtari, A., & Rubin, D. B. (2013). Bayesian
In particular, returning to the original application of data analysis (3rd ed.). Chapman and Hall. (The
fitting mixture models to student pauses while writing, third edition is revised and expanded and has
the hierarchical part of the model seems to be creat- material that the earlier editions lack.)
ing as many problems as it solves. It is not clearly Gelman, A., & Shirley, K. (2011). Inference from
better than fitting the non-hierarchical mixture model simulations and monitoring convergence. In
individually to each student. Another problem it has S. Brooks, A. Gelman, G. L. Jones, & X. Meng
for this application is the assumption that each pause (Eds.), Handbook of markov chain monte carlo
is independent. This does not correspond with what (pp. 163–174). CRC Press.
is known about the writing process: typically writ- Geyer, C. J. (2011). Introduction to Markov
ers spend long bursts of activities simply doing tran- chain Monte Carlo. In S. Brooks, A. Gelman,
scription (i.e., typing) and only rarely pause for higher G. L. Jones, & X. Meng (Eds.), Handbook of
level processing (e.g., spelling, grammar, word choice, markov chain monte carlo (pp. 3–48). CRC
planning, etc). In particular, rather than getting this Press.
model to work for this application, it may be better to Gruen, B., & Leisch, F. (2008). FlexMix version
look at hidden Markov models for individual student 2: Finite mixtures with concomitant variables
writings.2 and varying and constant parameters. Journal
of Statistical Software, 28 (4), 1–35. Retrieved
The slow mixing of the hierarchical mixture model from http://www.jstatsoft.org/v28/i04/
Li, T. (2013). A Bayesian hierarchical mixture ap-
2
Gary Feng and Paul Deane, private communication. proach to model timing data with application to
18
writing assessment (Unpublished master’s the-
sis). Florida State University.
McLachlan, G., & Krishnan, T. (2008). The EM al-
gorithm and extensions (2nd ed.). Wiley.
McLachlan, G., & Peel, D. (2000). Finite mixture
models. Wiley.
Mislevy, R. J., Almond, R. G., Yan, D., & Stein-
berg, L. S. (1999). Bayes nets in educational
assessment: Where the numbers come from. In
K. B. Laskey & H. Prade (Eds.), Uncertainty in
artificial intelligence ’99 (p. 437-446). Morgan-
Kaufmann.
Neal, R. M. (2011). MCMC using Hamiltonian dynam-
ics. In MCMCHandbook (Ed.), (pp. 113–162).
Plummer, M. (2012, May). JAGS version 3.2.0
user manual (3.2.0 ed.) [Computer software
manual]. Retrieved from http://mcmc-jags
.sourceforge.net/
Plummer, M., Best, N., Cowles, K., & Vines, K.
(2006). Coda: Convergence diagnosis and out-
put analysis for mcmc. R News, 6 (1), 7–11.
Retrieved from http://CRAN.R-project.org/
doc/Rnews/
R Core Team. (2014). R: A language and environ-
ment for statistical computing [Computer soft-
ware manual]. Vienna, Austria. Retrieved from
http://www.R-project.org/
Stan Development Team. (2013). Stan: A C++ library
for probability and sampling (Version 2.2.0 ed.)
[Computer software manual]. Retrieved from
http://mc-stan.org/
Thomas, A., O’Hara, B., Ligges, U., & Sturtz, S.
(2006). Making BUGS open. R News, 6 , 12-17.
Retrieved from http://www.rni.helsinki.fi/
~boh/publications/Rnews2006-1.pdf
Thomas, A., Spiegelhalter, D. J., & Gilks, W. R.
(1992). BUGS: A program to perform
Bayesian inference using Gibbs sampling. In
J. M. Bernardo, J. O. Berger, A. P. Dawid, &
A. F. M. Smith (Eds.), Bayesian statistics 4. (pp.
837–842). Clarendon Press.
Acknowledgments
Tingxuan Li wrote the first version of much of the
R code as part of her Master’s thesis. Paul Deane
and Gary Feng of Educational Testing Service have
been collaborators on the essay keystroke log analysis
work which is the inspiration for the present research.
Various members of the JAGS and Stan users groups
and development teams have been helpful in answer-
ing emails and postings related to these models, and
Kathy Laskey and several anonymous reviewers have
given me useful feedback (to which I should have paid
closer attention) on the draft of this paper.
19