depmixS4 : An R Package for Hidden Markov Models

This introduction to the R package depmixS4 is a (slightly) modiﬁed version of Visser and Speekenbrink (2010), published in the Journal of Statistical Software . Please refer to that article when using depmixS4 . The current version is 1.3-0; the version history and changes can be found in the NEWS ﬁle of the package. Below, the major versions are listed along with the most noteworthy changes. depmixS4 implements a general framework for deﬁning and estimating dependent mixture models in the R programming language. This includes standard Markov models, latent/hidden Markov models, and latent class and ﬁnite mixture distribution models. The models can be ﬁtted on mixed multivariate data with distributions from the glm family, the (logistic) multinomial, or the multivariate normal distribution. Other distributions can be added easily, and an example is provided with the exgaus distribution. Parameters are estimated by the expectation-maximization (EM) algorithm or, when (linear) constraints are imposed on the parameters, by direct numerical optimization with the Rsolnp or Rdonlp2 routines.


Introduction
Markov and latent Markov models are frequently used in the social sciences, in different areas and applications.In psychology, they are used for modelling learning processes; see Wickens depmixS4: An R Package for Hidden Markov Models (1982), for an overview, and e.g., Schmittmann, Visser, and Raijmakers (2006), for a recent application.In economics, latent Markov models are so-called regime switching models (see e.g., Kim 1994 andGhysels 1994).Further applications include speech recognition (Rabiner 1989), EEG analysis (Rainer and Miller 2000), and genetics (Krogh 1998).In these latter areas of application, latent Markov models are usually referred to as hidden Markov models.See for example Frühwirth-Schnatter (2006) for an overview of hidden Markov models with extensions.Further examples of applications can be found in e.g., Cappe, Moulines, and Ryden (2005, Chapter 1).A more gentle introduction into hidden Markov models with applications is the book by Zucchini and MacDonald (2009).
The depmixS4 package was motivated by the fact that while Markov models are used commonly in the social sciences, no comprehensive package was available for fitting such models.Existing software for estimating Markovian models include Panmark (van de Pol, Langeheine, and Jong 1996), and for latent class models Latent Gold (Vermunt and Magidson 2003).These programs lack a number of important features, besides not being freely available.There are currently some packages in R that handle hidden Markov models but they lack a number of features that we needed in our research.In particular, depmixS4 was designed to meet the following goals: 1. to be able to estimate parameters subject to general linear (in)equality constraints; 2. to be able to fit transition models with covariates, i.e., to have time-dependent transition matrices; 3. to be able to include covariates in the prior or initial state probabilities; 4. to be easily extensible, in particular, to allow users to easily add new uni-or multivariate response distributions and new transition models, e.g., continuous time observation models.
Although depmixS4 was designed to deal with longitudinal or time series data, for say T > 100, it can also handle the limit case when T = 1.In this case, there are no time dependencies between observed data and the model reduces to a finite mixture or latent class model.While there are specialized packages to deal with mixture data, as far as we know these do not allow the inclusion of covariates on the prior probabilities of class membership.The possibility to estimate the effects of covariates on prior and transition probabilities is a distinguishing feature of depmixS4.In Section 2, we provide an outline of the model and likelihood equations.
The depmixS4 package is implemented using the R system for statistical computing (R Development Core Team 2010) and is available from the Comprehensive R Archive Network at http://CRAN.R-project.org/package=depmixS4.

The dependent mixture model
The data considered here have the general form O 1:T = (O  relative reward for speeded and/or accurate responding.These variables are measured on 168, 134 and 137 occasions respectively (the first series of 168 trials is plotted in Figure 1).These data are more fully described in Dutilh, Wagenmakers, Visser, and van der Maas (2011), and in the next section a number of example models for these data is described.
The latent Markov model is usually associated with data of this type, in particular for multinomially distributed responses.However, commonly employed estimation procedures (e.g., van de Pol et al. 1996), are not suitable for long time series due to underflow problems.In contrast, the hidden Markov model is typically only used for 'long' univariate time series (Cappe et al. 2005, Chapter 1).We use the term "dependent mixture model" because one of the authors (Ingmar Visser) thought it was time for a new name to relate these models1 .
The fundamental assumption of a dependent mixture model is that at any time point, the observations are distributed as a mixture with n components (or states), and that timedependencies between the observations are due to time-dependencies between the mixture components (i.e., transition probabilities between the components).These latter dependencies are assumed to follow a first-order Markov process.In the models we are considering here, the mixture distributions, the initial mixture probabilities and transition probabilities can all depend on covariates z t .
In a dependent mixture model, the joint likelihood of observations O 1:T and latent states S 1:T = (S 1 , . . ., S T ), given model parameters θ and covariates z 1:T = (z 1 , . . ., z T ), can be depmixS4: An R Package for Hidden Markov Models written as: where we have the following elements: 1. S t is an element of S = {1 . . .n}, a set of n latent classes or states.
3. a ij (z t ) = P(S t+1 = j|S t = i, z t ), provides the probability of a transition from state i to state j with covariate z t .
4. b St is a vector of observation densities b k j (z t ) = P(O k t |S t = j, z t ) that provide the conditional densities of observations O k t associated with latent class/state j and covariate z t , j = 1, . . ., n, k = 1, . . ., m.
For the example data above, b k j could be a Gaussian distribution function for the response time variable, and a Bernoulli distribution for the accuracy variable.In the models considered here, both the transition probability functions a ij and the initial state probability functions π may depend on covariates as well as the response distributions b k j .

Likelihood
To obtain maximum likelihood estimates of the model parameters, we need the marginal likelihood of the observations.For hidden Markov models, this marginal (log-)likelihood is usually computed by the so-called forward-backward algorithm (Baum and Petrie 1966;Rabiner 1989), or rather by the forward part of this algorithm.Lystig and Hughes (2002) changed the forward algorithm in such a way as to allow computing the gradients of the log-likelihood at the same time.They start by rewriting the likelihood as follows (for ease of exposition the dependence on the model parameters and covariates is dropped here): where P(O 1 |O 0 ) := P(O 1 ).Note that for a simple, i.e., observed, Markov chain these probabilities reduce to P(O t |O 1 , . . ., O t−1 ) = P(O t |O t−1 ).The log-likelihood can now be expressed as: To compute the log-likelihood, Lystig and Hughes (2002) define the following (forward) recursion: where Φ t = n i=1 φ t (i).Combining Φ t = P(O t |O 1:(t−1) ), and equation (3) gives the following expression for the log-likelihood: (6)

Parameter estimation
Parameters are estimated in depmixS4 using the expectation-maximization (EM) algorithm or through the use of a general Newton-Raphson optimizer.In the EM algorithm, parameters are estimated by iteratively maximizing the expected joint log-likelihood of the parameters given the observations and states.Let θ = (θ 1 , θ 2 , θ 3 ) be the general parameter vector consisting of three subvectors with parameters for the prior model, transition model, and response models respectively.The joint log-likelihood can be written as: This likelihood depends on the unobserved states S 1:T .In the Expectation step, we replace these with their expected values given a set of (initial) parameters θ = (θ 1 , θ 2 , θ 3 ) and observations O 1:T .The expected log-likelihood: can be written as: where the expected values ξ t (j, k) = P (S t = k, S t−1 = j|O 1:T , z 1:T , θ ) and γ t (j) = P (S t = j|O 1:T , z 1:T , θ ) can be computed effectively by the forward-backward algorithm (see e.g., Rabiner 1989).The Maximization step consists of the maximization of (9) for θ.As the right hand side of (9) consists of three separate parts, we can maximize separately for θ 1 , θ 2 and θ 3 .In common models, maximization for θ 1 and θ 2 is performed by the nnet.defaultroutine in the nnet package (Venables and Ripley 2002), and maximization for θ 3 by the standard glm routine.Note that for the latter maximization, the expected values γ t (j) are used as prior weights of the observations O k t .The EM algorithm however has some drawbacks.First, it can be slow to converge towards the end of optimization.Second, applying constraints to parameters can be problematic; in particular, EM can lead to wrong parameter estimates when applying constraints.Hence, in depmixS4, EM is used by default in unconstrained models, but otherwise, direct optimization is used.Two options are available for direct optimization using package Rsolnp (Ghalanos and Theußl 2010;Ye 1987), or Rdonlp2 (Tamura 2009;Spellucci 2002).Both packages can handle general linear (in)equality constraints (and optionally also non-linear constraints).

Missing data
Missing data can be dealt with straightforwardly in computing the likelihood using the forward recursion in Equations (4-5).Assume we have observed O 1:(t−1) but that observation O t is missing.The key idea that, in this case, the filtering distribution, the probabilities φ t , should be identical to the state prediction distribution, as there is no additional information to estimate the current state.Thus, the forward variables φ t are now computed as: For later observations, we can then use this latter equation again, realizing that the filtering distribution is technically e.g.P(S t+1 |O 1:(t−1),t+1 ).Computationally, the easiest way to implement this is to simply set

Using depmixS4
Two steps are involved in using depmixS4 which are illustrated below with examples: 1. model specification with function depmix (or with mix for latent class and finite mixture models, see example below on adding covariates to prior probabilities); 2. model fitting with function fit.
We have separated the stages of model specification and model fitting because fitting large models can be fairly time-consuming and it is hence useful to be able to check the model specification before actually fitting the model.

Example data: speed
Throughout this article a data set called speed is used.As already indicated in the introduction, it consists of three time series with three variables: response time rt, accuracy corr, and a covariate, Pacc, which defines the relative pay-off for speeded versus accurate responding.Before describing some of the models that are fitted to these data, we provide a brief sketch of the reasons for gathering these data in the first place.
Response times are a very common dependent variable in psychological experiments and hence form the basis for inference about many psychological processes.A potential threat to such inference based on response times is formed by the speed-accuracy trade-off: different participants in an experiment may respond differently to typical instructions to 'respond as fast and accurate as possible'.A popular model which takes the speed-accuracy trade-off into account is the diffusion model (Ratcliff 1978), which has proven to provide accurate descriptions of response times in many different settings.
One potential problem with the diffusion model is that it predicts a continuous trade-off between speed and accuracy of responding, i.e., when participants are pressed to respond faster and faster, the diffusion model predicts that this would lead to a gradual decrease in accuracy.The speed data set that we analyze below was gathered to test this hypothesis versus the alternative hypothesis stating that there is a sudden transition from slow and accurate responding to fast responding at chance level.At each trial of the experiment, the participant is shown the current setting of the relative reward for speed versus accuracy.The bottom panel of Figure 1 shows the values of this variable.The experiment was designed to investigate what would happen when this reward variable changes from reward for accuracy only to reward for speed only.The speed data that we analyse here are from participant A in Experiment 1 in Dutilh et al. (2011), who provide a complete description of the experiment and the relevant theoretical background.
The central question regarding this data is whether it is indeed best described by two modes of responding rather than a single mode of responding with a continuous trade-off between speed and accuracy.The hallmark of a discontinuity between slow versus speeded responding is that switching between the two modes is asymmetric (see e.g.Van

A simple model
A dependent mixture model is defined by the number of states and the initial state, state transition, and response distribution functions.A dependent mixture model can be created with the depmix function as follows: R> library("depmixS4") R> data("speed") R> set.seed(1)R> mod <-depmix(response = rt ~1, data = speed, nstates = 2, + trstart = runif( 4)) The first line of code loads the depmixS4 package and data(speed) loads the speed data set.The line set.seed(1) is necessary to get starting values that will result in the right model, see more on starting values below.
The call to depmix specifies the model with a number of arguments.The response argument is used to specify the response variable, possibly with covariates, in the familiar format using formulae such as in lm or glm models.The second argument, data = speed, provides the data.frame in which the variables from response are interpreted.Third, the number of states of the model is given by nstates = 2.
Starting values.Note also that start values for the transition parameters are provided in this call by the trstart argument.Starting values for parameters can be provided using depmixS4: An R Package for Hidden Markov Models three arguments: instart for the parameters relating to the prior probabilities, trstart, relating the transition models, and respstart for the parameters of the response models.Note that the starting values for the initial and transition models as well as multinomial logit response models are interpreted as probabilities, and internally converted to multinomial logit parameters (if a logit link function is used).Start values can also be generated randomly within the EM algorithm by generating random uniform values for the values of γ t in (9) at iteration 0. Automatic generation of starting values is not available for constrained models (see below).In the call to fit below, the argument emc=em.control(rand=FALSE)controls the EM algorithm and specifies that random start values should not be generated2 .
Fitting the model and printing results.These statistics can also be extracted using logLik, AIC and BIC, respectively.By comparison, a 1-state model for these data, i.e., assuming there is no mixture, has a log-likelihood of −305.33, and 614.66, and 622.83 for the AIC and BIC respectively.Hence, the 2-state model fits the data much better than the 1-state model.Note that the 1-state model can be specified using mod <-depmix(rt ~1, data = speed, nstates = 1), although this model is trivial as it will simply return the mean and standard deviation of the rt variable.
The summary method of fitted models provides the parameter estimates, first for the prior probabilities model, second for the transition models, and third for the response models.

R> summary(fm)
Initial Since no further arguments were specified, the initial state, state transition and response distributions were set to their defaults (multinomial distributions for the first two, and a Gaussian distribution for the response).The resulting model indicates two well-separated states, one with slow and the second with fast responses.The transition probabilities indicate rather stable states, i.e., the probability of remaining in either of the states is around 0.9.The initial state probability estimates indicate that state 1 is the starting state for the process, with a negligible probability of starting in state 2.

Covariates on transition parameters
By default, the transition probabilities and the initial state probabilities are parameterized using a multinomial model with an identity link function.Using a multinomial logistic model allows one to include covariates on the initial state and transition probabilities.In this case, each row of the transition matrix is parameterized by a baseline category logistic multinomial, meaning that the parameter for the base category is fixed at zero (see Agresti 2002, p. 267 ff., for multinomial logistic models and various parameterizations).The default baseline category is the first state.Hence, for example, for a 3-state model, the initial state probability model would have three parameters of which the first is fixed at zero and the other two are freely estimated.Chung, Walls, and Park (2007) discuss a related latent transition model for repeated measurement data (T = 2) using logistic regression on the transition parameters; they rely on Bayesian methods of estimation.Covariates on the transition probabilities can be specified using a one-sided formula as in the following example: The summary provides all parameters of the model, also the (redundant) zeroes for the baseline category in the multinomial model.The summary also prints the transition probabilities at the zero value of the covariate.Note that scaling of the covariate is useful in this regard as it makes interpretation of these intercept probabilities easier.

Multivariate data
Multivariate data can be modelled by providing a list of formulae as well as a list of family objects for the distributions of the various responses.In above examples we have only used the response times which were modelled as a Gaussian distribution.The accuracy variable in the speed data can be modelled with a multinomial by specifying the following: As can be seen, state 1 has fast response times and accuracy is approximately at chance level (.4747), whereas state 2 corresponds with slower responding at higher accuracy levels (.9021).
Note that by specifying multivariate observations in terms of a list, the variables are considered conditionally independent (given the states).Conditionally dependent variables must be handled as a single element in the list.Effectively, this means specifying a multivariate response model.Currently, depmixS4 has one multivariate response model which is for multivariate normal variables.

Fixing and constraining parameters
Using package Rsolnp (Ghalanos and Theußl 2010) or Rdonlp2 (Tamura 2009), parameters may be fitted subject to general linear (in-)equality constraints.Constraining and fixing parameters is done using the conpat argument to the fit function, which specifies for each parameter in the model whether it's fixed (0) or free (1 or higher).Equality constraints can be imposed by giving two parameters the same number in the conpat vector.When only fixed values are required, the fixed argument can be used instead of conpat, with zeroes for fixed parameters and other values (e.g., ones) for non-fixed parameters.Fitting the models subject to these constraints is handled by the optimization routine solnp or, optionally, by donlp2.
To be able to construct the conpat and/or fixed vectors one needs the correct ordering of parameters which is briefly discussed next before proceeding with an example.
Parameter numbering.When using the conpat and fixed arguments, complete parameter vectors should be supplied, i.e., these vectors should have length equal to the number of parameters of the model, which can be obtained by calling npar(object).Note that this is not the same as the degrees of freedom used e.g., in the logLik function because npar also counts the baseline category zeroes from the multinomial logistic models.Parameters are numbered in the following order: 1. the prior model parameters; 2. the parameters for the transition models; 3. the response model parameters per state (and subsequently per response in the case of multivariate time series).

depmixS4: An R Package for Hidden Markov Models
To see the ordering of parameters use the following: R> setpars(mod, value = 1:npar(mod)) To see which parameters are fixed (by default only baseline parameters in the multinomial logistic models for the transition models and the initial state probabilities model): R> setpars(mod, getpars(mod, which = "fixed")) When fitting constraints it is useful to have good starting values for the parameters and hence we first fit the following model without constraints: R> trst <-c(0.9,0.1, 0, 0, 0.1, 0.9, 0, 0) R> mod <-depmix(list(rt ~1,corr ~1), data = speed, transition = ~Pacc, + nstates = 2, family = list(gaussian(), multinomial("identity")), + trstart = trst, instart = c(0.99,0.01)) R> fm1 <-fit(mod,verbose = FALSE, emc=em.control(rand=FALSE))converged at iteration 23 with logLik: -255.5 After this, we use the fitted values from this model to constrain the regression coefficients on the transition matrix (parameters number 6 and 10): Using summary on the fitted model shows that the regression coefficients are now estimated at the same value of 12.66.Setting elements 13 and 14 of conpat to zero resulted in fixing those parameters at their starting values of 0.5.This means that state 1 can now be interpreted as a 'guessing' state which is associated with comparatively fast responses.Similarly for elements 1 and 2, resulting in fixed initial probabilities.The function llratio computes the likelihood ratio χ 2 -statistic and the associated p-value with appropriate degrees of freedom for testing the tenability of constraints (Dannemann and Holzmann 2007).Note that these arguments (i.e., conpat and conrows) provide the possibility for arbitrary constraints, also between, e.g., a multinomial regression coefficient for the transition matrix and the mean of a Gaussian response model.Whether such constraints make sense is hence the responsibility of the user.

Adding covariates on the prior probabilities
To illustrate the use of covariates on the prior probabilities we have included another data set with depmixS4.The balance data consists of 4 binary items (correct-incorrect) on a balance scale task (Siegler 1981).The data form a subset of the data published in Jansen and van der Maas (2002).Before specifying specifying a model for these data, we briefly describe them.
The balance scale task is a famous task for testing cognitive strategies developed by Jean Piaget (see Siegler 1981).Figure 2 provides an example of a balance scale item.Participants' task is to say to which side the balance will tip when released, or alternatively, whether it will stay in balance.The item shown in Figure 2 is a so-called distance item: the number of weights placed on each side is equal, and only the distance of the weights to the fulcrum differs between each side.
Children in the lower grades of primary school are known to ignore the distance dimension, and base their answer only on the number of weights on each side.Hence, they would typically provide the wrong answer to these distance items.Slightly older children do take distance into account when responding to balance scale items, but they only do so when the number of weights is equal on each side.These two strategies that children employ are known as Rule I and Rule II.Other strategies can be teased apart by administering different items.The balance data set that we analyse here consists of 4 distance items on a balance scale task administered to 779 participants ranging from 5 to 19 years of age.The full set of items consisted of 25 items; other items in the test are used to detect other strategies that children and young adults employ in solving balance scale items (see Jansen and van der Maas 2002, for details).

R> fm
Convergence info: Log likelihood converged to within tol.(relative change) 'log Lik.' -917.5 (df=16) AIC: 1867BIC: 1942 Note here that we define a mix model instead of a depmix model as these data form independent observations.More formally, depmix models extend the class of 'mix' models by adding transition models.As for fitting mix models: as can be seen in Equation 9, the EM algorithm can be applied by simply dropping the second summand containing the transition parameters, and this is implemented as such in the EM algorithms in depmixS4.
As can be seen from the print of the fitted model above, the BIC for this model equals 1941.6.The similarly defined 2-class model for these data has a BIC of 1969.2, and the 4-class model has BIC equal to 1950.4.Hence, the 3-class seems to be adequate for describing these data.The intercept values of the multinomial logistic parameters provide the prior probabilities of each class when the covariate has value zero (note that in this case the value zero does not have much significance, scaling and/or centering of covariates may be useful in such cases).The summary function prints these values.As can be seen from those values, at age zero, the prior probability is overwhelmingly at the second class.Inspection of the response parameters reveals that class 2 is associated with incorrect responding, whereas class 1 is associated with correct responding; class 3 is an intermediate class with guessing behavior.Figure 3 depicts the prior class probabilities as function of age based on the fitted parameters.
As can be seen from Figure 3, at younger ages children predominantly apply Rule I, the wrong strategy for these items.According to the model, approximately 90 % of children at age 5 apply Rule I.The remaining 10 % are evenly split among the 2 other classes.At age 19, almost all participants belong to class 1.Interestingly, prior probability of the 'guess' class 2, first increases with age, and then decreases again.This suggests that children pass through a phase in which they are uncertain or possibly switch between applying different strategies.

Extending depmixS4
The depmixS4 package was designed with the aim of making it relatively easy to add new response distributions (as well as possibly new prior and transition models).To make this possible, the EM routine simply calls the fit methods of the separate response models without needing access to the internal workings of these routines.Referring to equation 9, the EM algorithm calls separate fit functions for each part of the model, the prior probability model, the transition models, and the response models.As a consequence, adding user-specified response models is straightforward.The currently implemented distributions are listed in Table 1.
User-defined distributions should extend the 'response' class and have the following slots: 1. y: The response variable.

x:
The design matrix, possibly only an intercept.

parameters:
A named list with the coefficients and possibly other parameters (e.g., the standard deviation in the normal response model).
4. fixed: A vector of logicals indicating whether parameters are fixed.
The fit method is defined as follows: R> setMethod("fit", "exgaus", + function(object, w) { + if(missing(w)) w <-NULL + y <-object@y + fit <-gamlss(y ~1, weights = w, family = exGAUS(), + control = gamlss.control(n.cyc= 100, trace = FALSE), + mu.start = object@parameters$mu, + sigma.start= exp(object@parameters$sigma), + nu.start = exp(object@parameters$nu)) + pars <-c(fit$mu.coefficients,fit$sigma.coefficients,+ fit$nu.coefficients)+ object <-setpars(object,pars) The fit method defines a gamlss model with only an intercept to be estimated and then sets the fitted parameters back into their respective slots in the 'exgaus' object.See the help for gamlss.distrfor interpretation of these parameters.After defining all the necessary methods for the new response model, we can now define the dependent mixture model using this response model.The function makeDepmix is included in depmixS4 to have full control over model specification, and we need it here.
We first create all the response models that we need as a double list: Next, we define the transition and prior probability models using the transInit function (which produces a transInit model, which also extends the 'response' class): R> trstart <-c(0.9,0.1, 0.1, 0.9) R> transition <-list() depmixS4: An R Package for Hidden Markov Models

Conclusions and future work
depmixS4 provides a flexible framework for fitting dependent mixture models for a large variety of response distributions.It can also fit latent class regression and finite mixture models, although it should be noted that more specialized packages are available for this such as flexmix (Leisch 2004;Grün and Leisch 2008).The package is intended for modelling (individual) time series data with the aim of characterizing the transition processes underlying the data.The possibility to use covariates on the transition matrix greatly enhances the flexibility of the model.The EM algorithm uses a very general interface that allows easy addition of new response models.
We are currently working on implementing the gradients for response and transition models with two goals in mind.First, to speed up (constrained) parameter optimization using Rdonlp2 or Rsolnp.Second, analytic gradients are useful in computing the Hessian of the estimated parameters so as to arrive at standard errors for those.We are also planning to implement goodness-of-fit statistics (Titman and Sharples 2008).

Figure 1 :
Figure 1: Response times (rt), accuracy (corr) and pay-off values (Pacc) for the first series of responses in dataset speed.

Figure 2 :
Figure 2: Balance scale item; this is a distance item (see the text for details).

Table 1 :
Response distribution available in depmixS4.