Package 'seqHMM'

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 [email protected] (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

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.

Description

The model estimation for each bootstrap sample uses the same method and tolerances as the original fit. If you want to change these, you can modify the elements of the input model such as model$estimation_results$method and model$controls before passing it to bootstrap_coefs().

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

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.

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.

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
)

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)

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)
  )
)

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)

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)

Description

Get Cluster Names from Mixture HMMs

Usage

cluster_names(object)

Arguments

Value

A character vector containing the cluster names.

Description

Set Cluster Names for Mixture Models

Usage

cluster_names(object) <- value

Arguments

Value

The modified object with updated cluster names.

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).

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]]

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 = "id", time = "time", responses, ...)

Arguments

Examples

# Convert multichannel stslist to long format data.table
biofam_seqs <- mhmm_biofam$observations # list of three stslist objects
d <- stslist_to_data(
  biofam_seqs, id = "Individual", time = "Age",
  responses = c("Marriage", "Children", "Residence")
)
head(d)

# Convert back to stslist
biofam_seqs2 <- data_to_stslist(
  d, id = "Individual", time = "Age", 
  responses = c("Marriage", "Children", "Residence")
)
all.equal(biofam_seqs, biofam_seqs2, check.attributes = FALSE)

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 = NULL,
  restarts = 0L,
  method = "EM-DNM",
  bound = Inf,
  control_restart = list(),
  control_mstep = list(),
  check_rank = NULL,
  ...
)

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 +
  N.Eastern + Southern + S.Eastern + Western,
  emission_formula = activity ~ male + catholic + gcse5eq,
  initial_formula = ~ 1, 
  transition_formula = ~ male + gcse5eq
  )

## End(Not run)

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 = NULL,
  restarts = 0L,
  method = "EM-DNM",
  bound = Inf,
  control_restart = list(),
  control_mstep = list(),
  check_rank = NULL,
  ...
)

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$ onward 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. Using argument save_all_solutions = TRUE allows to save the optimization results of all restarts, which is mostly useful for debugging purposes or if you are interested in seeing how much the parameter estimates vary due to variation in initial values.

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, 1e5 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.

By default, EM algorithm uses SQUAREM acceleration (Varadhan and Roland, 2008) to improve the convergence speed near optimum. This can be turned of by using control argument use_squarem = FALSE.

Value

Object of class nhmm or fanhmm.

References

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

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

Varadhan R and Roland C (2008). Simple and Globally Convergent Methods for Accelerating the Convergence of Any EM Algorithm. Scandinavian Journal of Statistics, 35: 335-353. https://doi.org/10.1111/j.1467-9469.2007.00585.x

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)

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)

Description

Function fit_model estimates the parameters of mixture hidden Markov models and its restricted variants using maximum likelihood. Initial values for estimation are taken from the corresponding elements 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 = FALSE,
  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 should 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

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)

Description

Prior cluster probability refers to the probability of sequence $i$ belonging to cluster $k$ without conditioning on the observed sequence, i.e., $P(i \in M_k | X_i)$, where $M_k$ is cluster $k$ and $X_i$ are the covariates for individual $i$. The posterior cluster probability refers to the probability of sequence $i$ belonging to cluster $k$ given the observed sequence, i.e., $P(i \in M_k | y_i, X_i)$, where $y_i$ is the observed sequence for individual $i$.

Usage

get_cluster_probs(model, ...)

## S3 method for class 'mnhmm'
get_cluster_probs(model, type = "prior", ids = NULL, ...)

## S3 method for class 'mhmm'
get_cluster_probs(model, type = "prior", ...)

Arguments

See Also

posterior_cluster_probabilities().

Description

Extract the Emission Probabilities of Hidden Markov Model

Usage

get_emission_probs(model, ...)

## S3 method for class 'nhmm'
get_emission_probs(model, ids = NULL, ...)

## S3 method for class 'mnhmm'
get_emission_probs(model, ids = NULL, ...)

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

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

Arguments

Description

Extract the Initial State Probabilities of Hidden Markov Model

Usage

get_initial_probs(model, ...)

## S3 method for class 'nhmm'
get_initial_probs(model, ids = NULL, ...)

## S3 method for class 'mnhmm'
get_initial_probs(model, ids = NULL, ...)

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

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

Arguments

Value

A data.table containing the initial state probabilities for each sequence defined by ids.

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

Description

Extract the State Transition Probabilities of Hidden Markov Model

Usage

get_transition_probs(model, ...)

## S3 method for class 'nhmm'
get_transition_probs(model, ids = NULL, ...)

## S3 method for class 'mnhmm'
get_transition_probs(model, ids = NULL, ...)

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

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

Arguments

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

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)

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)

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)

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)

Description

Log-likelihood of a Hidden Markov Model

Usage

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

## S3 method for class 'mhmm'
logLik(object, partials = FALSE, threads = 1, log_space = 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.

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.

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)

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)

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)
}

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)
}

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.

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,
  ...
)