HMM - general concept and strategy

Question

I am new to data science and am trying to deal with an HMM on time series - I ran into some problems and I feel I have to rethink my stategy in general. I would love to get some feedback on my thoughts!

My data: I have time series data with a sequence of measurements each hour. The data have a strong daily rhythm - low values in the night and a peak each day in the afternoon. In addition, I have additional information on the temp, rain etc. at the timepoints.

My goal: To train an HMM on the sequence and to generate synthetic data from the model. It is important, that the circadian rhythm in the original data is also captured by the synthetic data.

What I did so far: I extracted features from the time series (encoded the days as cosinus and sinus waves) and trained a multivariate HMM using the python library pomegranate. Mapping the states on the original time series looks good - some states are characteristic for the nights, some for the morning/afternoon and so on. I also generated sample data from that model.

My problem: The generated synthetic sample data do not show the typcial (and important) daily rhythm/pattern.

My understanding of the problem: An HMM is characterized by the initial state, the emission probability and the transition probability between the hidden states. Each stated n+1 only depends on the state of n. So - trying to nail the problem in very simple words - the model is not able to capture the 24 hour rhythm because it only relies on the ONE previous state.

My thoughts on possible improvements:

Instead of training a sequence of several years, I could slice the time series and retrain the model on hundreds of days. Then, I would generate days and add them to a sequence of the length needed.

Semi hidden markov model: I read that these exist to make sure, that one state can last over several time points. As I understand, this would not help to keep the 24 rhythm.

Harmonic HMM: I saw this in a paper (https://royalsocietypublishing.org/doi/10.1098/rsif.2017.0885). I don´t understand how it would be implemented… but as I understand, the purpose of it was to include a circadian oscillator in the transition matrix.

My question: Could you please give some feedback on how I could proceed? I would greatly appreciate any thoughts, ideas, explanations. Switching to another model is not an option - this problem is part of a bigger project where different approaches should be compared in the end. So my priority is to come up with the best solution that one can get using HMMs.

Thank you very much in advance!!

Answer 1

As you correctly diagnosed, the problem is that a simple HMM does not have any built-in mechanism to account for the daily cycle. One solution, used by Huang et al. (2018), is to use a non-homogeneous HMM, where the time of day is included as a covariate in the transition probability matrix. That is, you can allow the transition probabilities to depend on the time of day. Presumably, based on how you described your data, the probability of a transition from the "high" state to the "low" state would be highest in the evening, and the probability of a transition from "low" to "high" would be highest in the morning. This would capture the cycle, and should produce more realistic patterns in simulations.

An additional practical concern is that the relationship between transition probabilities and time of day should be modelled using a cyclical function, to make sure that it "loops" at midnight. You could check out the recent paper by Feldmann et al. (2023), which talks about this question for HMMs quite generally. (Detailed references are included at the bottom of the post.)

If you are not opposed to using a different software package, the R package hmmTMB allows for cyclical covariate effects (Michelot et al., 2022). As a disclaimer: I am the developer of hmmTMB. This can be used to reproduce something similar to the analysis of activity data presented by Huang et al. (2018). We can compare the observed data and data simulated from the fitted model to make sure that the cyclical patterns are captured. It seems to be the case here, although it looks like the variance in simulated activity is perhaps larger than in the real data. (Whether or not this is a problem depends on the aim of analysis.) The full code to fit the model and reproduce these plots is included below.

Code

library(hmmTMB)
library(ggplot2)
theme_set(theme_bw())

# Load Huang et al (2018) data from Github
path <- "https://raw.githubusercontent.com/huang1010/HMMS/master/"
data1 <- read.csv(url(paste0(path, "S9.csv")))[,-3]
data2 <- read.csv(url(paste0(path, "S20.csv")))[,-3]

# Format data and add "time of day" (tod) variable
data <- rbind(cbind(ID = "S9", data1), 
              cbind(ID = "S20", data2))
data$time <- as.POSIXct(data$time)
data$tod <- as.numeric(format(data$time, "%H")) + 
    as.numeric(format(data$time, "%M"))/60

# Define hidden state model, with transition probabilities
# modelled as cyclical splines of time of day
hid <- MarkovChain$new(data = data, 
                       n_states = 2, 
                       initial_state = "stationary",
                       formula = ~s(tod, k = 5, bs = "cc"))

# Define observation model for "activity" variable, using
# zero-inflated gamma distribution
dists <- list(activity = "zigamma2")
par0 <- list(activity = list(mean = c(20, 150), 
                             sd = c(20, 40), 
                             z = c(0.1, 0)))
obs <- Observation$new(data = data, 
                       dists = dists,
                       n_states = 2, 
                       par = par0)

# Fix zero inflation parameter to zero in state 2 
# (high activity state) to avoid numerical problems
fixpar <- list(obs = c("activity.z.state2.(Intercept)" = NA))

# Create and fit HMM
hmm <- HMM$new(obs = obs, hid = hid, fixpar = fixpar)
hmm$fit(silent = TRUE)

# Simulate from fitted model, using covariate values from original data
sim <- hmm$simulate(n = nrow(data), data = data)

# Plot observed activity data
ggplot(data, aes(time, activity)) +
    geom_point(size = 0.6) +
    geom_line() +
    facet_wrap("ID", scales = "free_x") +
    coord_cartesian(ylim = c(0, 300)) +
    labs(title = "Observed data")

# Plot simulated activity data
ggplot(foo, aes(time, activity)) +
    geom_point(size = 0.6) +
    geom_line() +
    facet_wrap("ID", scales = "free_x") +
    coord_cartesian(ylim = c(0, 300)) +
    labs(title = "Simulated data")

References

• Feldmann, C. C., Mews, S., Coculla, A., Stanewsky, R., & Langrock, R. (2023). Flexible modelling of diel and other periodic variation in hidden Markov models. Journal of Statistical Theory and Practice, 17(3), 45.

• Huang, Q., Cohen, D., Komarzynski, S., Li, X. M., Innominato, P., Lévi, F., & Finkenstädt, B. (2018). Hidden Markov models for monitoring circadian rhythmicity in telemetric activity data. Journal of The Royal Society Interface, 15(139), 20170885.

• Michelot, T. (2022). hmmTMB: Hidden Markov models with flexible covariate effects in R. arXiv preprint arXiv:2211.14139.