Arranging Plots from TraMineR and seqHMM

Question

I want to plot two sequence distribution plots with the TraMineR package (seqplot(type = 'd')) and below them two Hidden Markov Models which correspond to the sequence plots (those are clusters).

For the reprex, I just use the build in biofam data set. Using layout() produces the desired result up to the third plot. The fourth plot is drawn to a new page, so essentially the bottom right space is left empty, but the first three panes are used accordingly to the layout() specs.

library('seqHMM')
library('TraMineR')

data("biofam")
biofam <- biofam[sample(nrow(biofam),300),]
biofam.lab <- c("Parent", "Left", "Married", "Left+Marr",
                "Child", "Left+Child", "Left+Marr+Child", "Divorced")
biofam.seq <- seqdef(biofam, 10:25, labels=biofam.lab)
#>  [>] 8 distinct states appear in the data:
#>      1 = 0
#>      2 = 1
#>      3 = 2
#>      4 = 3
#>      5 = 4
#>      6 = 5
#>      7 = 6
#>      8 = 7
#>  [>] state coding:
#>        [alphabet]  [label]  [long label]
#>      1  0           0        Parent
#>      2  1           1        Left
#>      3  2           2        Married
#>      4  3           3        Left+Marr
#>      5  4           4        Child
#>      6  5           5        Left+Child
#>      7  6           6        Left+Marr+Child
#>      8  7           7        Divorced
#>  [>] 300 sequences in the data set
#>  [>] min/max sequence length: 16/16
### Plot the sequence
seqplot(biofam.seq,type='d',with.legend=FALSE)

### Fit a Hidden Markov Model
biofam.hmm <- build_hmm(observations = biofam.seq, n_states = 3)
### Plot the HMM once
plot(biofam.hmm, with.legend=FALSE)
### Use layout to plot two sequence distribution plots and below there HMM
layout(matrix(c(1,2,3,4), ncol = 2, byrow = TRUE))
layout.show(4)


par(cex = 0.65, mar = c(2.5,4.1,2.5,1.5))
### The Sequence plots
seqplot(biofam.seq, type = 'd', with.legend = FALSE, border = NA, ylab = NA)
seqplot(biofam.seq, type = 'd', with.legend = FALSE, border = NA, ylab = NA)
### Change the par() arguments
par(cex = 0.5, mar = c(1.5,3.5,3.5,1.5))
plot(biofam.hmm, with.legend = FALSE, vertex.size = 15)
plot(biofam.hmm, with.legend = FALSE, vertex.size = 15)

^{Created on 2022-11-20 by the reprex package (v2.0.1)}

As you can see, the layout.show() indicates the correct plotting dimensions, but the second HMM plot is drawn on the next pane. Is there anything I can do (besides from plotting one cluster after another with just the sequence distribution and the fitted HMM), to plot everything on one page?

Any help would be appreciated.

Kind regards

PS: Since I don't have enough reputation, I cannot include the rendered images. Hopefully you can recreate my issue with the given code by yourself.

Answer 1

The problem is that seqHMM:::plot.hmm―the method which is dispatched when calling plot on a "hmm" object―is messing around with the pars internally (bad practice). But we can modify the method by deactivating that behavior, and have a new plot_hmm function.

plot_hmm <- function(x, layout = "horizontal", pie = TRUE, vertex.size = 40, 
          vertex.label = "initial.probs", vertex.label.dist = "auto", 
          vertex.label.pos = "bottom", vertex.label.family = "sans", 
          loops = FALSE, edge.curved = TRUE, edge.label = "auto", edge.width = "auto", 
          cex.edge.width = 1, edge.arrow.size = 1.5, edge.label.family = "sans", 
          label.signif = 2, label.scientific = FALSE, label.max.length = 6, 
          trim = 1e-15, combine.slices = 0.05, combined.slice.color = "white", 
          combined.slice.label = "others", with.legend = "bottom", 
          ltext = NULL, legend.prop = 0.5, cex.legend = 1, ncol.legend = "auto", 
          cpal = "auto", cpal.legend = "auto", legend.order = TRUE, 
          main = NULL, withlegend, ...) {
  seqHMM:::check_deprecated_args(match.call())
  oldPar <- par(no.readonly = TRUE)
  if (is.null(main)) {
    par(mar = c(0.5, 0.5, 0.5, 0.5))
  }
  else {
    par(mai = c(0, 0, 1, 0))
  }
  # on.exit(par(oldPar), add = TRUE)
  # on.exit(par(mfrow = c(1, 1)), add = TRUE)
  dots <- list(...)
  labelprint <- function(z, labs) {
    if (labs == TRUE && (z > 0.001 || z == 0)) {
      labs <- FALSE
    }
    if (z < 10^-(label.max.length)) {
      z <- prettyNum(signif(round(z, digits = label.max.length), 
                            digits = label.signif), scientific = labs)
    }
    else {
      z <- prettyNum(signif(z, digits = label.signif), 
                     scientific = labs)
    }
  }
  if (!is.matrix(layout) && !is.function(layout)) {
    if (!(layout %in% c("horizontal", "vertical"))) {
      stop("Argument layout only accepts numerical matrices, igraph layout functions, or strings \"horizontal\" and \"vertical\".")
    }
  }
  if (!is.numeric(vertex.label.pos)) {
    choices <- c("bottom", "top", "left", "right")
    ind <- pmatch(vertex.label.pos, choices, duplicates.ok = TRUE)
    if (any(is.na(ind))) {
      stop("Argument vertex.label.pos only accepts values \"bottom\", \"top\", \"left\", \"right\" or a numerical vector.")
    }
    vertex.label.pos <- choices[ind]
  }
  choices <- c(TRUE, FALSE, "bottom", "top", "left", "right")
  ind <- pmatch(with.legend, choices)
  if (is.na(ind)) {
    stop("Argument with.legend must be one of TRUE, FALSE, \"bottom\", \"right\", \"top\", or \"left\".")
  }
  with.legend <- choices[ind]
  if (with.legend %in% c(TRUE, "auto")) {
    with.legend <- "bottom"
  }
  if (x$n_channels > 1) {
    x <- mc_to_sc(x)
  }
  if (pie == FALSE && with.legend != FALSE) {
    with.legend <- FALSE
  }
  if (!is.numeric(vertex.label.pos)) {
    vpos <- numeric(length(vertex.label.pos))
    for (i in 1:length(vertex.label.pos)) {
      if (vertex.label.pos[i] == "bottom") {
        vpos[i] <- pi/2
      }
      else if (vertex.label.pos[i] == "top") {
        vpos[i] <- -pi/2
      }
      else if (vertex.label.pos[i] == "left") {
        vpos[i] <- pi
      }
      else {
        vpos[i] <- 0
      }
    }
    vertex.label.pos <- vpos
  }
  if (length(vertex.label) == 1 && !is.na(vertex.label) && 
      vertex.label != FALSE) {
    if (vertex.label == "initial.probs") {
      vertex.label <- sapply(x$initial_probs, labelprint, 
                             labs = label.scientific)
    }
    else if (vertex.label == "names") {
      vertex.label <- x$state_names
    }
  }
  else if (length(vertex.label) != length(x$state_names)) {
    warning("The length of the vector provided for the argument \"vertex.label\" does not match the number of hidden states.")
    vertex.label <- rep(vertex.label, length.out = length(x$state_names))
  }
  if (is.character(vertex.label.dist)) {
    ind <- pmatch(vertex.label.dist, "auto")
    if (is.na(ind)) {
      stop("Argument vertex.label.dist only accepts the value \"auto\" or a numerical vector.")
    }
    vertex.label.dist <- vertex.size * 0.4/3.5
  }
  else if (length(vertex.label.dist) > 1 && length(vertex.label.dist) != 
           x$n_states) {
    warning("The length of the vector provided for the argument \"vertex.label.dist\" does not match the number of edges.")
    vertex.label.dist <- rep(vertex.label.dist, length.out = length(x$n_states))
  }
  transM <- x$transition_probs
  transM[transM < trim] <- 0
  edges <- transM
  edges[edges > 0] <- 1
  if (!loops) {
    diag(edges) <- 0
  }
  transitions <- transM
  if (loops == FALSE && length(transitions) > 1) {
    diag(transitions) <- 0
  }
  transitions <- t(transitions)[t(transitions) > 0]
  if (!is.na(edge.label) && edge.label != FALSE) {
    if (length(edge.label) == 1 && (edge.label == "auto" || 
                                    edge.label == TRUE)) {
      edge.label <- sapply(transitions, labelprint, labs = label.scientific)
    }
    else if (length(edge.label) > 1 && length(edge.label) != 
             length(transitions)) {
      warning("The length of the vector provided for the argument \"edge.label\" does not match the number of edges.")
      edge.label <- rep(edge.label, length.out = length(transitions))
    }
  }
  if (is.character(edge.width)) {
    ind <- pmatch(edge.width, "auto")
    if (is.na(ind)) {
      stop("Argument edge.width only accepts the value \"auto\" or a numerical vector.")
    }
    edge.width <- transitions * (7/max(transitions)) * cex.edge.width
  }
  else if (length(edge.width) > 1 && edge.width != length(transitions)) {
    warning("The length of the vector provided for the argument \"edge.width\" does not match the number of edges.")
    edge.width <- rep(edge.width, length.out = length(transitions))
  }
  g1 <- igraph::graph.adjacency(edges, mode = "directed")
  if (is.function(layout)) {
    glayout <- layout(g1)
  }
  else if (is.matrix(layout)) {
    glayout <- layout
  }
  else {
    if (layout == "horizontal") {
      glayout <- igraph::layout_on_grid(g1, width = x$n_states)
    }
    else if (layout == "vertical") {
      glayout <- igraph::layout_on_grid(g1, width = 1)
    }
  }
  if (length(cpal) == 1 && cpal == "auto") {
    pie.colors <- attr(x$observations, "cpal")
  }
  else if (length(cpal) != ncol(x$emiss)) {
    warning("The length of the vector provided for argument cpal does not match the number of observed states. Automatic color palette was used.")
    pie.colors <- attr(x$observations, "cpal")
  }
  else if (!all(isColor(cpal))) {
    stop(paste("Please provide a vector of colors for argument cpal or use value \"auto\" for automatic color palette."))
  }
  else {
    pie.colors <- cpal
  }
  if (with.legend != FALSE) {
    pie.colors.l <- pie.colors
  }
  if (with.legend != FALSE && pie == TRUE) {
    if (!is.null(ltext)) {
      if (legend.order) {
        if (length(ltext) != x$n_symbols) {
          stop(paste0("With legend.order = TRUE, the length of the argument ltext must match the number of (combined) observed states in the observed data (", 
                      x$n_symbols, ")."))
        }
      }
      else {
        if ((length(cpal) == 1 && cpal != "auto") && 
            length(ltext) != length(cpal.legend)) {
          stop(paste0("The number of colours in cpal.legend does not match the number of labels in ltext."))
        }
        ltext.orig <- ltext
        if (length(cpal) == 1 && cpal == "auto") {
          if (length(cpal.legend) == 1 && cpal.legend == 
              "auto") {
            cpal.legend <- attr(x$observations, "cpal")
          }
        }
        else {
          if (length(cpal.legend) == 1 && cpal.legend == 
              "auto") {
            cpal.legend <- cpal
          }
        }
      }
    }
    else {
      ltext <- ltext.orig <- x$symbol_names
      if (length(cpal) == 1 && cpal == "auto") {
        if (length(cpal.legend) == 1 && cpal.legend == 
            "auto") {
          cpal.legend <- attr(x$observations, "cpal")
        }
      }
      else {
        cpal.legend <- cpal
      }
    }
    if (with.legend == "bottom") {
      graphics::layout(matrix(1:2, nrow = 2), heights = c(1 - 
                                                            legend.prop, legend.prop))
    }
    else if (with.legend == "right") {
      graphics::layout(matrix(1:2, nrow = 1), widths = c(1 - 
                                                           legend.prop, legend.prop))
    }
    else if (with.legend == "left") {
      graphics::layout(matrix(2:1, nrow = 1), widths = c(legend.prop, 
                                                         1 - legend.prop))
    }
    else {
      graphics::layout(matrix(2:1, nrow = 2), widths = c(legend.prop, 
                                                         1 - legend.prop))
    }
    par(cex = 1)
  }
  if (!is.matrix(layout) && !is.function(layout)) {
    if (layout == "horizontal") {
      if (hasArg(rescale)) {
        rescale <- dots$rescale
      }
      else {
        rescale <- FALSE
      }
      if (hasArg(xlim)) {
        xlim <- dots$xlim
      }
      else {
        if (rescale == TRUE) {
          xlim <- c(-1, 1)
        }
        else {
          xlim <- c(-0.1, ncol(transM) - 1 + 0.1)
        }
      }
      if (hasArg(ylim)) {
        ylim <- dots$ylim
      }
      else {
        if (rescale == TRUE) {
          ylim <- c(-1, 1)
        }
        else {
          ylim <- c(-0.5, 0.5)
        }
      }
      dots[["xlim"]] <- NULL
      dots[["ylim"]] <- NULL
      dots[["rescale"]] <- NULL
    }
    else if (layout == "vertical") {
      if (hasArg(rescale)) {
        rescale <- dots$rescale
      }
      else {
        rescale <- FALSE
      }
      if (hasArg(xlim)) {
        xlim <- dots$xlim
      }
      else {
        if (rescale == TRUE) {
          xlim <- c(-1, 1)
        }
        else {
          xlim <- c(-0.5, 0.5)
        }
      }
      if (hasArg(ylim)) {
        ylim <- dots$ylim
      }
      else {
        if (rescale == TRUE) {
          ylim <- c(-1, 1)
        }
        else {
          ylim <- c(-0.1, ncol(transM) - 1 + 0.1)
        }
      }
      dots[["xlim"]] <- NULL
      dots[["ylim"]] <- NULL
      dots[["rescale"]] <- NULL
    }
  }
  if (pie == TRUE) {
    pie.values <- lapply(seq_len(nrow(transM)), function(i) x$emission_probs[i, 
    ])
    if (combine.slices > 0 && !all(unlist(pie.values)[unlist(pie.values) > 
                                                      0] > combine.slices)) {
      if (with.legend != FALSE) {
        pie.colors.l <- NULL
        lt <- NULL
      }
      for (i in 1:x$n_states) {
        cs.prob <- sum(pie.values[[i]][pie.values[[i]] < 
                                         combine.slices])
        pie.values[[i]][pie.values[[i]] < combine.slices] <- 0
        pie.values[[i]] <- c(pie.values[[i]], cs.prob)
        if (with.legend != FALSE) {
          pie.colors.l <- c(pie.colors.l, pie.colors[pie.values[[i]][1:(length(pie.values[[i]]) - 
                                                                          1)] >= combine.slices])
          lt <- c(lt, ltext[pie.values[[i]][1:(length(pie.values[[i]]) - 
                                                 1)] >= combine.slices])
        }
      }
      if (with.legend != FALSE) {
        ltext <- c(unique(lt), combined.slice.label)
        pie.colors.l <- c(unique(pie.colors.l), combined.slice.color)
      }
      if (ncol.legend == "auto") {
        if (with.legend == "bottom" || with.legend == 
            "top") {
          ncol.legend <- ceiling(length(pie.colors)/4)
        }
        else {
          ncol.legend <- 1
        }
      }
      pie.colors <- c(pie.colors, combined.slice.color)
    }
    else {
      if (ncol.legend == "auto") {
        if (with.legend == "bottom" || with.legend == 
            "top") {
          ncol.legend <- ceiling(ncol(x$emission_probs)/4)
        }
        else {
          ncol.legend <- 1
        }
      }
    }
    if (!is.matrix(layout) && !is.function(layout) && (layout == 
                                                       "horizontal" || layout == "vertical")) {
      do.call(igraph::plot.igraph, c(list(g1, layout = glayout, 
                                  vertex.shape = "pie", vertex.pie = pie.values, 
                                  vertex.pie.color = list(pie.colors), vertex.size = vertex.size, 
                                  vertex.label = vertex.label, vertex.label.dist = vertex.label.dist, 
                                  vertex.label.degree = vertex.label.pos, vertex.label.family = vertex.label.family, 
                                  edge.curved = edge.curved, edge.width = edge.width, 
                                  edge.label = edge.label, edge.label.family = edge.label.family, 
                                  edge.arrow.size = edge.arrow.size, xlim = xlim, 
                                  ylim = ylim, rescale = rescale, main = main), 
                             dots))
    }
    else {
      do.call(igraph::plot.igraph, c(list(g1, layout = glayout, 
                                  vertex.shape = "pie", vertex.pie = pie.values, 
                                  vertex.pie.color = list(pie.colors), vertex.size = vertex.size, 
                                  vertex.label = vertex.label, vertex.label.dist = vertex.label.dist, 
                                  vertex.label.degree = vertex.label.pos, vertex.label.family = vertex.label.family, 
                                  edge.curved = edge.curved, edge.width = edge.width, 
                                  edge.label = edge.label, edge.label.family = edge.label.family, 
                                  edge.arrow.size = edge.arrow.size, main = main), 
                             dots))
    }
  }
  else {
    if (!is.matrix(layout) && !is.function(layout) && (layout == 
                                                       "horizontal" || layout == "vertical")) {
      do.call(igraph::plot.igraph, c(list(g1, layout = glayout, 
                                  vertex.size = vertex.size, vertex.label = vertex.label, 
                                  vertex.label.dist = vertex.label.dist, vertex.label.degree = vertex.label.pos, 
                                  vertex.label.family = vertex.label.family, edge.curved = edge.curved, 
                                  edge.width = edge.width, edge.label = edge.label, 
                                  edge.label.family = edge.label.family, edge.arrow.size = edge.arrow.size, 
                                  xlim = xlim, ylim = ylim, rescale = rescale, 
                                  main = main), dots))
    }
    else {
      do.call(igraph::plot.igraph, c(list(g1, layout = glayout, 
                                  vertex.size = vertex.size, vertex.label = vertex.label, 
                                  vertex.label.dist = vertex.label.dist, vertex.label.degree = vertex.label.pos, 
                                  vertex.label.family = vertex.label.family, edge.curved = edge.curved, 
                                  edge.width = edge.width, edge.label = edge.label, 
                                  edge.label.family = edge.label.family, edge.arrow.size = edge.arrow.size, 
                                  main = main), dots))
    }
  }
  if (with.legend != FALSE && pie == TRUE) {
    if (legend.order) {
      seqlegend(x$observations, cpal = pie.colors.l, ltext = ltext, 
                position = "center", cex = cex.legend, ncol = ncol.legend, 
                with.missing = FALSE)
    }
    else {
      seqlegend(x$observations, cpal = cpal.legend, ltext = ltext.orig, 
                position = "center", cex = cex.legend, ncol = ncol.legend, 
                with.missing = FALSE)
    }
  }
  # par(mfrow = c(1, 1))
}

Now it works.

layout(matrix(c(1, 2, 3, 4), ncol=2, byrow=TRUE))
# layout.show(4)

par(cex=0.65, mar=c(2.5, 4.1, 2.5, 1.5))
seqplot(biofam.seq, type='d', with.legend=FALSE, bordplot_hmmer=NA, ylab=NA)
seqplot(biofam.seq, type='d', with.legend=FALSE, border=NA, ylab=NA)
par(cex=0.5, mar=c(1.5, 3.5, 3.5, 1.5))
plot_hmm(biofam.hmm, with.legend=FALSE, vertex.size=15)
plot_hmm(biofam.hmm, with.legend=FALSE, vertex.size=15)