How to learn a Hidden Markov Model with categorical responses in R?

Question

I am looking for a mature library to learn hidden markov models with categorical responses, and I want to be able to learn the HMM from several traces. I tried a few options, but I settled for the depmixS4 package.

I can learn a model with multinomial responses, but I do not understand the output I get from summary. Here is what I did:

library(depmixS4)

v = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3)

draws = data.frame(obs = v)

model = depmix(obs ~ 1, data = draws, nstates = 3, family = multinomial())
fm = fit(model)

Then, summary(fm) gives me:

converged at iteration 18 with logLik: -4.612935 
Initial state probabilities model 
pr1 pr2 pr3 
  0   1   0 

Transition matrix 
        toS1  toS2  toS3
fromS1 1.000 0.000 0.000
fromS2 0.000 0.833 0.167
fromS3 0.333 0.000 0.667

Response parameters 
Resp 1 : multinomial 
    Re1.(Intercept).1 Re1.(Intercept).2 Re1.(Intercept).3
St1                 0             4.315            16.664
St2                 0           -13.541           -16.039
St3                 0            14.078             1.221

The transition matrix is more or less what I expect, but the Response parameters aren't. I was expected a matrix with shape n_states x n_observations containing probabilities, and the first column is all zeroes, so this does not look like what I expected.

Now, the official documentation of the package states:

This [the dependent variable] is a binary matrix with N rows and Y columns, where Y is the total number of categories.

But even if I try using one-hot encoding for my variables, I get the same result. Here is the code:

v = c(1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 3)
v_fact = factor(v) # Factor is a list of categorical data
# One hot encoding, thanks to a magic one liner.
one_hot_v = model.matrix(~ 0 + v_fact)
draws = data.frame(one_hot_v)

model = depmix(v_fact ~ 1, data = draws, nstates = 3, family = multinomial())
fm = fit(model)

summary(fm)

Here is the summary.

converged at iteration 21 with logLik: -4.61294 
Initial state probabilities model 
pr1 pr2 pr3 
  0   1   0 

Transition matrix 
        toS1  toS2  toS3
fromS1 0.667 0.000 0.333
fromS2 0.167 0.833 0.000
fromS3 0.000 0.000 1.000

Response parameters 
Resp 1 : multinomial 
    Re1.(Intercept).1 Re1.(Intercept).2 Re1.(Intercept).3
St1                 0            12.598            -8.587
St2                 0           -12.960           -14.685
St3                 0             3.099            15.671

The summary looks the same as in the previous case.

So I have a few questions:

1. I know my responses (or observations) are discrete and categorical, am I on the right track?
2. How do I read these response parameters? What do they mean? (I would understand Intercept if I was doing linear regression, but not in this context).
3. I am interested in the likelihood of observing a value y when the chain is in state s. How can I see the relationship between a state and the response?

Answer 1

1. You are on the right track. A multinomial distribution makes sense here.
2. By default, the submodel for a multinomial response variable will be a multinomial logistic regression model (e.g. https://en.wikipedia.org/wiki/Multinomial_logistic_regression) with baseline coding such that the parameters for the first level of the categorical response are fixed to 0. Hence the values of 0 for "Re1.(Intercept).1" (which stands for "response variable 1"."intercept"."level 1") for all states (rows "St"). When you don't include any predictors/regressors/independent variables in this model (i.e. it is an intercept-only model), more directly interpretable parameters can be obtained by using an identity link function, changing your code slightly: family=multinomial("identity"). The parameters are then simply the probability of each level.
3. The state-dependent density of the response variables is given in the "dens" slot of the (fitted) model object, which is "Array of dimension sum(ntimes)nrespnstates providing the densities of the observed responses for each state." (see ?"depmix-class"). So e.g. fm@dens[,,1] will give you the likelihood of the observations assuming the chain is in state 1.

Hope this helps!