Title:	Mixture Hidden Markov Models for Social Sequence Data and Other Multivariate, Multichannel Categorical Time Series
Version:	2.0.0
Description:	Designed for estimating variants of hidden (latent) Markov models (HMMs), mixture HMMs, and non-homogeneous HMMs (NHMMs) for social sequence data and other categorical time series. Special cases include feedback-augmented NHMMs, Markov models without latent layer, mixture Markov models, and latent class models. The package supports models for one or multiple subjects with one or multiple parallel sequences (channels). External covariates can be added to explain cluster membership in mixture models as well as initial, transition and emission probabilities in NHMMs. The package provides functions for evaluating and comparing models, as well as functions for visualizing of multichannel sequence data and HMMs. For NHMMs, methods for computing average causal effects and marginal state and emission probabilities are available. Models are estimated using maximum likelihood via the EM algorithm or direct numerical maximization with analytical gradients. Documentation is available via several vignettes, and Helske and Helske (2019, <doi:10.18637/jss.v088.i03>). For methodology behind the NHMMs, see Helske (2025, <doi:10.48550/arXiv.2503.16014>).
LazyData:	true
LinkingTo:	nloptr, Rcpp (≥ 0.12.0), RcppArmadillo
Depends:	R (≥ 4.1.0)
Imports:	checkmate, cli, data.table, future.apply, ggplot2, ggseqplot, graphics, grDevices, grid, gridBase, igraph, lhs, Matrix, methods, nloptr, numDeriv, patchwork, progressr, Rcpp (≥ 0.12.0), RcppHungarian, rlang, stats, TraMineR (≥ 2.2-7), utils
Suggests:	covr, knitr, MASS, nnet, testthat (≥ 3.0.0),
License:	GPL-2 | GPL-3 [expanded from: GPL (≥ 2)]
Encoding:	UTF-8
BugReports:	https://github.com/helske/seqHMM/issues
VignetteBuilder:	knitr
RoxygenNote:	7.3.2
Config/testthat/edition:	3
SystemRequirements:	GNU make
Biarch:	true
NeedsCompilation:	yes
Packaged:	2025-05-16 11:54:34 UTC; jvhels
Author:	Jouni Helske [aut, cre], Satu Helske [aut]
Maintainer:	Jouni Helske <jouni.helske@iki.fi>
Repository:	CRAN
Date/Publication:	2025-05-17 00:10:02 UTC

The seqHMM package

Description

The seqHMM package is designed for fitting hidden (or latent) Markov models (HMMs) and mixture hidden Markov models (MHMMs) for social sequence data and other categorical time series. The package supports models for one or multiple subjects with one or multiple interdependent sequences (channels). External covariates can be added to explain cluster membership in mixture models. The package provides functions for evaluating and comparing models, as well as functions for easy plotting of multichannel sequences and hidden Markov models. Common restricted versions of (M)HMMs are also supported, namely Markov models, mixture Markov models, and latent class models.

Details

Maximum likelihood estimation via the EM algorithm and direct numerical maximization with analytical gradients is supported. All main algorithms are written in C++. Parallel computation is implemented via OpenMP for pre-2.0.0 functions, while estimation of non-homogenous models support parallelization via future package by parallelization of multistart optimizations and bootstrap sampling.

Author(s)

Maintainer: Jouni Helske jouni.helske@iki.fi (ORCID)

Authors:

• Satu Helske (ORCID)

References

Helske S. and Helske J. (2019). Mixture Hidden Markov Models for Sequence Data: The seqHMM Package in R, Journal of Statistical Software, 88(3), 1-32. doi:10.18637/jss.v088.i03

See Also

Useful links:

• Report bugs at https://github.com/helske/seqHMM/issues

Three-channel biofam data

Description

Biofam data from the TraMineR package converted into three channels.

Format

A list including three sequence data sets for 2000 individuals with 16 state variables, and a separate data frame with 1 id variable, 8 covariates, and 2 weight variables.

Details

This data is constructed from the TraMineR::biofam() data in the TraMineR package. Here the original state sequences are converted into three separate data sets: children, married, and left. These include the corresponding life states from age 15 to 30: childless or (having) children; single, married, or divorced; and (living) ⁠with parents⁠ or ⁠left home⁠.

Note that the divorced state does not give information on parenthood or residence, so a guess is made based on preceeding states.

The fourth data frame covariates is a collection of additional variables from the original data:

The data is loaded by calling data(biofam3c). It was built using following code:

data("biofam" , package = "TraMineR")
biofam3c <- with(biofam, {

## Building one channel per type of event left, children or married
bf <- as.matrix(biofam[, 10:25])
children <- bf == 4 | bf == 5 | bf == 6
married <- bf == 2 | bf == 3 | bf == 6
left <- bf == 1 | bf == 3 | bf == 5 | bf == 6 | bf == 7

children[children == TRUE] <- "children"
children[children == FALSE] <- "childless"
# Divorced parents
div <- bf[(rowSums(bf == 7) > 0 & rowSums(bf == 5) > 0) |
            (rowSums(bf == 7) > 0 & rowSums(bf == 6) > 0),]
children[rownames(bf) %in% rownames(div) & bf == 7] <- "children"

married[married == TRUE] <- "married"
married[married == FALSE] <- "single"
married[bf == 7] <- "divorced"

left[left == TRUE] <- "left home"
left[left == FALSE] <- "with parents"
# Divorced living with parents (before divorce)
wp <- bf[(rowSums(bf == 7) > 0 & rowSums(bf == 2) > 0 &
          rowSums(bf == 3) == 0 & rowSums(bf == 5) == 0 &
          rowSums(bf == 6) == 0) |
         (rowSums(bf == 7) > 0 & rowSums(bf == 4) > 0 &
          rowSums(bf == 3) == 0 & rowSums(bf == 5) == 0 &
          rowSums(bf == 6) == 0), ]
left[rownames(bf) %in% rownames(wp) & bf == 7] <- "with parents"

list("children" = children, "married" = married, "left" = left,
  "covariates" = biofam[, c(1:9, 26:27)])
})

Source

TraMineR::biofam() data constructed from the Swiss Household Panel https://forscenter.ch/projects/swiss-household-panel/

References

Müller, N. S., M. Studer, G. Ritschard (2007). Classification de parcours de vie à l'aide de l'optimal matching. In XIVe Rencontre de l a Société francophone de classification (SFC 2007), Paris, 5 - 7 septembre 2007, pp. 157–160.

Bootstrap Sampling of NHMM Coefficients

Description

It is possible to parallelize the bootstrap runs using the future package, e.g., by calling future::plan(multisession, workers = 2) before bootstrap_coefs(). See future::plan() for details.

Usage

bootstrap_coefs(model, ...)

## S3 method for class 'nhmm'
bootstrap_coefs(
  model,
  nsim,
  type = c("nonparametric", "parametric"),
  append = FALSE,
  ...
)

## S3 method for class 'mnhmm'
bootstrap_coefs(
  model,
  nsim,
  type = c("nonparametric", "parametric"),
  append = FALSE,
  ...
)

Arguments

Details

bootstrap_coefs() is compatible with progressr package, so you can use progressr::with_progress(bootstrap_coefs(fit)) to track the progress of bootstrapping.

Value

The original model with additional element model$boot.

Build a Hidden Markov Model

Description

Function build_hmm constructs a hidden Markov model object of class hmm.

Usage

build_hmm(
  observations,
  n_states,
  transition_probs,
  emission_probs,
  initial_probs,
  state_names = NULL,
  channel_names = NULL,
  ...
)

Arguments

Details

The returned model contains some attributes such as nobs and df, which define the number of observations in the model and the number of estimable model parameters, used in computing BIC. When computing nobs for a multichannel model with C channels, each observed value in a single channel amounts to 1/C observation, i.e. a fully observed time point for a single sequence amounts to one observation. For the degrees of freedom df, zero probabilities of the initial model are defined as structural zeroes.

Value

Object of class hmm with the following elements:

• observations
  State sequence object or a list of such objects containing the data.

• transition_probs
  A matrix of transition probabilities.

• emission_probs
  A matrix or a list of matrices of emission probabilities.

• initial_probs
  A vector of initial probabilities.

• state_names
  Names for hidden states.

• symbol_names
  Names for observed states.

• channel_names
  Names for channels of sequence data.

• length_of_sequences
  (Maximum) length of sequences.

• sequence_lengths
  A vector of sequence lengths.

• n_sequences
  Number of sequences.

• n_symbols
  Number of observed states (in each channel).

• n_states
  Number of hidden states.

• n_channels
  Number of channels.

See Also

fit_model() for estimating model parameters; and plot.hmm() for plotting hmm objects.

Examples

# Single-channel data

data("mvad", package = "TraMineR")

mvad_alphabet <- c(
  "employment", "FE", "HE", "joblessness", "school",
  "training"
)
mvad_labels <- c(
  "employment", "further education", "higher education",
  "joblessness", "school", "training"
)
mvad_scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad_seq <- seqdef(mvad, 15:86,
  alphabet = mvad_alphabet, states = mvad_scodes,
  labels = mvad_labels, xtstep = 6
)

# Initializing an HMM with 4 hidden states, random starting values
init_hmm_mvad1 <- build_hmm(observations = mvad_seq, n_states = 4)

# Starting values for the emission matrix
emiss <- matrix(NA, nrow = 4, ncol = 6)
emiss[1, ] <- seqstatf(mvad_seq[, 1:12])[, 2] + 1
emiss[2, ] <- seqstatf(mvad_seq[, 13:24])[, 2] + 1
emiss[3, ] <- seqstatf(mvad_seq[, 25:48])[, 2] + 1
emiss[4, ] <- seqstatf(mvad_seq[, 49:70])[, 2] + 1
emiss <- emiss / rowSums(emiss)

# Starting values for the transition matrix

tr <- matrix(
  c(
    0.80, 0.10, 0.05, 0.05,
    0.05, 0.80, 0.10, 0.05,
    0.05, 0.05, 0.80, 0.10,
    0.05, 0.05, 0.10, 0.80
  ),
  nrow = 4, ncol = 4, byrow = TRUE
)

# Starting values for initial state probabilities
init <- c(0.3, 0.3, 0.2, 0.2)

# HMM with own starting values
init_hmm_mvad2 <- build_hmm(
  observations = mvad_seq, transition_probs = tr,
  emission_probs = emiss, initial_probs = init
)

#########################################


# Multichannel data

# Three-state three-channel hidden Markov model
# See ?hmm_biofam for a five-state version

data("biofam3c")

# Building sequence objects
marr_seq <- seqdef(biofam3c$married,
  start = 15,
  alphabet = c("single", "married", "divorced"),
  cpal = c("violetred2", "darkgoldenrod2", "darkmagenta")
)
child_seq <- seqdef(biofam3c$children,
  start = 15,
  alphabet = c("childless", "children"),
  cpal = c("darkseagreen1", "coral3")
)
left_seq <- seqdef(biofam3c$left,
  start = 15,
  alphabet = c("with parents", "left home"),
  cpal = c("lightblue", "red3")
)

# You could also define the colors using cpal function from TraMineR
# cpal(marr_seq) <- c("violetred2", "darkgoldenrod2", "darkmagenta")
# cpal(child_seq) <- c("darkseagreen1", "coral3")
# cpal(left_seq) <- c("lightblue", "red3")

# Left-to-right HMM with 3 hidden states and random starting values
set.seed(1010)
init_hmm_bf1 <- build_hmm(
  observations = list(marr_seq, child_seq, left_seq),
  n_states = 3, left_right = TRUE, diag_c = 2
)


# Starting values for emission matrices

emiss_marr <- matrix(NA, nrow = 3, ncol = 3)
emiss_marr[1, ] <- seqstatf(marr_seq[, 1:5])[, 2] + 1
emiss_marr[2, ] <- seqstatf(marr_seq[, 6:10])[, 2] + 1
emiss_marr[3, ] <- seqstatf(marr_seq[, 11:16])[, 2] + 1
emiss_marr <- emiss_marr / rowSums(emiss_marr)

emiss_child <- matrix(NA, nrow = 3, ncol = 2)
emiss_child[1, ] <- seqstatf(child_seq[, 1:5])[, 2] + 1
emiss_child[2, ] <- seqstatf(child_seq[, 6:10])[, 2] + 1
emiss_child[3, ] <- seqstatf(child_seq[, 11:16])[, 2] + 1
emiss_child <- emiss_child / rowSums(emiss_child)

emiss_left <- matrix(NA, nrow = 3, ncol = 2)
emiss_left[1, ] <- seqstatf(left_seq[, 1:5])[, 2] + 1
emiss_left[2, ] <- seqstatf(left_seq[, 6:10])[, 2] + 1
emiss_left[3, ] <- seqstatf(left_seq[, 11:16])[, 2] + 1
emiss_left <- emiss_left / rowSums(emiss_left)

# Starting values for transition matrix
trans <- matrix(
  c(
    0.9, 0.07, 0.03,
    0, 0.9, 0.1,
    0, 0, 1
  ),
  nrow = 3, ncol = 3, byrow = TRUE
)

# Starting values for initial state probabilities
inits <- c(0.9, 0.09, 0.01)

# HMM with own starting values
init_hmm_bf2 <- build_hmm(
  observations = list(marr_seq, child_seq, left_seq),
  transition_probs = trans,
  emission_probs = list(emiss_marr, emiss_child, emiss_left),
  initial_probs = inits
)

Build a Latent Class Model

Description

Function build_lcm is a shortcut for constructing a latent class model as a restricted case of an mhmm object.

Usage

build_lcm(
  observations,
  n_clusters,
  emission_probs,
  formula = NULL,
  data = NULL,
  coefficients = NULL,
  cluster_names = NULL,
  channel_names = NULL
)

Arguments

Value

Object of class mhmm with the following elements:

• observations
  State sequence object or a list of such containing the data.

• transition_probs
  A matrix of transition probabilities.

• emission_probs
  A matrix or a list of matrices of emission probabilities.

• initial_probs
  A vector of initial probabilities.

• coefficients
  A matrix of parameter coefficients for covariates (covariates in rows, clusters in columns).

• X
  Covariate values for each subject.

• cluster_names
  Names for clusters.

• state_names
  Names for hidden states.

• symbol_names
  Names for observed states.

• channel_names
  Names for channels of sequence data

• length_of_sequences
  (Maximum) length of sequences.

• sequence_lengths
  A vector of sequence lengths.

• n_sequences
  Number of sequences.

• n_symbols
  Number of observed states (in each channel).

• n_states
  Number of hidden states.

• n_channels
  Number of channels.

• n_covariates
  Number of covariates.

• n_clusters
  Number of clusters.

See Also

fit_model() for estimating model parameters; summary.mhmm() for a summary of a mixture model; separate_mhmm() for organizing an mhmm object into a list of separate hmm objects; and plot.mhmm() for plotting mixture models.

Examples

# Simulate observations from two classes
set.seed(123)
obs <- seqdef(rbind(
  matrix(sample(letters[1:3], 500, TRUE, prob = c(0.1, 0.6, 0.3)), 50, 10),
  matrix(sample(letters[1:3], 200, TRUE, prob = c(0.4, 0.4, 0.2)), 20, 10)
))

# Initialize the model
set.seed(9087)
model <- build_lcm(obs, n_clusters = 2)

# Estimate model parameters
fit <- fit_model(model)

# How many of the observations were correctly classified:
sum(summary(fit$model)$most_probable_cluster == rep(c("Class 2", "Class 1"), 
  times = c(500, 200)))

############################################################
## Not run: 
# LCM for longitudinal data

# Define sequence data
data("mvad", package = "TraMineR")
mvad_alphabet <- c(
  "employment", "FE", "HE", "joblessness", "school",
  "training"
)
mvad_labels <- c(
  "employment", "further education", "higher education",
  "joblessness", "school", "training"
)
mvad_scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad_seq <- seqdef(mvad, 15:86,
  alphabet = mvad_alphabet, states = mvad_scodes,
  labels = mvad_labels, xtstep = 6
)

# Initialize the LCM with random starting values
set.seed(7654)
init_lcm_mvad1 <- build_lcm(
  observations = mvad_seq,
  n_clusters = 2, formula = ~male, data = mvad
)

# Own starting values for emission probabilities
emiss <- rbind(rep(1 / 6, 6), rep(1 / 6, 6))

# LCM with own starting values
init_lcm_mvad2 <- build_lcm(
  observations = mvad_seq,
  emission_probs = emiss, formula = ~male, data = mvad
)

# Estimate model parameters (EM algorithm with random restarts)
lcm_mvad <- fit_model(init_lcm_mvad1,
  control_em = list(restart = list(times = 5))
)$model

# Plot the LCM
plot(lcm_mvad, interactive = FALSE, ncol = 2)

###################################################################

# Binomial regression (comparison to glm)

require("MASS")
data("birthwt")

model <- build_lcm(
  observations = seqdef(birthwt$low), emission_probs = diag(2),
  formula = ~ age + lwt + smoke + ht, data = birthwt
)
fit <- fit_model(model)
summary(fit$model)
summary(glm(low ~ age + lwt + smoke + ht, binomial, data = birthwt))


# Multinomial regression (comparison to multinom)

require("nnet")

set.seed(123)
n <- 100
X <- cbind(1, x1 = runif(n, 0, 1), x2 = runif(n, 0, 1))
coefs <- cbind(0, c(-2, 5, -2), c(0, -2, 2))
pr <- exp(X %*% coefs) + stats::rnorm(n * 3)
pr <- pr / rowSums(pr)
y <- apply(pr, 1, which.max)
table(y)

model <- build_lcm(
  observations = seqdef(y), emission_probs = diag(3),
  formula = ~ x1 + x2, data = data.frame(X[, -1])
)
fit <- fit_model(model)
summary(fit$model)
summary(multinom(y ~ x1 + x2, data = data.frame(X[, -1])))

## End(Not run)

Build a Mixture Hidden Markov Model

Description

Function build_mhmm constructs a mixture hidden Markov model object of class mhmm.

Usage

build_mhmm(
  observations,
  n_states,
  transition_probs,
  emission_probs,
  initial_probs,
  formula = NULL,
  data = NULL,
  coefficients = NULL,
  cluster_names = NULL,
  state_names = NULL,
  channel_names = NULL,
  ...
)

Arguments

Details

The returned model contains some attributes such as nobs and df, which define the number of observations in the model and the number of estimable model parameters, used in computing BIC. When computing nobs for a multichannel model with C channels, each observed value in a single channel amounts to 1/C observation, i.e. a fully observed time point for a single sequence amounts to one observation. For the degrees of freedom df, zero probabilities of the initial model are defined as structural zeroes.

Value

Object of class mhmm with following elements:

• observations
  State sequence object or a list of such containing the data.

• transition_probs
  A matrix of transition probabilities.

• emission_probs
  A matrix or a list of matrices of emission probabilities.

• initial_probs
  A vector of initial probabilities.

• coefficients
  A matrix of parameter coefficients for covariates (covariates in rows, clusters in columns).

• X
  Covariate values for each subject.

• cluster_names
  Names for clusters.

• state_names
  Names for hidden states.

• symbol_names
  Names for observed states.

• channel_names
  Names for channels of sequence data

• length_of_sequences
  (Maximum) length of sequences.

• sequence_lengths
  A vector of sequence lengths.

• n_sequences
  Number of sequences.

• n_symbols
  Number of observed states (in each channel).

• n_states
  Number of hidden states.

• n_channels
  Number of channels.

• n_covariates
  Number of covariates.

• n_clusters
  Number of clusters.

References

Helske S. and Helske J. (2019). Mixture Hidden Markov Models for Sequence Data: The seqHMM Package in R, Journal of Statistical Software, 88(3), 1-32. doi:10.18637/jss.v088.i03

See Also

fit_model() for fitting mixture Hidden Markov models; summary.mhmm() for a summary of a MHMM; separate_mhmm() for reorganizing a MHMM into a list of separate hidden Markov models; and plot.mhmm() for plotting mhmm objects.

Examples

data("biofam3c")

## Building sequence objects
marr_seq <- seqdef(biofam3c$married,
  start = 15,
  alphabet = c("single", "married", "divorced"),
  cpal = c("#AB82FF", "#E6AB02", "#E7298A")
)
child_seq <- seqdef(biofam3c$children,
  start = 15,
  alphabet = c("childless", "children"),
  cpal = c("#66C2A5", "#FC8D62")
)
left_seq <- seqdef(biofam3c$left,
  start = 15,
  alphabet = c("with parents", "left home"),
  cpal = c("#A6CEE3", "#E31A1C")
)

## MHMM with random starting values, no covariates
set.seed(468)
init_mhmm_bf1 <- build_mhmm(
  observations = list(marr_seq, child_seq, left_seq),
  n_states = c(4, 4, 6),
  channel_names = c("Marriage", "Parenthood", "Residence")
)


## Starting values for emission probabilities

# Cluster 1
B1_marr <- matrix(
  c(
    0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.3, 0.6, 0.1, # High probability for married
    0.3, 0.3, 0.4
  ), # High probability for divorced
  nrow = 4, ncol = 3, byrow = TRUE
)

B1_child <- matrix(
  c(
    0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.9, 0.1,
    0.9, 0.1
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

B1_left <- matrix(
  c(
    0.9, 0.1, # High probability for living with parents
    0.1, 0.9, # High probability for having left home
    0.1, 0.9,
    0.1, 0.9
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

# Cluster 2

B2_marr <- matrix(
  c(
    0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.1, 0.8, 0.1, # High probability for married
    0.7, 0.2, 0.1
  ),
  nrow = 4, ncol = 3, byrow = TRUE
)

B2_child <- matrix(
  c(
    0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.9, 0.1,
    0.1, 0.9
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

B2_left <- matrix(
  c(
    0.9, 0.1, # High probability for living with parents
    0.1, 0.9,
    0.1, 0.9,
    0.1, 0.9
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

# Cluster 3
B3_marr <- matrix(
  c(
    0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.8, 0.1, 0.1,
    0.1, 0.8, 0.1, # High probability for married
    0.3, 0.4, 0.3,
    0.1, 0.1, 0.8
  ), # High probability for divorced
  nrow = 6, ncol = 3, byrow = TRUE
)

B3_child <- matrix(
  c(
    0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.5, 0.5,
    0.5, 0.5,
    0.5, 0.5,
    0.1, 0.9
  ),
  nrow = 6, ncol = 2, byrow = TRUE
)


B3_left <- matrix(
  c(
    0.9, 0.1, # High probability for living with parents
    0.1, 0.9,
    0.5, 0.5,
    0.5, 0.5,
    0.1, 0.9,
    0.1, 0.9
  ),
  nrow = 6, ncol = 2, byrow = TRUE
)

# Starting values for transition matrices
A1 <- matrix(
  c(
    0.80, 0.16, 0.03, 0.01,
    0, 0.90, 0.07, 0.03,
    0, 0, 0.90, 0.10,
    0, 0, 0, 1
  ),
  nrow = 4, ncol = 4, byrow = TRUE
)

A2 <- matrix(
  c(
    0.80, 0.10, 0.05, 0.03, 0.01, 0.01,
    0, 0.70, 0.10, 0.10, 0.05, 0.05,
    0, 0, 0.85, 0.01, 0.10, 0.04,
    0, 0, 0, 0.90, 0.05, 0.05,
    0, 0, 0, 0, 0.90, 0.10,
    0, 0, 0, 0, 0, 1
  ),
  nrow = 6, ncol = 6, byrow = TRUE
)

# Starting values for initial state probabilities
initial_probs1 <- c(0.9, 0.07, 0.02, 0.01)
initial_probs2 <- c(0.9, 0.04, 0.03, 0.01, 0.01, 0.01)

# Birth cohort
biofam3c$covariates$cohort <- cut(biofam3c$covariates$birthyr, c(1908, 1935, 1945, 1957))
biofam3c$covariates$cohort <- factor(
  biofam3c$covariates$cohort,
  labels = c("1909-1935", "1936-1945", "1946-1957")
)

## MHMM with own starting values and covariates
init_mhmm_bf2 <- build_mhmm(
  observations = list(marr_seq, child_seq, left_seq),
  initial_probs = list(initial_probs1, initial_probs1, initial_probs2),
  transition_probs = list(A1, A1, A2),
  emission_probs = list(
    list(B1_marr, B1_child, B1_left),
    list(B2_marr, B2_child, B2_left),
    list(B3_marr, B3_child, B3_left)
  ),
  formula = ~ sex + cohort, data = biofam3c$covariates,
  cluster_names = c("Cluster 1", "Cluster 2", "Cluster 3"),
  channel_names = c("Marriage", "Parenthood", "Residence"),
  state_names = list(
    paste("State", 1:4), paste("State", 1:4),
    paste("State", 1:6)
  )
)

Build a Markov Model

Description

Function build_mm() builds and automatically estimates a Markov model. It is also a shortcut for constructing a Markov model as a restricted case of an hmm object.

Usage

build_mm(observations)

Arguments

Details

Unlike the other build functions in seqHMM, the build_mm() function automatically estimates the model parameters. In case of no missing values, initial and transition probabilities are directly estimated from the observed initial state probabilities and transition counts. In case of missing values, the EM algorithm is run once.

Note that it is possible that the data contains a symbol from which there are no transitions anywhere (even to itself), which would lead to a row in transition matrix full of zeros. In this case the build_mm() (as well as the EM algorithm) assumes that the the state is absorbing in a way that probability of staying in this state is 1.

Value

Object of class hmm with following elements:

• observations
  State sequence object or a list of such containing the data.

• transition_probs
  A matrix of transition probabilities.

• emission_probs
  A matrix or a list of matrices of emission probabilities.

• initial_probs
  A vector of initial probabilities.

• state_names
  Names for hidden states.

• symbol_names
  Names for observed states.

• channel_names
  Names for channels of sequence data

• length_of_sequences
  (Maximum) length of sequences.

• sequence_lengths
  A vector of sequence lengths.

• n_sequences
  Number of sequences.

• n_symbols
  Number of observed states (in each channel).

• n_states
  Number of hidden states.

• n_channels
  Number of channels.

See Also

plot.hmm() for plotting the model.

Examples

# Construct sequence data
data("mvad", package = "TraMineR")

mvad_alphabet <-
  c("employment", "FE", "HE", "joblessness", "school", "training")
mvad_labels <- c(
  "employment", "further education", "higher education",
  "joblessness", "school", "training"
)
mvad_scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad_seq <- seqdef(mvad, 15:86,
  alphabet = mvad_alphabet,
  states = mvad_scodes, labels = mvad_labels, xtstep = 6,
  cpal = colorpalette[[6]]
)

# Estimate the Markov model
mm_mvad <- build_mm(observations = mvad_seq)

Build a Mixture Markov Model

Description

Function build_mmm() is a shortcut for constructing a mixture Markov model as a restricted case of an mhmm object.

Usage

build_mmm(
  observations,
  n_clusters,
  transition_probs,
  initial_probs,
  formula = NULL,
  data = NULL,
  coefficients = NULL,
  cluster_names = NULL,
  ...
)

Arguments

Value

Object of class mhmm with following elements:

• observations
  State sequence object or a list of such containing the data.

• transition_probs
  A matrix of transition probabilities.

• emission_probs
  A matrix or a list of matrices of emission probabilities.

• initial_probs
  A vector of initial probabilities.

• coefficients
  A matrix of parameter coefficients for covariates (covariates in rows, clusters in columns).

• X
  Covariate values for each subject.

• cluster_names
  Names for clusters.

• state_names
  Names for hidden states.

• symbol_names
  Names for observed states.

• channel_names
  Names for channels of sequence data

• length_of_sequences
  (Maximum) length of sequences.

• sequence_lengths
  A vector of sequence lengths.

• n_sequences
  Number of sequences.

• n_symbols
  Number of observed states (in each channel).

• n_states
  Number of hidden states.

• n_channels
  Number of channels.

• n_covariates
  Number of covariates.

• n_clusters
  Number of clusters.

See Also

fit_model() for estimating model parameters; summary.mhmm() for a summary of a mixture model; separate_mhmm() for organizing an mhmm object into a list of separate hmm objects; and plot.mhmm() for plotting mixture models.

Examples

# Define sequence data
data("mvad", package = "TraMineR")
mvad_alphabet <- c(
  "employment", "FE", "HE", "joblessness", "school",
  "training"
)
mvad_labels <- c(
  "employment", "further education", "higher education",
  "joblessness", "school", "training"
)
mvad_scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad_seq <- seqdef(mvad, 15:86,
  alphabet = mvad_alphabet, states = mvad_scodes,
  labels = mvad_labels, xtstep = 6
)

# Initialize the MMM
set.seed(123)
mmm_mvad <- build_mmm(
  observations = mvad_seq,
  n_clusters = 2,
  formula = ~male, data = mvad
)

## Not run: 
# Estimate model parameters
mmm_mvad <- fit_model(mmm_mvad)$model

# Plot model (both clusters in the same plot)
require(igraph)
plot(mmm_mvad,
  interactive = FALSE,
  # Modify legend position and properties
  with.legend = "right", legend.prop = 0.3, cex.legend = 1.2,
  # Define vertex layout
  layout = layout_as_star,
  # Modify edge properties
  edge.label = NA, edge.arrow.size = 0.8, edge.curved = 0.2,
  # Modify vertex label positions (initial probabilities)
  vertex.label.pos = c("left", "right", "right", "left", "left", "right")
)

# Summary of the MMM
summary(mmm_mvad)

## End(Not run)

Get Cluster Names from Mixture HMMs

Description

Get Cluster Names from Mixture HMMs

Usage

cluster_names(object)

Arguments

Value

A character vector containing the cluster names.

Set Cluster Names for Mixture Models

Description

Set Cluster Names for Mixture Models

Usage

cluster_names(object) <- value

Arguments

Value

The modified object with updated cluster names.

Get the Estimated Regression Coefficients of Non-Homogeneous Hidden Markov Models

Description

Get the Estimated Regression Coefficients of Non-Homogeneous Hidden Markov Models

Usage

## S3 method for class 'nhmm'
coef(object, probs = NULL, ...)

## S3 method for class 'mnhmm'
coef(object, probs = NULL, ...)

Arguments

Value

A list of data tables with the estimated coefficients for initial, transition, emission (separate data.table for each response), and cluster probabilities (in case of mixture model).

Color palettes

Description

A list containing ready defined color palettes with distinct colors using iWantHue. By default, seqHMM uses these palettes when assigning colors.

Format

A list with 200 color palettes.

Source

iWantHue web page https://medialab.github.io/iwanthue/

See Also

plot_colors() for visualization of color palettes. Implementations of iWantHue for R:

• https://github.com/hoesler/rwantshue

• https://github.com/johnbaums/hues

Examples

data("colorpalette")
# Color palette with 9 colors
colorpalette[[9]]
# Color palette with 24 colors
colorpalette[[24]]

Transform TraMineR's state sequence object to data.table and vice versa

Description

Transform TraMineR's state sequence object to data.table and vice versa

Usage

data_to_stslist(x, id, time, responses, seqdef_args = NULL, ...)

stslist_to_data(x, id, time, responses, ...)

Arguments

Estimate a Mixture Non-homogeneous Hidden Markov Model

Description

Function estimate_mnhmm estimates a mixture version of non-homogeneous hidden Markov model (MNHMM) where initial, transition, emission, and mixture probabilities can depend on covariates. See estimate_nhmm() for further details.

Usage

estimate_mnhmm(
  n_states,
  n_clusters,
  emission_formula,
  initial_formula = ~1,
  transition_formula = ~1,
  cluster_formula = ~1,
  data,
  time,
  id,
  lambda = 0,
  prior_obs = "fixed",
  state_names = NULL,
  cluster_names = NULL,
  inits = "random",
  init_sd = 2,
  restarts = 0L,
  method = "EM-DNM",
  bound = Inf,
  control_restart = list(),
  control_mstep = list(),
  ...
)

Arguments

Value

Object of class mnhmm.

See Also

estimate_nhmm() for further details.

Examples

data("mvad", package = "TraMineR")

d <- reshape(mvad, direction = "long", varying = list(15:86), 
  v.names = "activity")

## Not run: 
set.seed(1)
fit <- estimate_mnhmm(n_states = 3, n_clusters = 2,
  data = d, time = "time", id = "id", 
  cluster_formula = ~ male + catholic + gcse5eq + Grammar + 
    funemp + fmpr + livboth + Belfast +
  N.Eastern + Southern + S.Eastern + Western,
  emission_formula = activity ~ male + catholic + gcse5eq,
  initial_formula = ~ 1, 
  transition_formula = ~ male + gcse5eq
  )

## End(Not run)

Estimate a Non-homogeneous Hidden Markov Model

Description

Function estimate_nhmm estimates a non-homogeneous hidden Markov model (NHMM) where initial, transition, and emission probabilities can depend on covariates. Transition and emission probabilities can also depend on past responses, in which case the model is called feedback-augmented NHMM (FAN-HMM) (Helske, 2025).

Usage

estimate_nhmm(
  n_states,
  emission_formula,
  initial_formula = ~1,
  transition_formula = ~1,
  data,
  time,
  id,
  lambda = 0,
  prior_obs = "fixed",
  state_names = NULL,
  inits = "random",
  init_sd = 2,
  restarts = 0L,
  method = "EM-DNM",
  bound = Inf,
  control_restart = list(),
  control_mstep = list(),
  ...
)

Arguments

Details

In case of FAN-HMM with autoregressive dependency on the observational level, (i.e. response y_t depend on y_{t-1}), the emission probabilities at the first time point need special attention. By default, the model is initialized with fixed values for the first time point (prior_obs = "fixed"), meaning that if the input data consists of time points t=1, 2, \ldots, then the model is defined from t=2 onwards and the data on t=1 is used only for defining the emission probabilities at t=2. Note that in this case also the initial state probabilities correspond to t=2.

Alternatively, you can define prior_obs as a list of vectors, where the number of vectors is equal to the number of responses, and each vector gives the prior distribution for the response at t=0. For example, if you have response variables y and x, where y has 3 categories and x 2 categories, you can define prior_obs = list(y = c(0.5, 0.3, 0.2), x = c(0.7, 0.3)). These distributions are then used to marginalize out y_0 and x_0 in the relevant emission probabilities.

By default, the model parameters are estimated using EM-DNM algorithm which first runs some iterations (100 by default) of EM algorithm, and then switches to L-BFGS. Other options include any numerical optimization algorithm of nloptr::nloptr(), or plain EM algorithm where the M-step uses L-BFGS (provided by the NLopt library).

With multiple runs of optimization (by using the restarts argument), it is possible to parallelize these runs using the future package, e.g., by calling future::plan(multisession, workers = 2) before estimate_nhmm(). See future::plan() for details. This is compatible with progressr package, so you can use progressr::with_progress() to track the progress of these multiple runs.

During the estimation, the log-likelihood is scaled by the number of non-missing observations (nobs(model)), and the the covariate data is standardardized before optimization.

By default, the convergence is claimed when the relative change of the objective function is less than 1e-12, the absolute change is less than 1e-8 or the relative or absolute change of the working parameters eta is less than 1e-6. These can be changed by passing arguments ftol_rel, ftol_abs, xtol_rel, and xtol_abs via .... These, as well as, maxeval (maximum number of iterations, 1e4 by default), and print_level (default is 0, no console output, larger values are more verbose), are used by the chosen main optimization method. The number of initial EM iterations in EM-DNM can be set using argument maxeval_em_dnm (default is 100), and algorithm for direct numerical optimization can be defined using argument algorithm (see nloptr::nloptr() for possible options).

For controlling these stopping criteria for the multistart phase, argument control_restart takes a list such as list(ftol_rel = 0.01, print_level = 1). Default are as in the case of main optimization (which is always run once after the restarts, using best solution from restarts as initial value) Additionally, same options can be defined separately for the M-step of EM algorithm via list control_mstep. For control_mstep, the default values are ftol_rel = 1e-10, and maxeval = 1000, and otherwise identical to previous defaults above.

Value

Object of class nhmm or fanhmm.

References

Helske, J (2025). Feedback-augmented Non-homogeneous Hidden Markov Models for Longitudinal Causal Inference. arXiv preprint. doi:10.48550/arXiv.2503.16014.

Johnson, SG. The NLopt nonlinear-optimization package, http://github.com/stevengj/nlopt.

Examples

data("mvad", package = "TraMineR")

d <- reshape(mvad, direction = "long", varying = list(15:86), 
  v.names = "activity")

## Not run: 
set.seed(1)
fit <- estimate_nhmm(n_states = 3,
  data = d, time = "time", id = "id", 
  emission_formula = activity ~ gcse5eq, initial_formula = ~ 1, 
  transition_formula = ~ male + gcse5eq,
  method = "DNM", maxeval = 2 # very small number of iterations for CRAN
  )

## End(Not run)

A feedback-augmented non-homogeneuous hidden Markov Model for leaves data

Description

A FAN-HMM fitted for theleaes data.

Format

A model of class fanhmm with three hidden states

Details

The model is loaded by calling data(fanhmm_leaves). The code used to estimate the model is available on Github in data-raw folder.

Examples

data("fanhmm_leaves")

fanhmm_leaves

get_marginals(fanhmm_leaves)

Estimate Parameters of (Mixture) Hidden Markov Models and Their Restricted Variants

Description

Function fit_model estimates the parameters of mixture hidden Markov models and its restricted variants using maximimum likelihood. Initial values for estimation are taken from the corresponding components of the model with preservation of original zero probabilities.

Usage

fit_model(
  model,
  em_step = TRUE,
  global_step = FALSE,
  local_step = FALSE,
  control_em = list(),
  control_global = list(),
  control_local = list(),
  lb,
  ub,
  threads = 1,
  log_space = TRUE,
  constraints = NULL,
  fixed_inits = NULL,
  fixed_emissions = NULL,
  fixed_transitions = NULL,
  ...
)

Arguments

Details

The fitting function provides three estimation steps: 1) EM algorithm, 2) global optimization, and 3) local optimization. The user can call for one method or any combination of these steps, but should note that they are preformed in the above-mentioned order. The results from a former step are used as starting values in a latter, except for some of global optimization algorithms (such as MLSL and StoGO) which only use initial values for setting up the boundaries for the optimization.

It is possible to rerun the EM algorithm automatically using random starting values based on the first run of EM. Number of restarts is defined by the restart argument in control_em. As the EM algorithm is relatively fast, this method might be preferred option compared to the proper global optimization strategy of step 2.

The default global optimization method (triggered via global_step = TRUE) is the multilevel single-linkage method (MLSL) with the LDS modification (NLOPT_GD_MLSL_LDS as algorithm in control_global), with L-BFGS as the local optimizer. The MLSL method draws random starting points and performs a local optimization from each. The LDS modification uses low-discrepancy sequences instead of pseudo-random numbers as starting points and should improve the convergence rate. In order to reduce the computation time spent on non-global optima, the convergence tolerance of the local optimizer is set relatively large. At step 3, a local optimization (L-BFGS by default) is run with a lower tolerance to find the optimum with high precision.

There are some theoretical guarantees that the MLSL method used as the default optimizer in step 2 shoud find all local optima in a finite number of local optimizations. Of course, it might not always succeed in a reasonable time. The EM algorithm can help in finding good boundaries for the search, especially with good starting values, but in some cases it can mislead. A good strategy is to try a couple of different fitting options with different combinations of the methods: e.g. all steps, only global and local steps, and a few evaluations of EM followed by global and local optimization.

By default, the estimation time is limited to 60 seconds in global optimization step, so it is advisable to change the default settings for the proper global optimization.

Any algorithm available in the nloptr function can be used for the global and local steps.

Equality constraints for emission distributions can be defined using the argument constraints. This should be a vector with length equal to the number of states, with numbers starting from 1 and increasing for each unique row of the emission probability matrix. For example in case of five states with emissions of first and third states being equal, constraints = c(1, 2, 1, 3, 4). Similarly, some of the model parameters can be fixed to their initial values by using arguments fixed_inits, fixed_emissions, and fixed_transitions, where the structure of the arguments should be same as the corresponding model components, so that TRUE value means that the parameter should be fixed and FALSE otherwise (it is still treated as fixed if it is zero though). For both types of constrains, only numerical optimisation (local or global) is available, and currently the gradients are computed numerically (if needed) in these cases.

In a case where the is no transitions from one state to anywhere (even to itself), the state is defined as absorbing in a way that probability of staying in this state is fixed to 1. See also build_mm function.

Value

• logLik
  Log-likelihood of the estimated model.

• em_results
  Results after the EM step: log-likelihood (logLik), number of iterations (iterations), relative change in log-likelihoods between the last two iterations (change), and the log-likelihoods of the n_optimum best models after the EM step (best_opt_restart).

• global_results
  Results after the global step.

• local_results
  Results after the local step.

• call
  The matched function call.

References

Helske S. and Helske J. (2019). Mixture Hidden Markov Models for Sequence Data: The seqHMM Package in R, Journal of Statistical Software, 88(3), 1-32. doi:10.18637/jss.v088.i03

See Also

build_hmm(), build_mhmm(), build_mm(), build_mmm(), and build_lcm() for constructing different types of models; summary.mhmm() for a summary of a MHMM; separate_mhmm() for reorganizing a MHMM into a list of separate hidden Markov models; and plot.hmm() and plot.mhmm() for plotting model objects.

Examples

# Hidden Markov model for mvad data

data("mvad", package = "TraMineR")

mvad_alphabet <-
  c("employment", "FE", "HE", "joblessness", "school", "training")
mvad_labels <- c(
  "employment", "further education", "higher education",
  "joblessness", "school", "training"
)
mvad_scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad_seq <- seqdef(mvad, 15:86,
  alphabet = mvad_alphabet,
  states = mvad_scodes, labels = mvad_labels, xtstep = 6,
  cpal = colorpalette[[6]]
)

# Starting values for the emission matrix
emiss <- matrix(
  c(
    0.05, 0.05, 0.05, 0.05, 0.75, 0.05, # SC
    0.05, 0.75, 0.05, 0.05, 0.05, 0.05, # FE
    0.05, 0.05, 0.05, 0.4, 0.05, 0.4, # JL, TR
    0.05, 0.05, 0.75, 0.05, 0.05, 0.05, # HE
    0.75, 0.05, 0.05, 0.05, 0.05, 0.05
  ), # EM
  nrow = 5, ncol = 6, byrow = TRUE
)

# Starting values for the transition matrix
trans <- matrix(0.025, 5, 5)
diag(trans) <- 0.9

# Starting values for initial state probabilities
initial_probs <- c(0.2, 0.2, 0.2, 0.2, 0.2)

# Building a hidden Markov model
init_hmm_mvad <- build_hmm(
  observations = mvad_seq,
  transition_probs = trans, emission_probs = emiss,
  initial_probs = initial_probs
)

## Not run: 
set.seed(21)
fit_hmm_mvad <- fit_model(init_hmm_mvad, control_em = list(restart = list(times = 50)))
hmm_mvad <- fit_hmm_mvad$model

## End(Not run)

# save time, load the previously estimated model
data("hmm_mvad")

# Markov model
# Note: build_mm estimates model parameters from observations,
# no need for estimating with fit_model unless there are missing observations

mm_mvad <- build_mm(observations = mvad_seq)

# Comparing likelihoods, MM fits better
logLik(hmm_mvad)
logLik(mm_mvad)

## Not run: 
require("igraph") # for layout_in_circle

plot(mm_mvad,
  layout = layout_in_circle, legend.prop = 0.3,
  edge.curved = 0.3, edge.label = NA,
  vertex.label.pos = c(0, 0, pi, pi, pi, 0)
)

##############################################################


#' # Three-state three-channel hidden Markov model
# See ?hmm_biofam for five-state version

data("biofam3c")

# Building sequence objects
marr_seq <- seqdef(biofam3c$married,
  start = 15,
  alphabet = c("single", "married", "divorced"),
  cpal = c("violetred2", "darkgoldenrod2", "darkmagenta")
)
child_seq <- seqdef(biofam3c$children,
  start = 15,
  alphabet = c("childless", "children"),
  cpal = c("darkseagreen1", "coral3")
)
left_seq <- seqdef(biofam3c$left,
  start = 15,
  alphabet = c("with parents", "left home"),
  cpal = c("lightblue", "red3")
)

# Starting values for emission matrices

emiss_marr <- matrix(NA, nrow = 3, ncol = 3)
emiss_marr[1, ] <- seqstatf(marr_seq[, 1:5])[, 2] + 1
emiss_marr[2, ] <- seqstatf(marr_seq[, 6:10])[, 2] + 1
emiss_marr[3, ] <- seqstatf(marr_seq[, 11:16])[, 2] + 1
emiss_marr <- emiss_marr / rowSums(emiss_marr)

emiss_child <- matrix(NA, nrow = 3, ncol = 2)
emiss_child[1, ] <- seqstatf(child_seq[, 1:5])[, 2] + 1
emiss_child[2, ] <- seqstatf(child_seq[, 6:10])[, 2] + 1
emiss_child[3, ] <- seqstatf(child_seq[, 11:16])[, 2] + 1
emiss_child <- emiss_child / rowSums(emiss_child)

emiss_left <- matrix(NA, nrow = 3, ncol = 2)
emiss_left[1, ] <- seqstatf(left_seq[, 1:5])[, 2] + 1
emiss_left[2, ] <- seqstatf(left_seq[, 6:10])[, 2] + 1
emiss_left[3, ] <- seqstatf(left_seq[, 11:16])[, 2] + 1
emiss_left <- emiss_left / rowSums(emiss_left)

# Starting values for transition matrix
trans <- matrix(c(
  0.9, 0.07, 0.03,
  0, 0.9, 0.1,
  0, 0, 1
), nrow = 3, ncol = 3, byrow = TRUE)

# Starting values for initial state probabilities
inits <- c(0.9, 0.09, 0.01)

# Building hidden Markov model with initial parameter values
init_hmm_bf <- build_hmm(
  observations = list(marr_seq, child_seq, left_seq),
  transition_probs = trans,
  emission_probs = list(emiss_marr, emiss_child, emiss_left),
  initial_probs = inits
)

# Fitting the model with different optimization schemes

# Only EM with default values
hmm_1 <- fit_model(init_hmm_bf)
hmm_1$logLik # -24179.1

# Only L-BFGS
hmm_2 <- fit_model(init_hmm_bf, em_step = FALSE, local_step = TRUE)
hmm_2$logLik # -22267.75

# Global optimization via MLSL_LDS with L-BFGS as local optimizer and final polisher
# This can be slow, use parallel computing by adjusting threads argument
# (here threads = 1 for portability issues)
hmm_3 <- fit_model(
  init_hmm_bf,
  em_step = FALSE, global_step = TRUE, local_step = TRUE,
  control_global = list(maxeval = 5000, maxtime = 0), threads = 1
)
hmm_3$logLik # -21675.42

# EM with restarts, much faster than MLSL
set.seed(123)
hmm_4 <- fit_model(init_hmm_bf, control_em = list(restart = list(times = 5)))
hmm_4$logLik # -21675.4

# Global optimization via StoGO with L-BFGS as final polisher
# This can be slow, use parallel computing by adjusting threads argument
# (here threads = 1 for portability issues)
set.seed(123)
hmm_5 <- fit_model(
  init_hmm_bf,
  em_step = FALSE, global_step = TRUE, local_step = TRUE,
  lb = -50, ub = 50, control_global = list(
    algorithm = "NLOPT_GD_STOGO",
    maxeval = 2500, maxtime = 0
  ), threads = 1
)
hmm_5$logLik # -21675.4

##############################################################

# Mixture HMM

data("biofam3c")

# Building sequence objects
marr_seq <- seqdef(biofam3c$married,
  start = 15,
  alphabet = c("single", "married", "divorced"),
  cpal = c("violetred2", "darkgoldenrod2", "darkmagenta")
)
child_seq <- seqdef(biofam3c$children,
  start = 15,
  alphabet = c("childless", "children"),
  cpal = c("darkseagreen1", "coral3")
)
left_seq <- seqdef(biofam3c$left,
  start = 15,
  alphabet = c("with parents", "left home"),
  cpal = c("lightblue", "red3")
)

## Starting values for emission probabilities
# Cluster 1
B1_marr <- matrix(
  c(
    0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.3, 0.6, 0.1, # High probability for married
    0.3, 0.3, 0.4
  ), # High probability for divorced
  nrow = 4, ncol = 3, byrow = TRUE
)

B1_child <- matrix(
  c(
    0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.9, 0.1,
    0.9, 0.1
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

B1_left <- matrix(
  c(
    0.9, 0.1, # High probability for living with parents
    0.1, 0.9, # High probability for having left home
    0.1, 0.9,
    0.1, 0.9
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

# Cluster 2

B2_marr <- matrix(
  c(
    0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.1, 0.8, 0.1, # High probability for married
    0.7, 0.2, 0.1
  ),
  nrow = 4, ncol = 3, byrow = TRUE
)

B2_child <- matrix(
  c(
    0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.9, 0.1,
    0.1, 0.9
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

B2_left <- matrix(
  c(
    0.9, 0.1, # High probability for living with parents
    0.1, 0.9,
    0.1, 0.9,
    0.1, 0.9
  ),
  nrow = 4, ncol = 2, byrow = TRUE
)

# Cluster 3
B3_marr <- matrix(
  c(
    0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.8, 0.1, 0.1,
    0.1, 0.8, 0.1, # High probability for married
    0.3, 0.4, 0.3,
    0.1, 0.1, 0.8
  ), # High probability for divorced
  nrow = 6, ncol = 3, byrow = TRUE
)

B3_child <- matrix(
  c(
    0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.5, 0.5,
    0.5, 0.5,
    0.5, 0.5,
    0.1, 0.9
  ),
  nrow = 6, ncol = 2, byrow = TRUE
)


B3_left <- matrix(
  c(
    0.9, 0.1, # High probability for living with parents
    0.1, 0.9,
    0.5, 0.5,
    0.5, 0.5,
    0.1, 0.9,
    0.1, 0.9
  ),
  nrow = 6, ncol = 2, byrow = TRUE
)

# Starting values for transition matrices
A1 <- matrix(
  c(
    0.80, 0.16, 0.03, 0.01,
    0, 0.90, 0.07, 0.03,
    0, 0, 0.90, 0.10,
    0, 0, 0, 1
  ),
  nrow = 4, ncol = 4, byrow = TRUE
)

A2 <- matrix(
  c(
    0.80, 0.10, 0.05, 0.03, 0.01, 0.01,
    0, 0.70, 0.10, 0.10, 0.05, 0.05,
    0, 0, 0.85, 0.01, 0.10, 0.04,
    0, 0, 0, 0.90, 0.05, 0.05,
    0, 0, 0, 0, 0.90, 0.10,
    0, 0, 0, 0, 0, 1
  ),
  nrow = 6, ncol = 6, byrow = TRUE
)

# Starting values for initial state probabilities
initial_probs1 <- c(0.9, 0.07, 0.02, 0.01)
initial_probs2 <- c(0.9, 0.04, 0.03, 0.01, 0.01, 0.01)

# Birth cohort
biofam3c$covariates$cohort <- cut(biofam3c$covariates$birthyr, c(1908, 1935, 1945, 1957))
biofam3c$covariates$cohort <- factor(
  biofam3c$covariates$cohort,
  labels = c("1909-1935", "1936-1945", "1946-1957")
)

# Build mixture HMM
init_mhmm_bf <- build_mhmm(
  observations = list(marr_seq, child_seq, left_seq),
  initial_probs = list(initial_probs1, initial_probs1, initial_probs2),
  transition_probs = list(A1, A1, A2),
  emission_probs = list(
    list(B1_marr, B1_child, B1_left),
    list(B2_marr, B2_child, B2_left),
    list(B3_marr, B3_child, B3_left)
  ),
  formula = ~ sex + cohort, data = biofam3c$covariates,
  channel_names = c("Marriage", "Parenthood", "Residence")
)


# Fitting the model with different settings

# Only EM with default values
mhmm_1 <- fit_model(init_mhmm_bf)
mhmm_1$logLik # -12713.08

# Only L-BFGS
mhmm_2 <- fit_model(init_mhmm_bf, em_step = FALSE, local_step = TRUE)
mhmm_2$logLik # -12966.51

# Use EM with multiple restarts
set.seed(123)
mhmm_3 <- fit_model(init_mhmm_bf, control_em = list(restart = list(times = 5, transition = FALSE)))
mhmm_3$logLik # -12713.08

## End(Not run)

# Left-to-right HMM with equality constraint:

set.seed(1)

# Transition matrix
# Either stay or move to next state
A <- diag(c(0.9, 0.95, 0.95, 1))
A[1, 2] <- 0.1
A[2, 3] <- 0.05
A[3, 4] <- 0.05

# Emission matrix, rows 1 and 3 equal
B <- rbind(
  c(0.4, 0.2, 0.3, 0.1),
  c(0.1, 0.5, 0.1, 0.3),
  c(0.4, 0.2, 0.3, 0.1),
  c(0, 0.2, 0.4, 0.4)
)

# Start from first state
init <- c(1, 0, 0, 0)

# Simulate sequences
sim <- simulate_hmm(
  n_sequences = 100,
  sequence_length = 20, init, A, B
)

# initial model, use true values as inits for faster estimation here
model <- build_hmm(sim$observations, init = init, trans = A, emiss = B)

# estimate the model subject to constraints:
# First and third row of emission matrix are equal (see details)
fit <- fit_model(model,
  constraints = c(1, 2, 1, 3),
  em_step = FALSE, local_step = TRUE
)
fit$model

## Fix some emissions:

fixB <- matrix(FALSE, 4, 4)
fixB[2, 1] <- fixB[1, 3] <- TRUE # these are fixed to their initial values
fit <- fit_model(model,
  fixed_emissions = fixB,
  em_step = FALSE, local_step = TRUE
)
fit$model$emission_probs

Forward and Backward Probabilities for Hidden Markov Model

Description

The forward_backward function computes forward and backward probabilities of a hidden Markov model.

Usage

forward_backward(model, ...)

## S3 method for class 'hmm'
forward_backward(model, forward_only = FALSE, ...)

## S3 method for class 'mhmm'
forward_backward(model, forward_only = FALSE, ...)

## S3 method for class 'nhmm'
forward_backward(model, forward_only = FALSE, ...)

## S3 method for class 'mnhmm'
forward_backward(model, forward_only = FALSE, ...)

Arguments

Value

A data.frame with log-values of forward and backward probabilities.

Examples

# Load a pre-defined MHMM
data("mhmm_biofam")

# Compute forward and backward probabilities
fb <- forward_backward(mhmm_biofam)

head(fb)

Extract the Prior Cluster Probabilities of MHMM or MNHMM

Description

Extract the Prior Cluster Probabilities of MHMM or MNHMM

Usage

get_cluster_probs(model)

## S3 method for class 'mnhmm'
get_cluster_probs(model)

## S3 method for class 'mhmm'
get_cluster_probs(model)

Arguments

See Also

posterior_cluster_probabilities().

Extract the Emission Probabilities of Hidden Markov Model

Description

Extract the Emission Probabilities of Hidden Markov Model

Usage

get_emission_probs(model)

## S3 method for class 'nhmm'
get_emission_probs(model)

## S3 method for class 'mnhmm'
get_emission_probs(model)

## S3 method for class 'hmm'
get_emission_probs(model)

## S3 method for class 'mhmm'
get_emission_probs(model)

Arguments

Extract the Initial State Probabilities of Hidden Markov Model

Description

Extract the Initial State Probabilities of Hidden Markov Model

Usage

get_initial_probs(model)

## S3 method for class 'nhmm'
get_initial_probs(model)

## S3 method for class 'mnhmm'
get_initial_probs(model)

## S3 method for class 'hmm'
get_initial_probs(model)

## S3 method for class 'mhmm'
get_initial_probs(model)

Arguments

Compute the Marginal Probabilities from NHMMs

Description

get_marginals returns the marginal state, response, transition, and emission probabilities, optionally per grouping defined by condition. By default, the marginalization weights sequences by the corresponding posterior probabilities of the latent states, i.e., conditional probabilities of the latent states given all data (weighting = "posterior"). If weighting = "forward", marginalization is based on forward probabilities, i.e. state probabilities given data up to that point which allows you to compute, for example, state marginals of form P(state_t | data_1, \ldots, data_t) (whereas in posterior probability weighting the conditioning is on data_1,\ldots,data_T. If weighting = "none", all individuals and time points are treated equally, without accounting for the probability that individual is at particular state at particular time.

Usage

get_marginals(
  model,
  probs = NULL,
  condition = NULL,
  newdata = NULL,
  type = c("state", "response", "transition", "emission"),
  weighting = c("posterior", "forward", "none")
)

Arguments

Extract the State Transition Probabilities of Hidden Markov Model

Description

Extract the State Transition Probabilities of Hidden Markov Model

Usage

get_transition_probs(model)

## S3 method for class 'nhmm'
get_transition_probs(model)

## S3 method for class 'mnhmm'
get_transition_probs(model)

## S3 method for class 'hmm'
get_transition_probs(model)

## S3 method for class 'mhmm'
get_transition_probs(model)

Arguments

Plot Multidimensional Sequence Plots in a Grid

Description

Function gridplot plots multiple ssp objects to a grid.

Usage

gridplot(
  x,
  nrow = NA,
  ncol = NA,
  byrow = FALSE,
  with.legend = "auto",
  legend.pos = "auto",
  legend.pos2 = "center",
  title.legend = "auto",
  ncol.legend = "auto",
  with.missing.legend = "auto",
  row.prop = "auto",
  col.prop = "auto",
  cex.legend = 1
)

Arguments

Most Probable Paths of Hidden States

Description

Function hidden_paths computes the most probable path of hidden states of a (mixture) hidden Markov model given the observed sequences.

Usage

hidden_paths(model, ...)

## S3 method for class 'hmm'
hidden_paths(model, as_stslist = FALSE, ...)

## S3 method for class 'mhmm'
hidden_paths(model, as_stslist = FALSE, ...)

## S3 method for class 'nhmm'
hidden_paths(model, as_stslist = FALSE, ...)

## S3 method for class 'mnhmm'
hidden_paths(model, as_stslist = FALSE, ...)

Arguments

Value

The most probable paths of hidden states as an data.table. The log-probability is included as an attribute log_prop.

See Also

hmm_biofam for information on the model used in the example; and ggseqplot::ggseqiplot() and stacked_sequence_plot() for plotting hidden paths.

Examples

# Load a pre-defined HMM
data("hmm_biofam")

# Compute the most probable hidden state paths given the data and the model
mpp <- hidden_paths(hmm_biofam)
head(mpp)
# Plot hidden paths for the first 100 individuals
seqs <- data_to_stslist(mpp, "id", "time", "state")
stacked_sequence_plot(seqs, type = "i", ids = 1:100)

# Because the model structure is so sparse that the posterior probabilities are
# mostly peaked to single state at each time point, the joint probability of
# observations and most probable paths of hidden states is almost identical to
# log-likelihood:

sum(attr(mpp, "log_prob"))
logLik(hmm_biofam)

Hidden Markov model for the biofam data

Description

A five-state hidden Markov model (HMM) fitted for the TraMineR::biofam() data.

Format

A hidden Markov model of class hmm; a left-to-right model with four hidden states.

Details

The model is loaded by calling data(hmm_biofam). It was created with the following code:

data("biofam3c")

# Building sequence objects
marr_seq <- seqdef(biofam3c$married,
  start = 15,
  alphabet = c("single", "married", "divorced"),
  cpal = c("violetred2", "darkgoldenrod2", "darkmagenta")
)
child_seq <- seqdef(biofam3c$children,
  start = 15,
  alphabet = c("childless", "children"),
  cpal = c("darkseagreen1", "coral3")
)
left_seq <- seqdef(biofam3c$left,
  start = 15,
  alphabet = c("with parents", "left home"),
  cpal = c("lightblue", "red3")
)

init <- c(0.9, 0.05, 0.02, 0.02, 0.01)

# Starting values for transition matrix
trans <- matrix(
  c(0.8, 0.10, 0.05, 0.03, 0.02,
    0,    0.9, 0.05, 0.03, 0.02,
    0,      0,  0.9, 0.07, 0.03,
    0,      0,    0,  0.9,  0.1,
    0,      0,    0,    0,    1),
  nrow = 5, ncol = 5, byrow = TRUE)

# Starting values for emission matrices
emiss_marr <- matrix(
  c(0.9, 0.05, 0.05, # High probability for single
    0.9, 0.05, 0.05,
    0.05, 0.9, 0.05, # High probability for married
    0.05, 0.9, 0.05,
    0.3, 0.3, 0.4), # mixed group
  nrow = 5, ncol = 3, byrow = TRUE)

emiss_child <- matrix(
  c(0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.1, 0.9,
    0.1, 0.9,
    0.5, 0.5),
  nrow = 5, ncol = 2, byrow = TRUE)

emiss_left <- matrix(
  c(0.9, 0.1, # High probability for living with parents
    0.1, 0.9,
    0.1, 0.9,
    0.1, 0.9,
    0.5, 0.5),
  nrow = 5, ncol = 2, byrow = TRUE)

initmod <- build_hmm(
  observations = list(marr_seq, child_seq, left_seq),
  initial_probs = init, transition_probs = trans,
  emission_probs = list(emiss_marr, emiss_child,
    emiss_left),
  channel_names = c("Marriage", "Parenthood", "Residence"))

fit_biofam <- fit_model(initmod, em = FALSE, local = TRUE)
hmm_biofam <- fit_biofam$model

See Also

Examples of building and fitting HMMs in build_hmm() and fit_model(); and TraMineR::biofam() for the original data and biofam3c() for the three-channel version used in this model.

Examples

# Plotting the model
plot(hmm_biofam)

Hidden Markov model for the mvad data

Description

A hidden Markov model (MMM) fitted for the TraMineR::mvad() data.

Format

A hidden Markov model of class hmm; unrestricted model with six hidden states.

Details

Model was created with the following code:

data("mvad", package = "TraMineR")

mvad_alphabet <-
  c("employment", "FE", "HE", "joblessness", "school", "training")
mvad_labels <- c("employment", "further education", "higher education",
  "joblessness", "school", "training")
mvad_scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad_seq <- seqdef(mvad, 15:86, alphabet = mvad_alphabet,
  states = mvad_scodes, labels = mvad_labels, xtstep = 6,
  cpal = colorpalette[[6]])

# Starting values for the emission matrix
emiss <- matrix(
  c(0.05, 0.05, 0.05, 0.05, 0.75, 0.05, # SC
    0.05, 0.75, 0.05, 0.05, 0.05, 0.05, # FE
    0.05, 0.05, 0.05, 0.4,  0.05, 0.4,  # JL, TR
    0.05, 0.05, 0.75, 0.05, 0.05, 0.05, # HE
    0.75, 0.05, 0.05, 0.05, 0.05, 0.05),# EM
  nrow = 5, ncol = 6, byrow = TRUE)

# Starting values for the transition matrix
trans <- matrix(0.025, 5, 5)
diag(trans) <- 0.9

# Starting values for initial state probabilities
initial_probs <- c(0.2, 0.2, 0.2, 0.2, 0.2)

# Building a hidden Markov model
init_hmm_mvad <- build_hmm(observations = mvad_seq,
  transition_probs = trans, emission_probs = emiss,
  initial_probs = initial_probs)

set.seed(21)
fit_hmm_mvad <- fit_model(init_hmm_mvad, control_em = list(restart = list(times = 100)))
hmm_mvad <- fit_hmm_mvad$model

See Also

Examples of building and fitting HMMs in build_hmm() and fit_model(); and TraMineR::mvad() for more information on the data.

Examples

data("hmm_mvad")

# Plotting the model
plot(hmm_mvad)

Synthetic data on fathers' parental leaves in Finland

Description

Synthetic data on fathers' parental leaves in Finland

Format

A data.table with 9281 rows and 9 variables

Details

The leaves data is a synthetic version of the Finnish fathers' leave-taking data used in Helske et al. (2024) and Helske (2025). The data consists of variables

• workplace: Workplace ID.

• father: father ID within workplace. More accurately, this is the birth of a child, i.e. same father can have multiple entries in data, but each entry has separate ID.

• year: Year when the child was born.

• leave: Factor of leave-taking of the father.

• Occupation: Factor of skill level of the father's occupation

• reform2013: Factor indicating whether the father was eligible for the leave under the 2013 reform.

• same_occupation: Logical value, TRUE if father had same occupation as the previous father.

• lag_reform2013: Factor indicating whether the previous father was eligible for the reform.

• lag_occupation: Factor indiciting the occupation of previous father.

References

Helske S, Helske J, Chapman SN, Kotimäki S, Salin M, and Tikka S (2024). Heterogeneous workplace peer effects in fathers’ parental leave uptake in Finland. doi: 10.31235/osf.io/p3chf Helske J (2025). Feedback-augmented Non-homogeneous Hidden Markov Models for Longitudinal Causal Inference. ArXiv preprint. doi:10.48550/arXiv.2503.16014

Examples

data("leaves")
head(leaves)
# convert to stslist
leaves_sequences <- data_to_stslist(
  leaves, id = "workplace", time = "father", responses = "leave",
  seqdef_args = list(cpal = c("tomato", "navyblue", "goldenrod"))
)
stacked_sequence_plot(leaves_sequences)

Log-likelihood of a Hidden Markov Model

Description

Log-likelihood of a Hidden Markov Model

Usage

## S3 method for class 'hmm'
logLik(object, partials = FALSE, threads = 1, log_space = TRUE, ...)

## S3 method for class 'mhmm'
logLik(object, partials = FALSE, threads = 1, log_space = TRUE, ...)

Arguments

Value

Log-likelihood of the hidden Markov model. This is an object of class logLik with attributes nobs and df inherited from the model object.

Log-likelihood of a Non-homogeneous Hidden Markov Model

Description

Log-likelihood of a Non-homogeneous Hidden Markov Model

Usage

## S3 method for class 'nhmm'
logLik(object, partials = FALSE, ...)

## S3 method for class 'mnhmm'
logLik(object, partials = FALSE, ...)

Arguments

Value

Log-likelihood of the hidden Markov model. This is an object of class logLik with attributes nobs and df inherited from the model object.

Transform a Multichannel Hidden Markov Model into a Single Channel Representation

Description

Transforms data and parameters of a multichannel model into a single channel model. Observed states (symbols) are combined and parameters multiplied across channels.

Usage

mc_to_sc(model, combine_missing = TRUE, all_combinations = FALSE, cpal)

Arguments

Details

Note that in case of no missing observations, the log-likelihood of the original and transformed models are identical but the AIC and BIC can be different as the model attribute df is recomputed based on the single channel representation.

See Also

build_hmm() and fit_model() for building and fitting Hidden Markov models; and hmm_biofam() for information on the model used in the example.

Examples

# Loading a hidden Markov model of the biofam data (hmm object)
data("hmm_biofam")

# Convert the multichannel model to a single-channel model
sc <- mc_to_sc(hmm_biofam)

# Likelihoods of the single-channel and the multichannel model are the same
# (Might not be true if there are missing observations)
logLik(sc)
logLik(hmm_biofam)

Merge Multiple Sequence Objects into One (from Multichannel to Single Channel Data)

Description

Function mc_to_sc_data combines observed states of multiple sequence objects into one, time point by time point.

Usage

mc_to_sc_data(data, combine_missing = TRUE, all_combinations = FALSE, cpal)

Arguments

See Also

mc_to_sc() for transforming multichannel hmm or mhmm objects into single-channel representations; stacked_sequence_plot for plotting multiple sequence data sets in the same plot; and seqdef() for creating state sequence objects.

Examples

# Load three-channel sequence data
data("biofam3c")

# Building sequence objects
marr_seq <- seqdef(biofam3c$married,
  start = 15,
  alphabet = c("single", "married", "divorced"),
  cpal = c("violetred2", "darkgoldenrod2", "darkmagenta")
)
child_seq <- seqdef(biofam3c$children,
  start = 15,
  alphabet = c("childless", "children"),
  cpal = c("darkseagreen1", "coral3")
)
left_seq <- seqdef(biofam3c$left,
  start = 15,
  alphabet = c("with parents", "left home"),
  cpal = c("lightblue", "red3")
)

# Converting multichannel data to single-channel data
sc_data <- mc_to_sc_data(list(marr_seq, child_seq, left_seq))

# 10 combined states
alphabet(sc_data)

# Colors for combined states
attr(sc_data, "cpal") <- colorpalette[[14]][1:10]

# Plotting sequences for the first 10 subjects
stacked_sequence_plot(
  list(
    "Marriage" = marr_seq, "Parenthood" = child_seq,
    "Residence" = left_seq, "Combined" = sc_data
  ),
  type = "i",
  ids = 1:10
)


# Including all combinations (whether or not available in data)
sc_data_all <- mc_to_sc_data(list(marr_seq, child_seq, left_seq),
  all_combinations = TRUE
)

# 12 combined states, 2 with no observations in data
seqstatf(sc_data_all)

Mixture hidden Markov model for the biofam data

Description

A mixture hidden Markov model (MHMM) fitted for the TraMineR::biofam() data.

Format

A mixture hidden Markov model of class mhmm: three clusters with left-to-right models including 4, 4, and 6 hidden states. Two covariates, sex and cohort, explaining the cluster membership.

Details

The model was created with the following code:

data("biofam3c")

# Building sequence objects
marr_seq <- seqdef(biofam3c$married,
  start = 15,
  alphabet = c("single", "married", "divorced"),
  cpal = c("violetred2", "darkgoldenrod2", "darkmagenta")
)
child_seq <- seqdef(biofam3c$children,
  start = 15,
  alphabet = c("childless", "children"),
  cpal = c("darkseagreen1", "coral3")
)
left_seq <- seqdef(biofam3c$left,
  start = 15,
  alphabet = c("with parents", "left home"),
  cpal = c("lightblue", "red3")
)

## Starting values for emission probabilities
# Cluster 1
B1_marr <- matrix(
  c(0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.3, 0.6, 0.1, # High probability for married
    0.3, 0.3, 0.4), # High probability for divorced
  nrow = 4, ncol = 3, byrow = TRUE)

B1_child <- matrix(
  c(0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.9, 0.1,
    0.9, 0.1),
  nrow = 4, ncol = 2, byrow = TRUE)

B1_left <- matrix(
  c(0.9, 0.1, # High probability for living with parents
    0.1, 0.9, # High probability for having left home
    0.1, 0.9,
    0.1, 0.9),
  nrow = 4, ncol = 2, byrow = TRUE)

# Cluster 2

B2_marr <- matrix(
  c(0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.1, 0.8, 0.1, # High probability for married
    0.7, 0.2, 0.1),
  nrow = 4, ncol = 3, byrow = TRUE)

B2_child <- matrix(
  c(0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.9, 0.1,
    0.1, 0.9),
  nrow = 4, ncol = 2, byrow = TRUE)

B2_left <- matrix(
  c(0.9, 0.1, # High probability for living with parents
    0.1, 0.9,
    0.1, 0.9,
    0.1, 0.9),
  nrow = 4, ncol = 2, byrow = TRUE)

# Cluster 3
B3_marr <- matrix(
  c(0.8, 0.1, 0.1, # High probability for single
    0.8, 0.1, 0.1,
    0.8, 0.1, 0.1,
    0.1, 0.8, 0.1, # High probability for married
    0.3, 0.4, 0.3,
    0.1, 0.1, 0.8), # High probability for divorced
  nrow = 6, ncol = 3, byrow = TRUE)

B3_child <- matrix(
  c(0.9, 0.1, # High probability for childless
    0.9, 0.1,
    0.5, 0.5,
    0.5, 0.5,
    0.5, 0.5,
    0.1, 0.9),
  nrow = 6, ncol = 2, byrow = TRUE)


B3_left <- matrix(
  c(0.9, 0.1, # High probability for living with parents
    0.1, 0.9,
    0.5, 0.5,
    0.5, 0.5,
    0.1, 0.9,
    0.1, 0.9),
  nrow = 6, ncol = 2, byrow = TRUE)

# Starting values for transition matrices
A1 <- matrix(
  c(0.80, 0.16, 0.03, 0.01,
    0,    0.90, 0.07, 0.03,
    0,    0,    0.90, 0.10,
    0,    0,    0,       1),
  nrow = 4, ncol = 4, byrow = TRUE)

A2 <- matrix(
  c(0.80, 0.10, 0.05, 0.03, 0.01, 0.01,
    0,    0.70, 0.10, 0.10, 0.05, 0.05,
    0,    0,    0.85, 0.01, 0.10, 0.04,
    0,    0,    0,    0.90, 0.05, 0.05,
    0,    0,    0,    0,    0.90, 0.10,
    0,    0,    0,    0,    0,       1),
  nrow = 6, ncol = 6, byrow = TRUE)

# Starting values for initial state probabilities
initial_probs1 <- c(0.9, 0.07, 0.02, 0.01)
initial_probs2 <- c(0.9, 0.04, 0.03, 0.01, 0.01, 0.01)

# Birth cohort
biofam3c$covariates$cohort <- factor(cut(biofam3c$covariates$birthyr,
    c(1908, 1935, 1945, 1957)), labels = c("1909-1935", "1936-1945", "1946-1957"))

# Build mixture HMM
init_mhmm_bf <- build_mhmm(
  observations = list(marr_seq, child_seq, left_seq),
  initial_probs = list(initial_probs1, initial_probs1, initial_probs2),
  transition_probs = list(A1, A1, A2),
  emission_probs = list(list(B1_marr, B1_child, B1_left),
    list(B2_marr, B2_child, B2_left),
    list(B3_marr, B3_child, B3_left)),
  formula = ~sex + cohort, data = biofam3c$covariates,
  channel_names = c("Marriage", "Parenthood", "Residence"))

# Fitting the model
mhmm_biofam <- fit_model(init_mhmm_bf)$model

See Also

Examples of building and fitting MHMMs in build_mhmm() and fit_model(); and TraMineR::biofam() for the original data and biofam3c() for the three-channel version used in this model.

Examples

data("mhmm_biofam")

# use conditional_se = FALSE for more accurate standard errors
# (these are considerebly slower to compute)
summary(mhmm_biofam$model)

if (interactive()) {
  # Plotting the model for each cluster (change with Enter)
  plot(mhmm_biofam)
}

Mixture hidden Markov model for the mvad data

Description

A mixture hidden Markov model (MHMM) fitted for the TraMineR::mvad() data.

Format

A mixture hidden Markov model of class mhmm: two clusters including 3 and 4 hidden states. No covariates.

Details

The model is loaded by calling data(mhmm_mvad). It was created with the following code:

data("mvad", package = "TraMineR")

mvad_alphabet <-
  c("employment", "FE", "HE", "joblessness", "school", "training")
mvad_labels <- c("employment", "further education", "higher education",
  "joblessness", "school", "training")
mvad_scodes <- c("EM", "FE", "HE", "JL", "SC", "TR")
mvad_seq <- seqdef(mvad, 15:86, alphabet = mvad_alphabet,
  states = mvad_scodes, labels = mvad_labels, xtstep = 6,
  cpal = colorpalette[[6]])

# Starting values for the emission matrices
emiss_1 <- matrix(
  c(0.01, 0.01, 0.01, 0.01, 0.01, 0.95,
    0.95, 0.01, 0.01, 0.01, 0.01, 0.01,
    0.01, 0.01, 0.01, 0.95, 0.01, 0.01),
  nrow = 3, ncol = 6, byrow = TRUE)

emiss_2 <- matrix(
  c(0.01, 0.01, 0.01, 0.06, 0.90, 0.01,
    0.01, 0.95, 0.01, 0.01, 0.01, 0.01,
    0.01, 0.01, 0.95, 0.01, 0.01, 0.01,
    0.95, 0.01, 0.01, 0.01, 0.01, 0.01),
  nrow = 4, ncol = 6, byrow = TRUE)

# Starting values for the transition matrix

trans_1 <-  matrix(
  c(0.95, 0.03, 0.02,
    0.01, 0.98, 0.01,
    0.01, 0.01, 0.98),
  nrow = 3, ncol = 3, byrow = TRUE)

trans_2 <-  matrix(
  c(0.97, 0.01, 0.01, 0.01,
    0.01, 0.97, 0.01, 0.01,
    0.01, 0.01, 0.97, 0.01,
    0.01, 0.01, 0.01, 0.97),
  nrow = 4, ncol = 4, byrow = TRUE)

# Starting values for initial state probabilities
initial_probs_1 <- c(0.5, 0.25, 0.25)
initial_probs_2 <- c(0.4, 0.4, 0.1, 0.1)

# Building a hidden Markov model with starting values
init_mhmm_mvad <- build_mhmm(observations = mvad_seq,
  transition_probs = list(trans_1, trans_2),
  emission_probs = list(emiss_1, emiss_2),
  initial_probs = list(initial_probs_1, initial_probs_2))

# Fit the model
set.seed(123)
mhmm_mvad <- fit_model(init_mhmm_mvad, control_em = list(restart = list(times = 25)))$model

See Also

Examples of building and fitting MHMMs in build_mhmm() and fit_model(); and TraMineR::mvad() for more information on the data.

Examples

data("mhmm_mvad")

summary(mhmm_mvad)

if (interactive()) {
  # Plotting the model for each cluster (change with Enter)
  plot(mhmm_mvad)
}

Extract Most Probable Cluster for Each Sequence

Description

Extract Most Probable Cluster for Each Sequence

Usage

most_probable_cluster(x, type = "viterbi", hp = NULL)

Arguments

Value

A vector containing the most probable cluster for each sequence.

Interactive Stacked Plots of Multichannel Sequences and/or Most Probable Paths for Mixture Hidden Markov Models

Description

Function mssplot plots stacked sequence plots of observation sequences and/or most probable hidden state paths for each model of the mhmm object (model chosen according to the most probable path).

Usage

mssplot(
  x,
  ask = FALSE,
  which.plots = NULL,
  hidden.paths = NULL,
  plots = "obs",
  type = "d",
  tlim = 0,
  sortv = NULL,
  sort.channel = 1,
  dist.method = "OM",
  with.missing = FALSE,
  missing.color = NULL,
  title = NA,
  title.n = TRUE,
  cex.title = 1,
  title.pos = 1,
  with.legend = "auto",
  ncol.legend = "auto",
  with.missing.legend = "auto",
  legend.prop = 0.3,
  cex.legend = 1,
  hidden.states.colors = "auto",
  hidden.states.labels = "auto",
  xaxis = TRUE,
  xlab = NA,
  xtlab = NULL,
  xlab.pos = 1,
  ylab = "auto",
  hidden.states.title = "Hidden states",
  yaxis = FALSE,
  ylab.pos = "auto",
  cex.lab = 1,
  cex.axis = 1,
  respect_void = TRUE,
  ...
)

Arguments

See Also

build_mhmm() and fit_model() for building and fitting mixture hidden Markov models, hidden_paths() for computing the most probable paths (Viterbi paths) of hidden states, plot.mhmm() for plotting mhmm objects as directed graphs, and colorpalette() for default colors.

Number of Observations in Hidden Markov Model

Description

Extract the number of non-missing observations of HMM. When computing nobs for a multichannel model with $C$ channels, each observed value in a single channel amounts to $1/C$ observation, i.e. a fully observed time point for a single sequence amounts to one observation.

Usage

## S3 method for class 'hmm'
nobs(object, ...)

## S3 method for class 'mhmm'
nobs(object, ...)

## S3 method for class 'nhmm'
nobs(object, ...)

## S3 method for class 'mnhmm'
nobs(object, ...)

Arguments

Plot hidden Markov models

Description

Function plot.hmm plots a directed graph with pie charts of emission probabilities as vertices/nodes.

Usage

## S3 method for class 'hmm'
plot(
  x,
  layout = "horizontal",
  pie = TRUE,
  vertex.size = 40,
  vertex.label = "initial.probs",
  vertex.label.dist = "auto",
  vertex.label.pos = "bottom",
  vertex.label.family = "sans",
  loops = FALSE,
  edge.curved = TRUE,
  edge.label = "auto",
  edge.width = "auto",
  cex.edge.width = 1,
  edge.arrow.size = 1.5,
  edge.label.family = "sans",
  label.signif = 2,
  label.scientific = FALSE,
  label.max.length = 6,
  trim = 1e-15,
  combine.slices = 0.05,
  combined.slice.color = "white",
  combined.slice.label = "others",
  with.legend = "bottom",
  ltext = NULL,
  legend.prop = 0.5,
  cex.legend = 1,
  ncol.legend = "auto",
  cpal = "auto",
  cpal.legend = "auto",
  legend.order = TRUE,
  main = NULL,
  withlegend,
  ...
)

Arguments

See Also

build_hmm() and fit_model() for building and fitting Hidden Markov models, mc_to_sc() for transforming multistate hmm objects into single-channel objects, hmm_biofam() and hmm_mvad() for information on the models used in the examples, and igraph::plot.igraph() for the general plotting function of directed graphs.

Examples

# Multichannel data, left-to-right model

# Loading a HMM of the biofam data
data("hmm_biofam")

# Plotting hmm object
plot(hmm_biofam)

# Plotting HMM with
plot(hmm_biofam,
  # varying curvature of edges
  edge.curved = c(0, -0.7, 0.6, 0.7, 0, -0.7, 0),
  # legend with two columns and less space
  ncol.legend = 2, legend.prop = 0.4,
  # new label for combined slice
  combined.slice.label = "States with probability < 0.05"
)

# Plotting HMM with given coordinates
plot(hmm_biofam,
  # layout given in 2x5 matrix
  # x coordinates in the first column
  # y coordinates in the second column
  layout = matrix(c(
    1, 3, 3, 5, 3,
    0, 0, 1, 0, -1
  ), ncol = 2),
  # larger vertices
  vertex.size = 50,
  # straight edges
  edge.curved = FALSE,
  # thinner edges and arrows
  cex.edge.width = 0.5, edge.arrow.size = 1,
  # varying positions for vertex labels (initial probabilities)
  vertex.label.pos = c(pi, pi / 2, -pi / 2, 0, pi / 2),
  # different legend properties
  with.legend = "top", legend.prop = 0.3, cex.legend = 1.1,
  # Fix axes to the right scale
  xlim = c(0.5, 5.5), ylim = c(-1.5, 1.5), rescale = FALSE,
  # all states (not combining states with small probabilities)
  combine.slices = 0,
  # legend with two columns
  ncol.legend = 2
)

# Plotting HMM with own color palette
plot(hmm_biofam,
  cpal = 1:10,
  # States with emission probability less than 0.2 removed
  combine.slices = 0.2,
  # legend with two columns
  ncol.legend = 2
)

# Plotting HMM without pie graph and with a layout function
require("igraph")
# Setting the seed for a random layout
set.seed(1234)
plot(hmm_biofam,
  # Without pie graph
  pie = FALSE,
  # Using an automatic layout function from igraph
  layout = layout_nicely,
  vertex.size = 30,
  # Straight edges and probabilities of moving to the same state
  edge.curved = FALSE, loops = TRUE,
  # Labels with three significant digits
  label.signif = 3,
  # Fixed edge width
  edge.width = 1,
  # Remove edges with probability less than 0.01
  trim = 0.01,
  # Hidden state names as vertex labels
  vertex.label = "names",
  # Labels insidde vertices
  vertex.label.dist = 0,
  # Fix x-axis (more space on the right-hand side)
  xlim = c(-1, 1.3)
)


# Single-channel data, unrestricted model

# Loading a hidden Markov model of the mvad data (hmm object)
data("hmm_mvad")

# Plotting the HMM
plot(hmm_mvad)

# Checking the order of observed states (needed for the next call)
require(TraMineR)
alphabet(hmm_mvad$observations)

# Plotting the HMM with own legend (note: observation "none" nonexistent in the observations)
plot(hmm_mvad,
  # Override the default order in the legend
  legend.order = FALSE,
  # Colours in the pies (ordered by the alphabet of observations)
  cpal = c("purple", "pink", "brown", "lightblue", "orange", "green"),
  # Colours in the legend (matching to ltext)
  cpal.legend = c("orange", "pink", "brown", "green", "lightblue", "purple", "gray"),
  # Labels in the legend (matching to cpal.legend)
  ltext = c("school", "further educ", "higher educ", "training", "jobless", "employed", "none")
)

require("igraph")
plot(hmm_mvad,
  # Layout in circle (layout function from igraph)
  layout = layout_in_circle,
  # Less curved edges with smaller arrows, no labels
  edge.curved = 0.2, edge.arrow.size = 0.9, edge.label = NA,
  # Positioning vertex labels (initial probabilities)
  vertex.label.pos = c("right", "right", "left", "left", "right"),
  # Less space for the legend
  legend.prop = 0.3
)

Interactive Plotting for Mixed Hidden Markov Model (mhmm)

Description

Function plot.mhmm plots a directed graph of the parameters of each model with pie charts of emission probabilities as vertices/nodes.

Usage

## S3 method for class 'mhmm'
plot(
  x,
  interactive = TRUE,
  ask = FALSE,
  which.plots = NULL,
  nrow = NA,
  ncol = NA,
  byrow = FALSE,
  row.prop = "auto",
  col.prop = "auto",
  layout = "horizontal",
  pie = TRUE,
  vertex.size = 40,
  vertex.label = "initial.probs",
  vertex.label.dist = "auto",
  vertex.label.pos = "bottom",
  vertex.label.family = "sans",
  loops = FALSE,
  edge.curved = TRUE,
  edge.label = "auto",
  edge.width = "auto",
  cex.edge.width = 1,
  edge.arrow.size = 1.5,
  edge.label.family = "sans",
  label.signif = 2,
  label.scientific = FALSE,
  label.max.length = 6,
  trim = 1e-15,
  combine.slices = 0.05,
  combined.slice.color = "white",
  combined.slice.label = "others",
  with.legend = "bottom",
  ltext = NULL,
  legend.prop = 0.5,
  cex.legend = 1,
  ncol.legend = "auto",
  cpal = "auto",
  main = "auto",
  withlegend,
  ...
)

Arguments

References

Helske S. and Helske J. (2019). Mixture Hidden Markov Models for Sequence Data: The seqHMM Package in R, Journal of Statistical Software, 88(3), 1-32. doi:10.18637/jss.v088.i03

See Also

build_mhmm() and fit_model() for building and fitting mixture hidden Markov models; igraph::plot.igraph() for plotting directed graphs; and mhmm_biofam() and mhmm_mvad() for the models used in examples.

Examples

# Loading mixture hidden Markov model (mhmm object)
# of the biofam data
data("mhmm_biofam")

# Plotting only the first cluster
plot(mhmm_biofam, which.plots = 1)

if (interactive()) {
  # Plotting each cluster (change with Enter)
  plot(mhmm_biofam)

  # Choosing the cluster (one at a time)
  plot(mhmm_biofam, ask = TRUE)

  # Loading MHMM of the mvad data
  data("mhmm_mvad")

  # Plotting models in the same graph (in a grid)
  # Note: the plotting window must be high enough!
  set.seed(123)
  plot(mhmm_mvad,
    interactive = FALSE,
    # automatic layout, legend on the right-hand side
    layout = layout_nicely, with.legend = "right",
    # Smaller and less curved edges
    edge.curved = 0.2, cex.edge.width = 0.5, edge.arrow.size = 0.7,
    vertex.label.pos = -4 * pi / 5, vertex.label.dist = 5
  )
}

Stack Multichannel Sequence Plots and/or Most Probable Paths Plots from Hidden Markov Models

Description

Function plot.ssp plots stacked sequence plots from ssp objects defined with ssp().

Usage

## S3 method for class 'ssp'
plot(x, ...)

Arguments

References

Helske S. and Helske J. (2019). Mixture Hidden Markov Models for Sequence Data: The seqHMM Package in R, Journal of Statistical Software, 88(3), 1-32. doi:10.18637/jss.v088.i03

See Also

ssp() for more examples and information on defining the plot before using plot.ssp; ssplot() for straight plotting of ssp objects; and gridplot() for plotting multiple ssp objects.

Plot Colorpalettes

Description

Function plot_colors plots colors and their labels for easy visualization of a colorpalette.

Usage

plot_colors(x, labels = NULL)

Arguments

See Also

See e.g. the colorpalette() data and RColorBrewer package for ready-made color palettes.

Examples

plot_colors(colorpalette[[5]], labels = c("one", "two", "three", "four", "five"))

plot_colors(colorpalette[[10]])

plot_colors(1:7)

plot_colors(c("yellow", "orange", "red", "purple", "blue", "green"))

plot_colors(grDevices::rainbow(15))

Extract Posterior Cluster Probabilities

Description

Extract Posterior Cluster Probabilities

Usage

posterior_cluster_probabilities(x)

Arguments

Value

a data.frame of posterior cluster probabilities for each sequence and cluster.

Posterior Probabilities for Hidden Markov Models

Description

Function posterior_probs computes the posterior probabilities of hidden states of a (mixture) hidden Markov model.

Usage

posterior_probs(model, ...)

## S3 method for class 'hmm'
posterior_probs(model, ...)

## S3 method for class 'mhmm'
posterior_probs(model, ...)

## S3 method for class 'nhmm'
posterior_probs(model, ...)

## S3 method for class 'mnhmm'
posterior_probs(model, ...)

Arguments

Value

A data frame of posterior probabilities for each state and sequence.

Examples

# Load a pre-defined MHMM
data("mhmm_biofam")

# Compute posterior probabilities
pb <- posterior_probs(mhmm_biofam)

Predictions from Non-homogeneous Hidden Markov Models

Description

This function computes the marginal forward predictions for NHMMs and MNHMMs, where the marginalization is (by default) over individuals and time points, weighted by the latent state probabilities.

Usage

## S3 method for class 'nhmm'
predict(
  object,
  newdata,
  newdata2 = NULL,
  condition = NULL,
  type = c("state", "response", "transition", "emission"),
  probs = c(0.025, 0.975),
  boot_idx = FALSE,
  ...
)

## S3 method for class 'mnhmm'
predict(
  object,
  newdata,
  newdata2 = NULL,
  condition = NULL,
  type = c("state", "response", "transition", "emission"),
  probs = c(0.025, 0.975),
  boot_idx = FALSE,
  ...
)

Arguments

Print Method for a Hidden Markov Model

Description

Prints the parameters of a (mixture) hidden Markov model.

Usage

## S3 method for class 'hmm'
print(x, digits = 3, ...)

## S3 method for class 'mhmm'
print(x, digits = 3, ...)

## S3 method for class 'nhmm'
print(x, digits = 3, ...)

## S3 method for class 'mnhmm'
print(x, digits = 3, ...)

## S3 method for class 'summary_mhmm'
print(x, digits = 3, ...)

Arguments

See Also

build_hmm() and fit_model() for building and fitting hidden Markov models.

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

TraMineR

alphabet, seqdef, seqstatf

Convert return code from estimate_nhmm and estimate_mnhmm to text

Description

Convert return code from estimate_nhmm and estimate_mnhmm to text

Usage

return_msg(code)

Arguments

Value

Code translated to informative message.

Reorganize a mixture hidden Markov model to a list of separate hidden Markov models (covariates ignored)

Description

The separate_mhmm function reorganizes the parameters of a mhmm object into a list where each list component is an object of class hmm consisting of the parameters of the corresponding cluster.

Usage

separate_mhmm(model)

Arguments

Value

List with components of class hmm.

See Also

build_mhmm() and fit_model() for building and fitting MHMMs; and mhmm_biofam() for more information on the model used in examples.

Examples

# Loading mixture hidden Markov model (mhmm object)
# of the biofam data
data("mhmm_biofam")

# Separate models for clusters
sep_hmm <- separate_mhmm(mhmm_biofam)

# Plotting the model for the first cluster
plot(sep_hmm[[1]])

Deprecated function(s) in the seqHMM package

Description

These functions still work but will be removed (defunct) in the next version of seqHMM.

Details

• ssplot, ssp, mssplot, plot.ssp. Use stacked_sequence_plot() instead.

• gridplot Use stacked_sequence_plot(), ggseqplot, and patchwork packages instead.

Simulate hidden Markov models

Description

Simulate sequences of observed and hidden states given parameters of a hidden Markov model.

Usage

simulate_hmm(
  n_sequences,
  initial_probs,
  transition_probs,
  emission_probs,
  sequence_length
)

Arguments

Value

A list of state sequence objects of class stslist.

See Also

build_hmm() and fit_model() for building and fitting hidden Markov models; stacked_sequence_plot() for plotting multiple sequence data sets; seqdef() for more information on state sequence objects; and simulate_mhmm() for simulating mixture hidden Markov models.

Examples

# Parameters for the HMM
emission_probs <- matrix(c(0.5, 0.2, 0.5, 0.8), 2, 2)
transition_probs <- matrix(c(5 / 6, 1 / 6, 1 / 6, 5 / 6), 2, 2)
initial_probs <- c(1, 0)

# Setting the seed for simulation
set.seed(1)

# Simulating sequences
sim <- simulate_hmm(
  n_sequences = 10, initial_probs = initial_probs,
  transition_probs = transition_probs,
  emission_probs = emission_probs,
  sequence_length = 20
)

stacked_sequence_plot(sim, sort_by = "mds", type = "i")

Simulate Parameters of Hidden Markov Models

Description

These are helper functions for quick construction of initial values for various model building functions. Mostly useful for global optimization algorithms which do not depend on initial values.

Usage

simulate_initial_probs(n_states, n_clusters = 1, alpha = 1)

simulate_transition_probs(
  n_states,
  n_clusters = 1,
  left_right = FALSE,
  diag_c = 0,
  alpha = 1
)

simulate_emission_probs(n_states, n_symbols, n_clusters = 1, alpha = 1)

Arguments

Simulate Mixture Hidden Markov Models

Description

Simulate sequences of observed and hidden states given the parameters of a mixture hidden Markov model.

Usage

simulate_mhmm(
  n_sequences,
  initial_probs,
  transition_probs,
  emission_probs,
  sequence_length,
  formula = NULL,
  data = NULL,
  coefficients = NULL
)

Arguments

Value

A list of state sequence objects of class stslist.

See Also

build_mhmm() and fit_model() for building and fitting mixture hidden Markov models.

Examples

emission_probs_1 <- matrix(c(0.75, 0.05, 0.25, 0.95), 2, 2)
emission_probs_2 <- matrix(c(0.1, 0.8, 0.9, 0.2), 2, 2)
colnames(emission_probs_1) <- colnames(emission_probs_2) <-
  c("heads", "tails")

transition_probs_1 <- matrix(c(9, 0.1, 1, 9.9) / 10, 2, 2)
transition_probs_2 <- matrix(c(35, 1, 1, 35) / 36, 2, 2)
rownames(emission_probs_1) <- rownames(transition_probs_1) <-
  colnames(transition_probs_1) <- c("coin 1", "coin 2")
rownames(emission_probs_2) <- rownames(transition_probs_2) <-
  colnames(transition_probs_2) <- c("coin 3", "coin 4")

initial_probs_1 <- c(1, 0)
initial_probs_2 <- c(1, 0)

n <- 30
set.seed(123)
covariate_1 <- runif(n)
covariate_2 <- sample(c("A", "B"),
  size = n, replace = TRUE,
  prob = c(0.3, 0.7)
)
dataf <- data.frame(covariate_1, covariate_2)

coefs <- cbind(cluster_1 = c(0, 0, 0), cluster_2 = c(-1.5, 3, -0.7))
rownames(coefs) <- c("(Intercept)", "covariate_1", "covariate_2B")

sim <- simulate_mhmm(
  n = n, initial_probs = list(initial_probs_1, initial_probs_2),
  transition_probs = list(transition_probs_1, transition_probs_2),
  emission_probs = list(emission_probs_1, emission_probs_2),
  sequence_length = 20, formula = ~ covariate_1 + covariate_2,
  data = dataf, coefficients = coefs
)

stacked_sequence_plot(sim, 
  sort_by = "start", sort_channel = "states", type = "i"
)

hmm <- build_mhmm(sim$observations,
  initial_probs = list(initial_probs_1, initial_probs_2),
  transition_probs = list(transition_probs_1, transition_probs_2),
  emission_probs = list(emission_probs_1, emission_probs_2),
  formula = ~ covariate_1 + covariate_2,
  data = dataf
)

fit <- fit_model(hmm)
fit$model

paths <- hidden_paths(fit$model, as_stslist = TRUE)

stacked_sequence_plot(
  list(
    "estimated paths" = paths, 
    "true (simulated)" = sim$states
  ),
  sort_by = "start",
  sort_channel = "true (simulated)",
  type = "i"
)

Simulate Mixture Non-homogeneous Hidden Markov Models

Description

Simulate sequences of observed and hidden states given the parameters of a mixture non-homogeneous hidden Markov model.

Usage

simulate_mnhmm(
  n_states,
  n_clusters,
  emission_formula,
  initial_formula = ~1,
  transition_formula = ~1,
  cluster_formula = ~1,
  data,
  id,
  time,
  coefs = NULL,
  init_sd = 2 * is.null(coefs)
)

Arguments

Value

A list with the model used in simulation as well as the simulated hidden state sequences.

Simulate Non-homogeneous Hidden Markov Models

Description

Simulate sequences of observed and hidden states given the parameters of a non-homogeneous hidden Markov model.

Usage

simulate_nhmm(
  n_states,
  emission_formula,
  initial_formula = ~1,
  transition_formula = ~1,
  data,
  id,
  time,
  coefs = NULL,
  init_sd = 2 * is.null(coefs)
)

Arguments

Value

A list with the model used in simulation as well as the simulated hidden state sequences.

Sort sequences in a sequence object

Description

Sort sequences in a sequence object

Usage

sort_sequences(x, sort_by = "start", sort_channel = 1, dist_method = "OM")

Arguments

Define Arguments for Plotting Multichannel Sequences and/or Most Probable Paths from Hidden Markov Models

Description

Function ssp defines the arguments for plotting with plot.ssp() or gridplot().

Usage

ssp(
  x,
  hidden.paths = NULL,
  plots = "obs",
  type = "d",
  tlim = 0,
  sortv = NULL,
  sort.channel = 1,
  dist.method = "OM",
  with.missing = FALSE,
  missing.color = NULL,
  title = NA,
  title.n = TRUE,
  cex.title = 1,
  title.pos = 1,
  with.legend = "auto",
  ncol.legend = "auto",
  with.missing.legend = "auto",
  legend.prop = 0.3,
  cex.legend = 1,
  hidden.states.colors = "auto",
  hidden.states.labels = "auto",
  xaxis = TRUE,
  xlab = NA,
  xtlab = NULL,
  xlab.pos = 1,
  ylab = "auto",
  hidden.states.title = "Hidden states",
  yaxis = FALSE,
  ylab.pos = "auto",
  cex.lab = 1,
  cex.axis = 1,
  withlegend,
  respect_void = TRUE,
  ...
)

Arguments

Details

This function is deprecated, use stacked_sequence_plot() instead.

Value

Object of class ssp.

Stacked Plots of Multichannel Sequences and/or Most Probable Paths from Hidden Markov Models

Description

Function ssplot plots stacked sequence plots of sequence object created with the seqdef() function or observations and/or most probable paths of hmm objects.

Usage

ssplot(
  x,
  hidden.paths = NULL,
  plots = "obs",
  type = "d",
  tlim = 0,
  sortv = NULL,
  sort.channel = 1,
  dist.method = "OM",
  with.missing = FALSE,
  missing.color = NULL,
  title = NA,
  title.n = TRUE,
  cex.title = 1,
  title.pos = 1,
  with.legend = "auto",
  ncol.legend = "auto",
  with.missing.legend = "auto",
  legend.prop = 0.3,
  cex.legend = 1,
  hidden.states.colors = "auto",
  hidden.states.labels = "auto",
  xaxis = TRUE,
  xlab = NA,
  xtlab = NULL,
  xlab.pos = 1,
  ylab = "auto",
  hidden.states.title = "Hidden states",
  yaxis = FALSE,
  ylab.pos = "auto",
  cex.lab = 1,
  cex.axis = 1,
  respect_void = TRUE,
  ...
)

Arguments

Details

This function is deprecated and will be removed in future versions of seqHMM.

Stacked Sequence Plots of Multichannel Sequences and/or Most Probable Paths from Hidden Markov Models

Description

Function stacked_sequence_plot draws stacked sequence plots of sequence object created with the TraMineR::seqdef function or observations and/or most probable paths of model objects of seqHMM (e.g., hmm and mhmm).

Usage

stacked_sequence_plot(
  x,
  plots = "obs",
  type = "distribution",
  ids,
  sort_by = "none",
  sort_channel,
  dist_method = "OM",
  group = NULL,
  legend_position = "right",
  ...
)

Arguments

Examples

p <- stacked_sequence_plot(
  mhmm_biofam, 
  plots = "both", 
  type = "d", 
  legend_position = c("right", "right", "right", "none")
)
library("ggplot2")
p & theme(plot.margin = unit(c(1, 1, 0, 2), "mm"))

Get State Names of Hidden Markov Model

Description

Get State Names of Hidden Markov Model

Set State Names of Hidden Markov Model

Usage

state_names(object)

## S3 method for class 'hmm'
state_names(object)

## S3 method for class 'mhmm'
state_names(object)

## S3 method for class 'nhmm'
state_names(object)

## S3 method for class 'mnhmm'
state_names(object)

state_names(object) <- value

## S3 replacement method for class 'hmm'
state_names(object) <- value

## S3 replacement method for class 'mhmm'
state_names(object) <- value

## S3 replacement method for class 'nhmm'
state_names(object) <- value

## S3 replacement method for class 'mnhmm'
state_names(object) <- value

Arguments

Value

A character vector containing the state names, or a list of such vectors in case of mixture models.

The original object with updated state names.

Summary method for mixture hidden Markov models

Description

Function summary.mhmm gives a summary of a mixture hidden Markov model.

Usage

## S3 method for class 'mhmm'
summary(object, parameters = FALSE, conditional_se = TRUE, ...)

Arguments

Details

The summary.mhmm function computes features from a mixture hidden Markov model and stores them as a list. A print method prints summaries of these: log-likelihood and BIC, coefficients and standard errors of covariates, means of prior cluster probabilities, and information on most probable clusters.

Value

• transition_probs
  Transition probabilities. Only returned if parameters = TRUE.

• emission_probs
  Emission probabilities. Only returned if parameters = TRUE.

• initial_probs
  Initial state probabilities. Only returned if parameters = TRUE.

• logLik
  Log-likelihood.

• BIC
  Bayesian information criterion.

• most_probable_cluster
  The most probable cluster according to posterior probabilities.

• coefficients
  Coefficients of covariates.

• vcov
  Variance-covariance matrix of coefficients.

• prior_cluster_probabilities
  Prior cluster probabilities (mixing proportions) given the covariates.

• posterior_cluster_probabilities
  Posterior cluster membership probabilities.

• classification_table
  Cluster probabilities (columns) by the most probable cluster (rows).

See Also

build_mhmm() and fit_model() for building and fitting mixture hidden Markov models; and mhmm_biofam() for information on the model used in examples.

Examples

# Loading mixture hidden Markov model (mhmm object)
# of the biofam data
data("mhmm_biofam")

# Model summary
summary(mhmm_biofam)

Summary method for mixture non-homogenous hidden Markov models

Description

Summary method for mixture non-homogenous hidden Markov models

Usage

## S3 method for class 'mnhmm'
summary(object, ...)

Arguments

Trim Small Probabilities of Hidden Markov Model

Description

Function trim_model tries to set small insignificant probabilities to zero without decreasing the likelihood.

Usage

trim_model(
  model,
  maxit = 0,
  return_loglik = FALSE,
  zerotol = 1e-08,
  verbose = TRUE,
  ...
)

Arguments

See Also

build_hmm() and fit_model() for building and fitting hidden Markov models; and hmm_biofam() for information on the model used in the example.

Examples

data("hmm_biofam")

# Testing if changing parameter values smaller than 1e-03 to zero
# leads to improved log-likelihood.
hmm_trim <- trim_model(hmm_biofam, zerotol = 1e-03, maxit = 10)

Update Covariate Values of NHMM

Description

This function can be used to replace original covariate values of NHMMs. The responses, model formulae and estimated coefficients are not altered.

Usage

## S3 method for class 'nhmm'
update(object, newdata, ...)

## S3 method for class 'mnhmm'
update(object, newdata, ...)

Arguments

Variance-Covariance Matrix for Coefficients of Covariates of Mixture Hidden Markov Model

Description

Returns the asymptotic covariances matrix of maximum likelihood estimates of the coefficients corresponding to the explanatory variables of the model.

Usage

## S3 method for class 'mhmm'
vcov(object, conditional = TRUE, threads = 1, log_space = TRUE, ...)

Arguments

Details

The conditional standard errors are computed using analytical formulas by assuming that the coefficient estimates are not correlated with other model parameter estimates (or that the other parameters are assumed to be fixed). This often underestimates the true standard errors, but is substantially faster approach for preliminary analysis. The non-conditional standard errors are based on the numerical approximation of the full Hessian of the coefficients and the model parameters corresponding to nonzero probabilities. Computing the non-conditional standard errors can be slow for large models as the Jacobian of analytical gradients is computed using finite difference approximation.

Alternatively, by using the non-homogeneous model via estimate_mnhmm you can compute the standard errors of the coefficients using the bootstrap method.

Value

Matrix containing the variance-covariance matrix of coefficients.