The purpose of this post is to give a workflow that one could follow if your goal is to fit Non-Homogeneous Hidden Markov Models (NHMMs) and to simulate sequences of states. Julie Quinn wrote a series of great posts on fitting Hidden Markov Models (HMMs) here but the goal of this post is to discuss the non-homogeneous counterparts of HMM. NHMMs can be distinguished from the former because they involve non-stationary transition probabilities between states. These dynamic transition probability matrices are conditioned on one or more external covariates that influences transitions between states. One example of the use of NHMMs is to model and simulate precipitation. Precipitation models that simulate without using atmospheric information cannot be expected to perform well on conditions that deviate from those on which the model was fit. Thus using a covariate that shows changes in atmospheric circulation (such as geopotential height), can help capture some of the nonstationarity that a solely precipitation-based model alone could not capture.

I have had a lot of time to work with NHMMs in R; however at the beginning, I found overwhelmingly few examples to work off of, so this post is meant to outline a workflow that I have settled on that you can apply to your own application. First off, there are tons of packages available for HMMs, but very few that handle the dynamic transition matrices of NHMMs. My favorite is depmixS4 found here. This package has all the tools needed to fit your NHMM, but it takes some experimenting to understand the syntax and to piece together the functions to create a full workflow.

First we want to fit a depmix model. In the first line of the code, I am creating a model called *modNHMM*. The data that I am fitting the NHMM with are the first 9 principal components of geopotential height. It is important that you list these components in this syntax and they should be the column titles of your dataset, which is *synoptic.pcs* in my case. The number of states is how many states you want to fit your model with. For a precipitation/streamflow application, you could fit a two-state NHMM to represent wet/dry conditions, or if you are trying to identify multiple regimes, you may have more.

```
library(depmixs4)
modNHMM <- depmix(list(PC1~1,PC2~1,PC3~1, PC4~1, PC5~1, PC6~1, PC7~1, PC8~1,
PC9~1),
nstates = num.states,
family=list(gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),gaussian(),
gaussian()),
ntimes = nrow(synoptic.pcs),
data = data.frame(synoptic.pcs),transition = ~predictor$PC1+predictor$PC2+predictor$PC3+predictor$PC4)
fit.modNHMM.depmix.paleo.check <- fit(modNHMM)
synoptic.state.assignments<- posterior(fit.modHMMs.depmix.paleo.check)$state # state sequence using the Viterbi algorithm
```

Then we choose a distribution for our responses, which I choose to be a Gaussian distribution. The argument, *ntimes*, refers to the length of our time series, which is the number of rows in our dataset. Next we specify the dataset that contains the PC data and the transition which is dictated by our external covariates. In this case, my covariates are in a dataframe called *predictor* and each column corresponds to the 4 covariates (which happen to be PCs of a different dataset) that will influence my transitions. Then we fit the model with our prescribed attributes using the *fit* function. Finally, we want to calculate the Viterbi sequence, which is the most probable path or sequence of states over the period that the model is fit. This step (last line of code) will return a vector of the length of the dataset with each day classified into one of *num.states*.

Now we need to locate the transition probability matrices. We use the following command:

`fit.modNHMM.depmix.paleo.check@transition`

If we were fitting an HMM, we would get one transition probability matrix. However, we get an output that looks like this:

Now instead of having constant values for transitions, we have equations of the form: *Intercept+b1*PC1+b2*PC2+b3*PC3+b4*PC4* where the *b* values are the coefficient values listed in the table. If were were to look at the first block [[1]], this block dictates transitions from State 1 to the other 5 states. The transitions from the states are as follows:

State 1 to State 1: 0+0+0+0+0 (State 1 is used as a reference variable in this case and the probability would be found by subtracting the other probabilities from 1 at the end)

State 1 to State 2: -3.11+-0.22*PC1+0.22*PC2+0.014*PC3-0.13*PC4

and so on. You have to remember to divide each value by the total across the row in order to return probabilities.

Once you have these symbolic equations for the transition probability matrices, you can create a list of matrices which will allow you to simulate sequences of states for new sets of the PC1, PC2, PC3, PC4 covariates. You can get a sense how you might create *n* different transition matrices if you have times series of length *n* of the covariates. Below I am just representing those symbolic equations in code using the *getpars* function to acquire the coefficients and store the resulting daily matrices in a list called *mm*. Depending on the number of covariates or states, you will need to adjust the indices accordingly.

```
n=dim(df)[[1]]
mm<-matrix(list(), 1,n)
for (j in 1:n){
transition_matrix=matrix(, nrow = 5, ncol = 5)
for (i in 6:10){
transition_matrix[1,i-5]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[1,]))
for (i in 6:10){
transition_matrix[1,i-5]=exp(transition_matrix[1,i-5])/denominator
}
for (i in 31:35){
transition_matrix[2,i-30]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[2,]))
for (i in 31:35){
transition_matrix[2,i-30]=exp(transition_matrix[2,i-30])/denominator
}
for (i in 56:60){
transition_matrix[3,i-55]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[3,]))
for (i in 56:60){
transition_matrix[3,i-55]=exp(transition_matrix[3,i-55])/denominator
}
for (i in 81:85){
transition_matrix[4,i-80]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[4,]))
for (i in 81:85){
transition_matrix[4,i-80]=exp(transition_matrix[4,i-80])/denominator
}
for (i in 106:110){
transition_matrix[5,i-105]=getpars(fit.modHMMs.depmix.paleo.check)[i]+(getpars(fit.modHMMs.depmix.paleo.check)[i+5]*df$PC1[j])+ (getpars(fit.modHMMs.depmix.paleo.check)[i+10]*df$PC2[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+15]*df$PC3[j])+(getpars(fit.modHMMs.depmix.paleo.check)[i+20]*df$PC4[j])
}
denominator=sum(exp(transition_matrix[5,]))
for (i in 106:110){
transition_matrix[5,i-105]=exp(transition_matrix[5,i-105])/denominator
}
mm[[j]]=transition_matrix
}
```

Once we have these matrices, we can then simulate state sequences that can result from the chain of transition matrices. For this part, we need to create markov lists with our transition matrices:

```
library(markovchain)
mcObject <- mclapply(X=1:iter,mc.preschedule=TRUE,mc.cores=1,FUN=function(j){
mcObject.time.varying <- mclapply(X=1:n.sim,mc.preschedule=TRUE,mc.cores=1,FUN=function(t){
tr.prob=as.matrix(mm[[t]])
mcObject.time.varying.out <- new("markovchain", states = c("1","2","3","4","5"),
transitionMatrix = tr.prob, name = paste("McObject",t,sep=""))
return(McObject.time.varying.out)
}
)
mcObject.final <- new("markovchainList",markovchains = mcObject.time.varying, name = "mcObject.nh")
return(
mcObject.final
)
}
```

Finally we simulate using the following:

```
simulate.mc <- function(mcObject,num.states,dates.sim,last.month,last.day,n.sim,iter) {
#this function will simulate the Markov chain iter times
#Arguments:
#mcObject = a Markov chain object from the markovchain package
#num.states = the number of states
#date.sim = a time series of dates for the simulation period
#last.month = last month of the season
#last.day = last day of the last month of the season
#iter = the number of iterations
#day and month sequences for the simulation period
days.sim <- as.numeric(format(dates.sim,"%d"))
months.sim <- as.numeric(format(dates.sim,"%m"))
n.sim <- length(dates.sim) #length of simulation
final.mc.sim <- mclapply(X=1:iter,mc.preschedule=TRUE,mc.cores=1,FUN=function(i){
mc.sim <- as.numeric(rmarkovchain(n=1,object=mcObject[[i]][[1]]))
end.state <- paste(mc.sim[length(mc.sim)])
for (t in 1:n.sim) {
mc.sim <- c(mc.sim,as.numeric(rmarkovchain(n=1,object=mcObject[[i]][[t]],t0=end.state)))
end.state <- paste(mc.sim[length(mc.sim)])
if(months.sim[t]==last.month & days.sim[t]==last.day) {end.state <- paste(sample(1:num.states,size=1))}
}
#here is the final mc simulation
final.mc.sim.iter <- mc.sim[2:(n.sim+1)]
return(final.mc.sim.iter)
}
)
return(final.mc.sim)
}
simulations=matrix(list(), 1,1000)
for (i in 1:1000){
simulations[[i]] <- simulate.mc(mcObject=mcWeather.Regime,num.states=num.states,
dates.sim=dates.sim,last.month=last.month,last.day=last.day,iter=iter)
}
```

And that’s it! You can simulate for many different iterations (1000 in my case) and you will be returned a large list with your 1000 sequence of states over the simulation period.

This is a great post! Thanks for sharing. I am trying to test my data set by using this code but I got an error (invalid ‘trim’ argument) at the last step (simulations[[1]]<-). My data is yearly based (precipitation and PMDI as covariate) and I am not using principle component analysis prior to fitting the model. Why do you think I am getting this error? Any help will be appreciated! Thank you!

Hi Burcu, I haven’t run into that error before, but it seems like that error is a formatting one. I would double check that your date is in the form of a Date object. You can convert a numeric/character date sequence into a Date object using the as.Date() function in R. Hope that helps! Otherwise, feel free to deconstruct what I have into something that suits your application. The key is really just to simulate using rmarkovchain once you have the markov chain lists from the prior block of code. You can create a different piece of code that will line up the dates for you.

Hello Rohini

This is a great post, and extremely helpful. I am so glad you shared this. I am running into a issue, in which sometimes the fitting the model gives me an error – initial value in ‘vmmin’ is not finite and sometimes it says- initial values not feasible. How do I address these issues?

Thank you so much for your help.

Hi @sonam, this is not an error I have encountered before, but it looks to be related to the fitting which utilizes a likelihood approach. For example, 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. It looks like the initial value when the log-likelihood of the parameters might lead to an undefined/infinite likelihood. This may occur due an unstable fit- so if your number of states is not ideal or covariates are not meaningful. There may be a way to manually change the initial condition, or you can consider changing your formulation to something a bit more well-posed.

Great post! A few points for all:

1) HMM models are notorious for getting stuck in local maxima. Irrespective of R package and model configuration, it is better to fit a few (like 10) models with random starting points (by just changing seed() in R, for example). Then one should pick the model with some set criteria, for example, lowest BIC.

2) Directly simulate N traces with starting seed K from fitted model using depmixS4::simulate(fit.modHMMs.depmix.paleo.check, nsim = N ,seed = K)

3) One can avoid the for loop to get parameters of transition model by: lapply(fit.modHMMs.depmix.paleo.check@transition, function(x) x@parameters$coefficients)

You are the true expert on this Sudarshana! Thanks for sharing your comments here!