diff --git a/NAMESPACE b/NAMESPACE
index 2d24976a..50bd5147 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -65,8 +65,10 @@ export(gridplot)
 export(hidden_paths)
 export(mc_to_sc)
 export(mc_to_sc_data)
+export(most_probable_cluster)
 export(mssplot)
 export(plot_colors)
+export(posterior_cluster_probabilities)
 export(posterior_probs)
 export(separate_mhmm)
 export(seqdef)
diff --git a/R/most_probable_cluster.R b/R/most_probable_cluster.R
index cc6e07b0..71d6bc10 100644
--- a/R/most_probable_cluster.R
+++ b/R/most_probable_cluster.R
@@ -1,6 +1,23 @@
-most_probable_cluster <- function(x, type = "viterbi", hp) {
+#' Extract Most Probable Cluster for Each Sequence
+#' 
+#' @param x An object of class `mhmm` or `mnhmm`.
+#' @param type A character string specifying the method to use. Either
+#' `"viterbi"` (default) or `"posterior"`. Former uses the most probable hidden 
+#' path to determine the cluster membership for each sequence, while the latter 
+#' finds the cluster which has the largest sum of posterior probabilities of 
+#' states of that cluster. 
+#' @param hp An output from [hidden_paths()] function. Only used in case of 
+#' `type = "viterbi"`. If missing, hidden paths will be computed using `x`.
+#' @return A vector containing the most probable cluster for each sequence.
+#' @export
+most_probable_cluster <- function(x, type = "viterbi", hp = NULL) {
+  stopifnot_(
+    inherits(x, "mhmm") || inherits(x, "mnhmm"),
+    "Argument {.arg x} must be a {.cls mhmm} or {.cls mnhmm} object."
+  )
   type <- match.arg(type, c("viterbi", "posterior"))
   if (type == "viterbi") {
+    if (is.null(hp)) hp <- hidden_paths(x)
     mm <- NULL
     state_names <- unlist(x$state_names)
     clusters <- numeric(x$n_sequences)
@@ -11,15 +28,24 @@ most_probable_cluster <- function(x, type = "viterbi", hp) {
       } else {
         idx <- which(grepl(paste0(x$cluster_names[i], ": "), hp[, 1]))
       }
-      clusters[idx] <- x$cluster_names[i]
+      clusters[idx] <- i
     }
   } else {
     clusters <- apply(posterior_cluster_probabilities(x), 1, which.max)
   }
   factor(clusters, levels = seq_len(x$n_clusters), labels = x$cluster_names)
 }
-
+#' Extract Posterior Cluster Probabilities
+#' 
+#' @param x An object of class `mhmm` or `mnhmm`.
+#' @return matrix of posterior cluster probabilities for each sequence and 
+#' cluster.
+#' @export
 posterior_cluster_probabilities <- function(x) {
+  stopifnot_(
+    inherits(x, "mhmm") || inherits(x, "mnhmm"),
+    "Argument {.arg x} must be a {.cls mhmm} or {.cls mnhmm} object."
+  )
   pp <- posterior_probs(x, as_data_frame = FALSE)
   posterior_cluster_probabilities <- matrix(0, x$n_sequences, x$n_clusters)
   n_states <- rep(x$n_states, length.out = x$n_clusters)
diff --git a/R/mssplot-deprecated.R b/R/mssplot-deprecated.R
index 530a8260..8fc10f6c 100644
--- a/R/mssplot-deprecated.R
+++ b/R/mssplot-deprecated.R
@@ -151,23 +151,9 @@
 #' @param respect_void  If `TRUE` (default), states at the time points
 #' corresponding to TraMineR's void in the observed sequences are set to void
 #' in the hidden state sequences as well.
-#'
 #' @param ... Other arguments to be passed on to
 #'   [TraMineR::seqplot()].
 #'
-#' @examples
-#' # Loading mixture hidden Markov model (mhmm object)
-#' # of the biofam data
-#' data("mhmm_biofam")
-#'
-#' # Plotting the first cluster only
-#' mssplot(mhmm_biofam, which.plots = 1)
-#'
-#' if (interactive()) {
-#'   # Interactive plot
-#'   mssplot(mhmm_biofam)
-#' }
-#'
 #' @seealso [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,
diff --git a/R/plot.ssp-deprecated.R b/R/plot.ssp-deprecated.R
index ab9e2c2f..4b0d17d8 100644
--- a/R/plot.ssp-deprecated.R
+++ b/R/plot.ssp-deprecated.R
@@ -14,25 +14,6 @@
 #'
 #' @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
-#'
-#' @examples
-#'
-#' data("biofam3c")
-#'
-#' ## Building sequence objects
-#' child_seq <- seqdef(biofam3c$children, start = 15)
-#' marr_seq <- seqdef(biofam3c$married, start = 15)
-#' left_seq <- seqdef(biofam3c$left, start = 15)
-#'
-#' ## Choosing colors
-#' attr(child_seq, "cpal") <- c("#66C2A5", "#FC8D62")
-#' attr(marr_seq, "cpal") <- c("#AB82FF", "#E6AB02", "#E7298A")
-#' attr(left_seq, "cpal") <- c("#A6CEE3", "#E31A1C")
-#'
-#'
-#' # Plotting state distribution plots of observations
-#' ssp1 <- ssp(list(child_seq, marr_seq, left_seq))
-#' plot(ssp1)
 plot.ssp <- function(x, ...) {
   .Deprecated("stacked_sequence_plot")
   plot.new()
diff --git a/R/simulate_mhmm.R b/R/simulate_mhmm.R
index f338ecfe..3a156244 100644
--- a/R/simulate_mhmm.R
+++ b/R/simulate_mhmm.R
@@ -69,9 +69,8 @@
 #'   data = dataf, coefficients = coefs
 #' )
 #'
-#' ssplot(sim$observations,
-#'   hidden.paths = sim$states, plots = "both",
-#'   sortv = "from.start", sort.channel = 0, type = "I"
+#' stacked_sequence_plot(sim, 
+#'   sort_by = "start", sort_channel = "states", type = "i"
 #' )
 #'
 #' hmm <- build_mhmm(sim$observations,
@@ -87,10 +86,14 @@
 #'
 #' paths <- hidden_paths(fit$model)
 #'
-#' ssplot(list(estimates = paths, true = sim$states),
-#'   sortv = "from.start",
-#'   sort.channel = 2, ylab = c("estimated paths", "true (simulated)"),
-#'   type = "I"
+#' stacked_sequence_plot(
+#'   list(
+#'     "estimated paths" = paths, 
+#'     "true (simulated)" = sim$states
+#'   ),
+#'   sort_by = "start",
+#'   sort_channel = "true (simulated)",
+#'   type = "i"
 #' )
 #'
 simulate_mhmm <- function(
diff --git a/R/ssp-deprecated.R b/R/ssp-deprecated.R
index 0e572119..0f54485b 100644
--- a/R/ssp-deprecated.R
+++ b/R/ssp-deprecated.R
@@ -154,91 +154,7 @@
 #' corresponding to TraMineR's void in the observed sequences are set to void
 #' in the hidden state sequences as well.
 #' @param ... Other arguments to be passed on to [TraMineR::seqplot()].
-#'
 #' @return Object of class `ssp`.
-#'
-#' @seealso [plot.ssp()] for plotting objects created with
-#'   the `ssp` function; [gridplot()] for plotting multiple `ssp`
-#'   objects; [build_hmm()] and [fit_model()] for building and
-#'   fitting hidden Markov models; [hidden_paths()] for
-#'   computing the most probable paths of hidden states; and [biofam3c()] and
-#'   [hmm_biofam()] for information on the data and model used in the example.
-#'
-#' @examples
-#' data("biofam3c")
-#'
-#' ## Building sequence objects
-#' child_seq <- seqdef(biofam3c$children, start = 15)
-#' marr_seq <- seqdef(biofam3c$married, start = 15)
-#' left_seq <- seqdef(biofam3c$left, start = 15)
-#'
-#' ## Choosing colors
-#' attr(child_seq, "cpal") <- c("#66C2A5", "#FC8D62")
-#' attr(marr_seq, "cpal") <- c("#AB82FF", "#E6AB02", "#E7298A")
-#' attr(left_seq, "cpal") <- c("#A6CEE3", "#E31A1C")
-#'
-#'
-#' # Defining the plot for state distribution plots of observations
-#' ssp1 <- ssp(list(
-#'   "Parenthood" = child_seq, "Marriage" = marr_seq,
-#'   "Residence" = left_seq
-#' ))
-#' # Plotting ssp1
-#' plot(ssp1)
-#'
-#' \dontrun{
-#' # Defining the plot for sequence index plots of observations
-#' ssp2 <- ssp(
-#'   list(child_seq, marr_seq, left_seq),
-#'   type = "I", plots = "obs",
-#'   # Sorting subjects according to the beginning of the 2nd channel (marr_seq)
-#'   sortv = "from.start", sort.channel = 2,
-#'   # Controlling the size, positions, and names for channel labels
-#'   ylab.pos = c(1, 2, 1), cex.lab = 1, ylab = c("Children", "Married", "Residence"),
-#'   # Plotting without legend
-#'   with.legend = FALSE
-#' )
-#' plot(ssp2)
-#'
-#' # Plotting hidden Markov models
-#'
-#' # Loading data
-#' data("hmm_biofam")
-#'
-#' # Plotting observations and most probable hidden states paths
-#' ssp3 <- ssp(
-#'   hmm_biofam,
-#'   type = "I", plots = "both",
-#'   # Sorting according to multidimensional scaling of hidden states paths
-#'   sortv = "mds.hidden",
-#'   # Controlling title
-#'   title = "Biofam", cex.title = 1.5,
-#'   # Labels for x axis and tick marks
-#'   xtlab = 15:30, xlab = "Age"
-#' )
-#' plot(ssp3)
-#'
-#' # Computing the most probable paths of hidden states
-#' hid <- hidden_paths(hmm_biofam)
-#' # Giving names for hidden states
-#' library(TraMineR)
-#' alphabet(hid) <- paste("Hidden state", 1:5)
-#'
-#' # Plotting observations and hidden state paths
-#' ssp4 <- ssp(
-#'   hmm_biofam,
-#'   type = "I", plots = "hidden.paths",
-#'   # Sequence object of most probable paths
-#'   hidden.paths = hid,
-#'   # Sorting according to the end of hidden state paths
-#'   sortv = "from.end", sort.channel = 0,
-#'   # Contolling legend position, type, and proportion
-#'   with.legend = "bottom.combined", legend.prop = 0.15,
-#'   # Plotting without title and y label
-#'   title = FALSE, ylab = FALSE
-#' )
-#' plot(ssp4)
-#' }
 ssp <- function(
     x, hidden.paths = NULL,
     plots = "obs", type = "d", tlim = 0,
diff --git a/R/ssplot-deprecated.R b/R/ssplot-deprecated.R
index db68513f..629a7cd5 100644
--- a/R/ssplot-deprecated.R
+++ b/R/ssplot-deprecated.R
@@ -155,82 +155,6 @@
 #'
 #' @param ... Other arguments to be passed on to
 #'   [TraMineR::seqplot()].
-#'
-#' @examples
-#' data("biofam3c")
-#'
-#' # Creating sequence objects
-#' child_seq <- seqdef(biofam3c$children, start = 15)
-#' marr_seq <- seqdef(biofam3c$married, start = 15)
-#' left_seq <- seqdef(biofam3c$left, start = 15)
-#'
-#' ## Choosing colors
-#' attr(child_seq, "cpal") <- c("#66C2A5", "#FC8D62")
-#' attr(marr_seq, "cpal") <- c("#AB82FF", "#E6AB02", "#E7298A")
-#' attr(left_seq, "cpal") <- c("#A6CEE3", "#E31A1C")
-#'
-#'
-#' # Plotting state distribution plots of observations
-#' ssplot(list(
-#'   "Children" = child_seq, "Marriage" = marr_seq,
-#'   "Residence" = left_seq
-#' ))
-#'
-#' \dontrun{
-#' # Plotting sequence index plots of observations
-#' ssplot(
-#'   list(child_seq, marr_seq, left_seq),
-#'   type = "I",
-#'   # Sorting subjects according to the beginning of the 2nd channel (marr_seq)
-#'   sortv = "from.start", sort.channel = 2,
-#'   # Controlling the size, positions, and names for channel labels
-#'   ylab.pos = c(1, 2, 1), cex.lab = 1, ylab = c("Children", "Married", "Residence"),
-#'   # Plotting without legend
-#'   with.legend = FALSE
-#' )
-#'
-#' # Plotting hidden Markov models
-#'
-#' # Loading a ready-made HMM for the biofam data
-#' data("hmm_biofam")
-#'
-#' # Plotting observations and hidden states paths
-#' ssplot(
-#'   hmm_biofam,
-#'   type = "I", plots = "both",
-#'   # Sorting according to multidimensional scaling of hidden states paths
-#'   sortv = "mds.hidden",
-#'   ylab = c("Children", "Married", "Left home"),
-#'   # Controlling title
-#'   title = "Biofam", cex.title = 1.5,
-#'   # Labels for x axis and tick marks
-#'   xtlab = 15:30, xlab = "Age"
-#' )
-#'
-#' # Computing the most probable paths of hidden states
-#' hidden.paths <- hidden_paths(hmm_biofam)
-#' hidden.paths_seq <- seqdef(hidden.paths, labels = paste("Hidden state", 1:5))
-#'
-#' # Plotting observations and hidden state paths
-#' ssplot(
-#'   hmm_biofam,
-#'   type = "I", plots = "hidden.paths",
-#'   # Sequence object of most probable paths
-#'   hidden.paths = hidden.paths_seq,
-#'   # Sorting according to the end of hidden state paths
-#'   sortv = "from.end", sort.channel = 0,
-#'   # Contolling legend position, type, and proportion
-#'   with.legend = "bottom", legend.prop = 0.15,
-#'   # Plotting without title and y label
-#'   title = FALSE, ylab = FALSE
-#' )
-#' }
-#' @seealso [ssp()] for creating `ssp` objects and [plot.ssp()]
-#' and [gridplot()] for plotting these;
-#' [build_hmm()] and [fit_model()] for building and
-#' fitting hidden Markov models; [hidden_paths()] for
-#' computing the most probable paths of hidden states; and [biofam3c()]
-#' [hmm_biofam()] for information on the data and model used in the example.
 ssplot <- function(x, hidden.paths = NULL,
                    plots = "obs", type = "d", tlim = 0,
                    sortv = NULL, sort.channel = 1, dist.method = "OM",
diff --git a/man/most_probable_cluster.Rd b/man/most_probable_cluster.Rd
new file mode 100644
index 00000000..932f7f35
--- /dev/null
+++ b/man/most_probable_cluster.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/most_probable_cluster.R
+\name{most_probable_cluster}
+\alias{most_probable_cluster}
+\title{Extract Most Probable Cluster for Each Sequence}
+\usage{
+most_probable_cluster(x, type = "viterbi", hp = NULL)
+}
+\arguments{
+\item{x}{An object of class \code{mhmm} or \code{mnhmm}.}
+
+\item{type}{A character string specifying the method to use. Either
+\code{"viterbi"} (default) or \code{"posterior"}. Former uses the most probable hidden
+path to determine the cluster membership for each sequence, while the latter
+finds the cluster which has the largest sum of posterior probabilities of
+states of that cluster.}
+
+\item{hp}{An output from \code{\link[=hidden_paths]{hidden_paths()}} function. Only used in case of
+\code{type = "viterbi"}. If missing, hidden paths will be computed using \code{x}.}
+}
+\value{
+A vector containing the most probable cluster for each sequence.
+}
+\description{
+Extract Most Probable Cluster for Each Sequence
+}
diff --git a/man/mssplot.Rd b/man/mssplot.Rd
index 134abe93..7cfa2ad9 100644
--- a/man/mssplot.Rd
+++ b/man/mssplot.Rd
@@ -197,20 +197,6 @@ in the hidden state sequences as well.}
 Function \code{mssplot} plots stacked sequence plots of observation sequences
 and/or most probable hidden state paths for each model of the \code{mhmm}
 object (model chosen according to the most probable path).
-}
-\examples{
-# Loading mixture hidden Markov model (mhmm object)
-# of the biofam data
-data("mhmm_biofam")
-
-# Plotting the first cluster only
-mssplot(mhmm_biofam, which.plots = 1)
-
-if (interactive()) {
-  # Interactive plot
-  mssplot(mhmm_biofam)
-}
-
 }
 \seealso{
 \code{\link[=build_mhmm]{build_mhmm()}} and \code{\link[=fit_model]{fit_model()}} for building and
diff --git a/man/plot.ssp.Rd b/man/plot.ssp.Rd
index b2662bd8..995d394b 100644
--- a/man/plot.ssp.Rd
+++ b/man/plot.ssp.Rd
@@ -16,25 +16,6 @@ Models}
 Function \code{plot.ssp} plots stacked sequence plots from \code{ssp} objects defined with
 \code{\link[=ssp]{ssp()}}.
 }
-\examples{
-
-data("biofam3c")
-
-## Building sequence objects
-child_seq <- seqdef(biofam3c$children, start = 15)
-marr_seq <- seqdef(biofam3c$married, start = 15)
-left_seq <- seqdef(biofam3c$left, start = 15)
-
-## Choosing colors
-attr(child_seq, "cpal") <- c("#66C2A5", "#FC8D62")
-attr(marr_seq, "cpal") <- c("#AB82FF", "#E6AB02", "#E7298A")
-attr(left_seq, "cpal") <- c("#A6CEE3", "#E31A1C")
-
-
-# Plotting state distribution plots of observations
-ssp1 <- ssp(list(child_seq, marr_seq, left_seq))
-plot(ssp1)
-}
 \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
diff --git a/man/posterior_cluster_probabilities.Rd b/man/posterior_cluster_probabilities.Rd
new file mode 100644
index 00000000..b63dc259
--- /dev/null
+++ b/man/posterior_cluster_probabilities.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/most_probable_cluster.R
+\name{posterior_cluster_probabilities}
+\alias{posterior_cluster_probabilities}
+\title{Extract Posterior Cluster Probabilities}
+\usage{
+posterior_cluster_probabilities(x)
+}
+\arguments{
+\item{x}{An object of class \code{mhmm} or \code{mnhmm}.}
+}
+\value{
+matrix of posterior cluster probabilities for each sequence and
+cluster.
+}
+\description{
+Extract Posterior Cluster Probabilities
+}
diff --git a/man/simulate_mhmm.Rd b/man/simulate_mhmm.Rd
index 518178c4..e8056014 100644
--- a/man/simulate_mhmm.Rd
+++ b/man/simulate_mhmm.Rd
@@ -89,9 +89,8 @@ sim <- simulate_mhmm(
   data = dataf, coefficients = coefs
 )
 
-ssplot(sim$observations,
-  hidden.paths = sim$states, plots = "both",
-  sortv = "from.start", sort.channel = 0, type = "I"
+stacked_sequence_plot(sim, 
+  sort_by = "start", sort_channel = "states", type = "i"
 )
 
 hmm <- build_mhmm(sim$observations,
@@ -107,10 +106,14 @@ fit$model
 
 paths <- hidden_paths(fit$model)
 
-ssplot(list(estimates = paths, true = sim$states),
-  sortv = "from.start",
-  sort.channel = 2, ylab = c("estimated paths", "true (simulated)"),
-  type = "I"
+stacked_sequence_plot(
+  list(
+    "estimated paths" = paths, 
+    "true (simulated)" = sim$states
+  ),
+  sort_by = "start",
+  sort_channel = "true (simulated)",
+  type = "i"
 )
 
 }
diff --git a/man/ssp.Rd b/man/ssp.Rd
index 66b43699..911c2810 100644
--- a/man/ssp.Rd
+++ b/man/ssp.Rd
@@ -200,87 +200,3 @@ Object of class \code{ssp}.
 Function \code{ssp} defines the arguments for plotting with
 \code{\link[=plot.ssp]{plot.ssp()}} or \code{\link[=gridplot]{gridplot()}}.
 }
-\examples{
-data("biofam3c")
-
-## Building sequence objects
-child_seq <- seqdef(biofam3c$children, start = 15)
-marr_seq <- seqdef(biofam3c$married, start = 15)
-left_seq <- seqdef(biofam3c$left, start = 15)
-
-## Choosing colors
-attr(child_seq, "cpal") <- c("#66C2A5", "#FC8D62")
-attr(marr_seq, "cpal") <- c("#AB82FF", "#E6AB02", "#E7298A")
-attr(left_seq, "cpal") <- c("#A6CEE3", "#E31A1C")
-
-
-# Defining the plot for state distribution plots of observations
-ssp1 <- ssp(list(
-  "Parenthood" = child_seq, "Marriage" = marr_seq,
-  "Residence" = left_seq
-))
-# Plotting ssp1
-plot(ssp1)
-
-\dontrun{
-# Defining the plot for sequence index plots of observations
-ssp2 <- ssp(
-  list(child_seq, marr_seq, left_seq),
-  type = "I", plots = "obs",
-  # Sorting subjects according to the beginning of the 2nd channel (marr_seq)
-  sortv = "from.start", sort.channel = 2,
-  # Controlling the size, positions, and names for channel labels
-  ylab.pos = c(1, 2, 1), cex.lab = 1, ylab = c("Children", "Married", "Residence"),
-  # Plotting without legend
-  with.legend = FALSE
-)
-plot(ssp2)
-
-# Plotting hidden Markov models
-
-# Loading data
-data("hmm_biofam")
-
-# Plotting observations and most probable hidden states paths
-ssp3 <- ssp(
-  hmm_biofam,
-  type = "I", plots = "both",
-  # Sorting according to multidimensional scaling of hidden states paths
-  sortv = "mds.hidden",
-  # Controlling title
-  title = "Biofam", cex.title = 1.5,
-  # Labels for x axis and tick marks
-  xtlab = 15:30, xlab = "Age"
-)
-plot(ssp3)
-
-# Computing the most probable paths of hidden states
-hid <- hidden_paths(hmm_biofam)
-# Giving names for hidden states
-library(TraMineR)
-alphabet(hid) <- paste("Hidden state", 1:5)
-
-# Plotting observations and hidden state paths
-ssp4 <- ssp(
-  hmm_biofam,
-  type = "I", plots = "hidden.paths",
-  # Sequence object of most probable paths
-  hidden.paths = hid,
-  # Sorting according to the end of hidden state paths
-  sortv = "from.end", sort.channel = 0,
-  # Contolling legend position, type, and proportion
-  with.legend = "bottom.combined", legend.prop = 0.15,
-  # Plotting without title and y label
-  title = FALSE, ylab = FALSE
-)
-plot(ssp4)
-}
-}
-\seealso{
-\code{\link[=plot.ssp]{plot.ssp()}} for plotting objects created with
-the \code{ssp} function; \code{\link[=gridplot]{gridplot()}} for plotting multiple \code{ssp}
-objects; \code{\link[=build_hmm]{build_hmm()}} and \code{\link[=fit_model]{fit_model()}} for building and
-fitting hidden Markov models; \code{\link[=hidden_paths]{hidden_paths()}} for
-computing the most probable paths of hidden states; and \code{\link[=biofam3c]{biofam3c()}} and
-\code{\link[=hmm_biofam]{hmm_biofam()}} for information on the data and model used in the example.
-}
diff --git a/man/ssplot.Rd b/man/ssplot.Rd
index f1cb2c42..8c67ab84 100644
--- a/man/ssplot.Rd
+++ b/man/ssplot.Rd
@@ -198,81 +198,3 @@ probable paths of \code{hmm} objects.
 \details{
 This function is deprecated and will be removed in future versions of \code{seqHMM}.
 }
-\examples{
-data("biofam3c")
-
-# Creating sequence objects
-child_seq <- seqdef(biofam3c$children, start = 15)
-marr_seq <- seqdef(biofam3c$married, start = 15)
-left_seq <- seqdef(biofam3c$left, start = 15)
-
-## Choosing colors
-attr(child_seq, "cpal") <- c("#66C2A5", "#FC8D62")
-attr(marr_seq, "cpal") <- c("#AB82FF", "#E6AB02", "#E7298A")
-attr(left_seq, "cpal") <- c("#A6CEE3", "#E31A1C")
-
-
-# Plotting state distribution plots of observations
-ssplot(list(
-  "Children" = child_seq, "Marriage" = marr_seq,
-  "Residence" = left_seq
-))
-
-\dontrun{
-# Plotting sequence index plots of observations
-ssplot(
-  list(child_seq, marr_seq, left_seq),
-  type = "I",
-  # Sorting subjects according to the beginning of the 2nd channel (marr_seq)
-  sortv = "from.start", sort.channel = 2,
-  # Controlling the size, positions, and names for channel labels
-  ylab.pos = c(1, 2, 1), cex.lab = 1, ylab = c("Children", "Married", "Residence"),
-  # Plotting without legend
-  with.legend = FALSE
-)
-
-# Plotting hidden Markov models
-
-# Loading a ready-made HMM for the biofam data
-data("hmm_biofam")
-
-# Plotting observations and hidden states paths
-ssplot(
-  hmm_biofam,
-  type = "I", plots = "both",
-  # Sorting according to multidimensional scaling of hidden states paths
-  sortv = "mds.hidden",
-  ylab = c("Children", "Married", "Left home"),
-  # Controlling title
-  title = "Biofam", cex.title = 1.5,
-  # Labels for x axis and tick marks
-  xtlab = 15:30, xlab = "Age"
-)
-
-# Computing the most probable paths of hidden states
-hidden.paths <- hidden_paths(hmm_biofam)
-hidden.paths_seq <- seqdef(hidden.paths, labels = paste("Hidden state", 1:5))
-
-# Plotting observations and hidden state paths
-ssplot(
-  hmm_biofam,
-  type = "I", plots = "hidden.paths",
-  # Sequence object of most probable paths
-  hidden.paths = hidden.paths_seq,
-  # Sorting according to the end of hidden state paths
-  sortv = "from.end", sort.channel = 0,
-  # Contolling legend position, type, and proportion
-  with.legend = "bottom", legend.prop = 0.15,
-  # Plotting without title and y label
-  title = FALSE, ylab = FALSE
-)
-}
-}
-\seealso{
-\code{\link[=ssp]{ssp()}} for creating \code{ssp} objects and \code{\link[=plot.ssp]{plot.ssp()}}
-and \code{\link[=gridplot]{gridplot()}} for plotting these;
-\code{\link[=build_hmm]{build_hmm()}} and \code{\link[=fit_model]{fit_model()}} for building and
-fitting hidden Markov models; \code{\link[=hidden_paths]{hidden_paths()}} for
-computing the most probable paths of hidden states; and \code{\link[=biofam3c]{biofam3c()}}
-\code{\link[=hmm_biofam]{hmm_biofam()}} for information on the data and model used in the example.
-}
diff --git a/man/stacked_sequence_plot.Rd b/man/stacked_sequence_plot.Rd
index 25ac48b6..ba12aea1 100644
--- a/man/stacked_sequence_plot.Rd
+++ b/man/stacked_sequence_plot.Rd
@@ -69,7 +69,9 @@ of length 1, or of length matching the number of channels to be plotted.}
 \code{\link[ggseqplot:ggseqdplot]{ggseqplot::ggseqdplot()}}.}
 }
 \description{
-Function \verb{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\verb{(e.g.,}hmm\code{and}mhmm`).
+Function \code{stacked_sequence_plot} draws stacked sequence plots of sequence
+object created with the \link[TraMineR:seqdef]{TraMineR::seqdef} function or observations and/or most
+probable paths of model objects of \code{seqHMM} (e.g., \code{hmm} and \code{mhmm}).
 }
 \examples{
 p <- stacked_sequence_plot(
diff --git a/src/EM.cpp b/src/EM.cpp
index dc692515..763bed6c 100644
--- a/src/EM.cpp
+++ b/src/EM.cpp
@@ -1,7 +1,7 @@
 // EM algorithm for non-mixture hidden Markov models
 
 #include "forward_backward.h"
-
+#include <omp.h>
 // [[Rcpp::export]]
 
 Rcpp::List EM(const arma::mat& transition_, const arma::cube& emission_, const arma::vec& init_,
diff --git a/src/EMx.cpp b/src/EMx.cpp
index 10cc48c9..deece6b7 100644
--- a/src/EMx.cpp
+++ b/src/EMx.cpp
@@ -2,6 +2,7 @@
 #include "forward_backward.h"
 #include "optcoef.h"
 #include "reparma.h"
+#include <omp.h>
 
 // [[Rcpp::export]]
 Rcpp::List EMx(const arma::mat& transition_, const arma::cube& emission_, const arma::vec& init_,
diff --git a/src/backward_nhmm.cpp b/src/backward_nhmm.cpp
index 5c6b95d1..1054280c 100644
--- a/src/backward_nhmm.cpp
+++ b/src/backward_nhmm.cpp
@@ -1,5 +1,4 @@
 // backward algorithm for NHMM
-#include <RcppArmadillo.h>
 #include "backward_nhmm.h"
 #include "get_parameters.h"
 #include "logsumexp.h"
diff --git a/src/forward_nhmm.cpp b/src/forward_nhmm.cpp
index 1b42369c..57448708 100644
--- a/src/forward_nhmm.cpp
+++ b/src/forward_nhmm.cpp
@@ -1,5 +1,4 @@
 // forward algorithm for NHMM
-#include <RcppArmadillo.h>
 #include "forward_nhmm.h"
 #include "get_parameters.h"
 #include "logsumexp.h"
diff --git a/src/forwardbackward.cpp b/src/forwardbackward.cpp
index 0af8ff7e..2c68080e 100644
--- a/src/forwardbackward.cpp
+++ b/src/forwardbackward.cpp
@@ -1,6 +1,6 @@
 // Forward-backward algorithm for non-mixture hidden Markov models
 #include "forward_backward.h"
-
+#include <omp.h>
 // [[Rcpp::export]]
 Rcpp::List forwardbackward(const arma::mat& transition, const arma::cube& emission,
   const arma::vec& init, const arma::ucube& obs, bool forwardonly, unsigned int threads) {
diff --git a/src/forwardbackwardx.cpp b/src/forwardbackwardx.cpp
index e747e242..98d30c49 100644
--- a/src/forwardbackwardx.cpp
+++ b/src/forwardbackwardx.cpp
@@ -1,9 +1,9 @@
 // Forward-backward algorithm for mixture hidden Markov models
 #include "forward_backward.h"
 #include "reparma.h"
+#include <omp.h>
 
 // [[Rcpp::export]]
-
 Rcpp::List forwardbackwardx(const arma::mat& transition, const arma::cube& emission,
   const arma::vec& init, const arma::ucube obs, const arma::mat& coef, const arma::mat& X,
   const arma::uvec& numberOfStates, bool forwardonly, unsigned int threads) {
diff --git a/src/internalBackward.cpp b/src/internalBackward.cpp
index 1189e0e2..82626057 100644
--- a/src/internalBackward.cpp
+++ b/src/internalBackward.cpp
@@ -1,5 +1,6 @@
 // Internal backward algorithms for HMMs and MHMMs
 #include "forward_backward.h"
+#include <omp.h>
 
 void internalBackward(const arma::mat& transition, const arma::cube& emission,
     const arma::ucube& obs, arma::cube& beta, const arma::mat& scales, 
diff --git a/src/internalForward.cpp b/src/internalForward.cpp
index f38087f1..95436aa8 100644
--- a/src/internalForward.cpp
+++ b/src/internalForward.cpp
@@ -1,6 +1,7 @@
 // Internal forward algorithms for HMMs and MHMMs
 
 #include "forward_backward.h"
+#include <omp.h>
 
 void internalForward(const arma::mat& transition_t, const arma::cube& emission, const arma::vec& init,
   const arma::ucube& obs, arma::cube& alpha, arma::mat& scales, unsigned int threads) {
diff --git a/src/logLikHMM.cpp b/src/logLikHMM.cpp
index fc22d59c..916cf0b8 100644
--- a/src/logLikHMM.cpp
+++ b/src/logLikHMM.cpp
@@ -1,5 +1,7 @@
 // log-likelihood of HMM
 #include <RcppArmadillo.h>
+#include <omp.h>
+
 // [[Rcpp::export]]
 Rcpp::NumericVector logLikHMM(const arma::mat& transition, const arma::cube& emission,
     const arma::vec& init, const arma::ucube& obs, unsigned int threads) {
diff --git a/src/logLikMixHMM.cpp b/src/logLikMixHMM.cpp
index aa0c3085..9784f280 100644
--- a/src/logLikMixHMM.cpp
+++ b/src/logLikMixHMM.cpp
@@ -1,9 +1,9 @@
 // log-likelihood of MHMM using log-space
 #include <RcppArmadillo.h>
 #include "reparma.h"
+#include <omp.h>
 
 // [[Rcpp::export]]
-
 Rcpp::NumericVector logLikMixHMM(const arma::mat& transition, const arma::cube& emission,
   const arma::vec& init, const arma::ucube& obs, const arma::mat& coef, const arma::mat& X,
   const arma::uvec& numberOfStates, unsigned int threads) {
diff --git a/src/log_EM.cpp b/src/log_EM.cpp
index e813fce3..db1c7873 100644
--- a/src/log_EM.cpp
+++ b/src/log_EM.cpp
@@ -2,8 +2,8 @@
 
 #include "log_forward_backward.h"
 #include "logsumexp.h"
+#include <omp.h>
 // [[Rcpp::export]]
-
 Rcpp::List log_EM(const arma::mat& transition_, const arma::cube& emission_, 
   const arma::vec& init_, const arma::ucube& obs, const arma::uvec& nSymbols, 
   int itermax, double tol, int trace, unsigned int threads) {
diff --git a/src/log_EMx.cpp b/src/log_EMx.cpp
index 225cbd91..6bee6d47 100644
--- a/src/log_EMx.cpp
+++ b/src/log_EMx.cpp
@@ -4,8 +4,9 @@
 #include "optcoef.h"
 #include "logsumexp.h"
 #include "reparma.h"
-// [[Rcpp::export]]
+#include <omp.h>
 
+// [[Rcpp::export]]
 Rcpp::List log_EMx(const arma::mat& transition_, const arma::cube& emission_, 
   const arma::vec& init_, const arma::ucube& obs, const arma::uvec& nSymbols, 
   const arma::mat& coef_, const arma::mat& X, const arma::uvec& numberOfStates, 
diff --git a/src/log_forwardbackward.cpp b/src/log_forwardbackward.cpp
index a007010c..67619f88 100644
--- a/src/log_forwardbackward.cpp
+++ b/src/log_forwardbackward.cpp
@@ -1,8 +1,9 @@
 // Forward-backward algorithm for non-mixture hidden Markov models using log-space
 
 #include "log_forward_backward.h"
-// [[Rcpp::export]]
+#include <omp.h>
 
+// [[Rcpp::export]]
 Rcpp::List log_forwardbackward(const arma::mat& transition_, const arma::cube& emission_, 
   const arma::vec& init_, const arma::ucube& obs, bool forwardonly, unsigned int threads) {
 
diff --git a/src/log_forwardbackwardx.cpp b/src/log_forwardbackwardx.cpp
index a1c04677..9bc97778 100644
--- a/src/log_forwardbackwardx.cpp
+++ b/src/log_forwardbackwardx.cpp
@@ -2,8 +2,9 @@
 
 #include "log_forward_backward.h"
 #include "reparma.h"
-// [[Rcpp::export]]
+#include <omp.h>
 
+// [[Rcpp::export]]
 Rcpp::List log_forwardbackwardx(const arma::mat& transition_, 
   const arma::cube& emission_, const arma::vec& init_,
   const arma::ucube& obs, const arma::mat& coef, const arma::mat& X,
diff --git a/src/log_internalBackward.cpp b/src/log_internalBackward.cpp
index 5e70b26e..6afc0768 100644
--- a/src/log_internalBackward.cpp
+++ b/src/log_internalBackward.cpp
@@ -2,6 +2,8 @@
 
 #include "forward_backward.h"
 #include "logsumexp.h"
+#include <omp.h>
+
 void log_internalBackward(const arma::mat& transition, const arma::cube& emission,
   const arma::ucube& obs, arma::cube& beta, unsigned int threads) {
 
diff --git a/src/log_internalForward.cpp b/src/log_internalForward.cpp
index c9977f6e..e319d206 100644
--- a/src/log_internalForward.cpp
+++ b/src/log_internalForward.cpp
@@ -2,6 +2,7 @@
 
 #include "forward_backward.h"
 #include "logsumexp.h"
+#include <omp.h>
 
 void log_internalForward(const arma::mat& transition, const arma::cube& emission,
   const arma::vec& init, const arma::ucube& obs, arma::cube& alpha, unsigned int threads) {
diff --git a/src/log_logLikHMM.cpp b/src/log_logLikHMM.cpp
index b90c9434..f143d56c 100644
--- a/src/log_logLikHMM.cpp
+++ b/src/log_logLikHMM.cpp
@@ -1,8 +1,9 @@
 // log-likelihood of HMM using log-space
 
 #include "logsumexp.h"
-// [[Rcpp::export]]
+#include <omp.h>
 
+// [[Rcpp::export]]
 Rcpp::NumericVector log_logLikHMM(const arma::mat& transition_, const arma::cube& emission_, 
   const arma::vec& init_, const arma::ucube& obs, unsigned int threads) {
   
diff --git a/src/log_logLikMixHMM.cpp b/src/log_logLikMixHMM.cpp
index 3309a0a4..a53f401b 100644
--- a/src/log_logLikMixHMM.cpp
+++ b/src/log_logLikMixHMM.cpp
@@ -2,8 +2,9 @@
 
 #include "logsumexp.h"
 #include "reparma.h"
-// [[Rcpp::export]]
+#include <omp.h>
 
+// [[Rcpp::export]]
 Rcpp::NumericVector log_logLikMixHMM(arma::mat transition, arma::cube emission, arma::vec init,
   const arma::ucube& obs, const arma::mat& coef, const arma::mat& X,
   const arma::uvec& numberOfStates, unsigned int threads) {
diff --git a/src/log_objective.cpp b/src/log_objective.cpp
index 29561e7b..b755e5eb 100644
--- a/src/log_objective.cpp
+++ b/src/log_objective.cpp
@@ -2,6 +2,7 @@
 
 #include "log_forward_backward.h"
 #include "logsumexp.h"
+#include <omp.h>
 
 // [[Rcpp::export]]
 Rcpp::List log_objective(const arma::mat& transition, const arma::cube& emission,
diff --git a/src/log_objectivex.cpp b/src/log_objectivex.cpp
index 4c5fec2b..731fa947 100644
--- a/src/log_objectivex.cpp
+++ b/src/log_objectivex.cpp
@@ -2,6 +2,8 @@
 #include "log_forward_backward.h"
 #include "logsumexp.h"
 #include "reparma.h"
+
+#include <omp.h>
 // [[Rcpp::export]]
 Rcpp::List log_objectivex(const arma::mat& transition, const arma::cube& emission,
   const arma::vec& init, const arma::ucube& obs, const arma::umat& ANZ,
diff --git a/src/objective.cpp b/src/objective.cpp
index b6d69e5e..225f8f90 100644
--- a/src/objective.cpp
+++ b/src/objective.cpp
@@ -1,6 +1,7 @@
 // log-likelihood and gradients of HMM
 #include "optcoef.h"
 #include "forward_backward.h"
+#include <omp.h>
 
 // [[Rcpp::export]]
 Rcpp::List objective(const arma::mat& transition, const arma::cube& emission,
diff --git a/src/objectivex.cpp b/src/objectivex.cpp
index 69aaf783..70eb2ea9 100644
--- a/src/objectivex.cpp
+++ b/src/objectivex.cpp
@@ -2,6 +2,8 @@
 #include "optcoef.h"
 #include "forward_backward.h"
 #include "reparma.h"
+#include <omp.h>
+
 // [[Rcpp::export]]
 Rcpp::List objectivex(const arma::mat& transition, const arma::cube& emission,
                 const arma::vec& init, const arma::ucube& obs, const arma::umat& ANZ,
diff --git a/src/optcoef.h b/src/optcoef.h
index c8078dcc..8f6809a2 100644
--- a/src/optcoef.h
+++ b/src/optcoef.h
@@ -1,7 +1,7 @@
 #ifndef OPTCOEF_H
 #define OPTCOEF_H
 
-
+#define ARMA_WARN_LEVEL 1
 #include <RcppArmadillo.h>
 unsigned int optCoef(arma::mat& weights, const arma::ucube& obs, const arma::cube& emission,
                      const arma::mat& bsi, arma::mat& coef, const arma::mat& X,
diff --git a/src/seqHMM.h b/src/seqHMM.h
deleted file mode 100644
index 9243d8eb..00000000
--- a/src/seqHMM.h
+++ /dev/null
@@ -1,14 +0,0 @@
-#ifndef SEQHMM_H
-#define SEQHMM_H
-
-#ifdef _OPENMP
-// [[Rcpp::plugins(openmp)]]
-#include <omp.h>
-#endif
-
-#define ARMA_NO_DEBUG
-#define ARMA_WARN_LEVEL 1
-#include <RcppArmadillo.h>
-// [[Rcpp::depends(RcppArmadillo)]]
-
-#endif