R/TS_plots.R

In weecology/LDATS: Latent Dirichlet Allocation Coupled with Time Series Analyses

#' @title Plot an LDATS TS model
#'
#' @description Generalization of the \code{\link[graphics]{plot}} function to 
#'   work on fitted TS model objects (class \code{TS_fit}) returned from 
#'   \code{\link{TS}}. 
#' 
#' @param x A \code{TS_fit} object of a multinomial time series model fit by
#'   \code{\link{TS}}.
#' 
#' @param ... Additional arguments to be passed to subfunctions. Not currently
#'   used, just retained for alignment with \code{\link[graphics]{plot}}.
#'
#' @param plot_type "diagnostic" or "summary".
#'
#' @param bin_width Width of the bins used in the histograms of the summary 
#'   time series plot, in units of the x-axis (the time variable used to fit 
#'   the model).
#'
#' @param xname Label for the x-axis in the summary time series plot. Defaults
#'   to \code{NULL}, which results in usage of the \code{timename} element
#'   of the control list (held in\code{control$TS_control$timename}). To have
#'   no label printed, set \code{xname = ""}.
#'
#' @param border Border for the histogram, default is \code{NA}.
#'
#' @param selection Indicator of the change points to use in the time series
#'   summary plot. Currently only defined for \code{"median"} and 
#'   \code{"mode"}.
#'
#' @param cols \code{list} of elements used to define the colors for the two
#'   panels of the summary plot, as generated simply using 
#'   \code{\link{set_TS_summary_plot_cols}}. \code{cols} has two elements 
#'   \code{rho} and \code{gamma}, each corresponding to the related panel, 
#'   and each containing default values for entries named \code{cols},
#'   \code{option}, and \code{alpha}. See \code{\link{set_gamma_colors}} and 
#'   \code{\link{set_rho_hist_colors}} for details on usage.
#'
#' @param LDATS \code{logical} indicating if the plot is part of a larger 
#'   LDATS plot output.
#'
#' @param interactive \code{logical} input, should be code{TRUE} unless
#'   testing.
#' 
#' @return \code{NULL}.
#' 
#' @examples 
#' \donttest{
#'   data(rodents)
#'   document_term_table <- rodents$document_term_table
#'   document_covariate_table <- rodents$document_covariate_table
#'   LDA_models <- LDA_set(document_term_table, topics = 2)[[1]]
#'   data <- document_covariate_table
#'   data$gamma <- LDA_models@gamma
#'   weights <- document_weights(document_term_table)
#'   TSmod <- TS(data, gamma ~ 1, nchangepoints = 1, "newmoon", weights)
#'   plot(TSmod)
#' }
#'
#' @export 
#'
plot.TS_fit <- function(x, ..., plot_type = "summary", interactive = FALSE,
                        cols = set_TS_summary_plot_cols(), bin_width = 1, 
                        xname = NULL, border = NA, selection = "median", 
                        LDATS = FALSE){
  if (plot_type == "diagnostic"){
    TS_diagnostics_plot(x, interactive = interactive)
  } else if (plot_type == "summary"){
    TS_summary_plot(x, cols, bin_width, xname, border, selection, LDATS)
  }
}

#' @title Plot the diagnostics of the parameters fit in a TS model
#'
#' @description Plot 4-panel figures (showing trace plots, posterior ECDF, 
#'   posterior density, and iteration autocorrelation) for each of the 
#'   parameters (change point locations and regressors) fitted within a 
#'   multinomial time series model (fit by \code{\link{TS}}). \cr \cr
#'   \code{eta_diagnostics_plots} creates the diagnostic plots
#'   for the regressors (etas) of a time series model. \cr \cr
#'   \code{rho_diagnostics_plots} creates the diagnostic plots
#'   for the change point locations (rho) of a time series model.
#'
#' @param interactive \code{logical} input, should be code{TRUE} unless
#'   testing.
#'
#' @param x Object of class \code{TS_fit}, generated by \code{\link{TS}} to
#'   have its diagnostics plotted.
#' 
#' @return \code{NULL}.
#' 
#' @examples 
#' \donttest{
#'   data(rodents)
#'   document_term_table <- rodents$document_term_table
#'   document_covariate_table <- rodents$document_covariate_table
#'   LDA_models <- LDA_set(document_term_table, topics = 2)[[1]]
#'   data <- document_covariate_table
#'   data$gamma <- LDA_models@gamma
#'   weights <- document_weights(document_term_table)
#'   TSmod <- TS(data, gamma ~ 1, nchangepoints = 1, "newmoon", weights)
#'   TS_diagnostics_plot(TSmod)
#' }
#'
#' @export 
#'
TS_diagnostics_plot <- function(x, interactive = TRUE){
  rho_diagnostics_plots(x, interactive)
  eta_diagnostics_plots(x, interactive)
}

#' @rdname TS_diagnostics_plot
#'
#' @export 
#'
eta_diagnostics_plots <- function(x, interactive){
  on.exit(devAskNewPage(FALSE))
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  etas <- x$etas
  if (is.null(etas)){
    return()
  }
  netas <- ncol(etas) 
  for (i in 1:netas){
    devAskNewPage(interactive)
    par(mfrow = c(2, 2), mar = c(5, 5, 1, 1))
    chead <- colnames(etas)[i]
    spl1 <- strsplit(chead, "_")[[1]]
    seglab <- paste0("Segment ", spl1[1])
    spl2 <- strsplit(spl1[2], ":")[[1]]
    toplab <- paste0(" Topic ", spl2[1])
    coflab <- gsub("\\(Intercept\\)", "Intercept", spl2[2])
    lab <- paste0(seglab, toplab, " ", coflab)
    trace_plot(etas[ , i], lab)
    ecdf_plot(etas[ , i], lab)
    posterior_plot(etas[ , i], lab)
    autocorr_plot(etas[ , i])
  }  
}

#' @rdname TS_diagnostics_plot
#'
#' @export 
#'
rho_diagnostics_plots <- function(x, interactive){
  on.exit(devAskNewPage(FALSE))
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  rhos <- x$rhos
  if (is.null(rhos)){
    return()
  }
  nrhos <- ncol(rhos) 
  for (i in 1:nrhos){
    devAskNewPage(interactive)
    par(mfrow = c(2, 2), mar = c(5, 5, 1, 1))
    lab <- paste0("Change point ", i, " location")
    trace_plot(rhos[ , i], lab)
    ecdf_plot(rhos[ , i], lab)
    posterior_plot(rhos[ , i], lab)
    autocorr_plot(rhos[ , i])
  }
}

#' @title Produce the trace plot panel for the TS diagnostic plot of a 
#'   parameter
#' 
#' @description Produce a trace plot for the parameter of interest (rho or 
#'   eta) as part of \code{\link{TS_diagnostics_plot}}. A horizontal line 
#'   is added to show the median of the posterior.
#'
#' @param x Vector of parameter values drawn from the posterior distribution,
#'   indexed to the iteration by the order of the vector.
#'
#' @param ylab \code{character} value used to label the y axis.
#' 
#' @return \code{NULL}.
#'
#' @examples
#'  trace_plot(rnorm(100, 0, 1))
#'
#' @export
#'
trace_plot <- function(x, ylab = "parameter value"){
  plot(x, type = "l", lwd = 1, col = 0,
       xlab = "Iteration", ylab = ylab, las = 1, bty = "L")
  ext <- 0.01 * length(x)
  points(c(-ext, length(x) + ext), rep(median(x), 2), type = "l", lwd = 2, 
        lty = 2)
  points(seq_along(x), x, type = "l", col = rgb(0.4, 0.4, 0.4, alpha = 0.9))
}

#' @title Produce the posterior distribution ECDF panel for the TS 
#'   diagnostic plot of a parameter
#' 
#' @description Produce a vanilla ECDF (empirical cumulative distribution
#'   function) plot using \code{ecdf} for the parameter of interest (rho or 
#'   eta) as part of \code{\link{TS_diagnostics_plot}}. A horizontal line 
#'   is added to show the median of the posterior.
#'
#' @param x Vector of parameter values drawn from the posterior distribution,
#'   indexed to the iteration by the order of the vector.
#'
#' @param xlab \code{character} value used to label the x axis.
#' 
#' @return \code{NULL}.
#'
#' @examples
#'  ecdf_plot(rnorm(100, 0, 1))
#'
#' @export
#'
ecdf_plot <- function(x, xlab = "parameter value"){
  ECDF <- ecdf(x)
  plot(ECDF, main = "", xlab = xlab, ylab = "%", las = 1, bty = "L")
  abline(a = 0.5, b = 0, lwd = 2, lty = 2)
}

#' @title Produce the posterior distribution histogram panel for the TS 
#'   diagnostic plot of a parameter
#' 
#' @description Produce a vanilla histogram plot using \code{hist} for the 
#'   parameter of interest (rho or eta) as part of 
#'   \code{\link{TS_diagnostics_plot}}. A vertical line is added to show the 
#'   median of the posterior.
#'
#' @param x Vector of parameter values drawn from the posterior distribution,
#'   indexed to the iteration by the order of the vector.
#'
#' @param xlab \code{character} value used to label the x axis.
#' 
#' @return \code{NULL}.
#'
#' @examples
#'  posterior_plot(rnorm(100, 0, 1))
#'
#' @export
#'
posterior_plot <- function(x, xlab = "parameter value"){
  hist(x, las = 1, main = "", xlab = xlab)
  points(rep(median(x), 2), c(0, 1e5), type = "l", lwd = 2, lty = 2)
}

#' @title Produce the autocorrelation panel for the TS diagnostic plot of
#'   a parameter
#' 
#' @description Produce a vanilla ACF plot using \code{\link[stats]{acf}} for
#'   the parameter of interest (rho or eta) as part of 
#'   \code{\link{TS_diagnostics_plot}}.
#'
#' @param x Vector of parameter values drawn from the posterior distribution,
#'   indexed to the iteration by the order of the vector.
#' 
#' @return \code{NULL}.
#'
#' @examples
#'  autocorr_plot(rnorm(100, 0, 1))
#'
#' @export
#'
autocorr_plot <- function(x){
  acf(x, las = 1, ylab = "Autocorrelation")
}

#' @title Create the list of colors for the TS summary plot
#'
#' @description A default list generator function that produces the options
#'   for the colors controlling the panels of the TS summary plots, so needed
#'   because the panels should be in different color schemes. See 
#'   \code{\link{set_gamma_colors}} and \code{\link{set_rho_hist_colors}} for
#'   specific details on usage.
#'
#' @param rho_cols Colors to be used to plot the histograms of change points.
#'   Any valid color values (\emph{e.g.}, see \code{\link[grDevices]{colors}},
#'   \code{\link[grDevices]{rgb}}) can be input as with a standard plot. 
#'   The default (\code{rho_cols = NULL}) triggers use of 
#'   \code{\link[viridis]{viridis}} color options (see \code{rho_option}).
#'
#' @param rho_option A \code{character} string indicating the color option
#'   from \code{\link[viridis]{viridis}} to use if `rho_cols == NULL`. Four 
#'   options are available: "magma" (or "A"), "inferno" (or "B"), "plasma" 
#'   (or "C"), "viridis" (or "D", the default option) and "cividis" (or "E").
#'
#' @param rho_alpha Numeric value [0,1] that indicates the transparency of the 
#'   colors used. Supported only on some devices, see 
#'   \code{\link[grDevices]{rgb}}.
#'
#' @param gamma_cols Colors to be used to plot the LDA topic proportions,
#'   time series of observed topic proportions, and time series of fitted 
#'   topic proportions. Any valid color values (\emph{e.g.}, see 
#'   \code{\link[grDevices]{colors}}, \code{\link[grDevices]{rgb}}) can be 
#'   input as with a standard plot. The default (\code{gamma_cols = NULL})
#'   triggers use of \code{\link[viridis]{viridis}} color options (see 
#'   \code{gamma_option}).
#'
#' @param gamma_option A \code{character} string indicating the color option
#'   from \code{\link[viridis]{viridis}} to use if gamma_cols == NULL`. Four 
#'   options are available: "magma" (or "A"), "inferno" (or "B"), "plasma" 
#'   (or "C"), "viridis" (or "D", the default option) and "cividis" (or "E").
#'
#' @param gamma_alpha Numeric value [0,1] that indicates the transparency of 
#'   the colors used. Supported only on some devices, see 
#'   \code{\link[grDevices]{rgb}}.
#'
#' @return \code{list} of elements used to define the colors for the two
#'   panels. Contains two elements \code{rho} and \code{gamma}, each 
#'   corresponding to the related panel, and each containing default values 
#'   for entries named \code{cols}, \code{option}, and \code{alpha}. 
#'
#' @examples
#'   set_TS_summary_plot_cols()
#'
#' @export
#'
set_TS_summary_plot_cols <- function(rho_cols = NULL, rho_option = "D", 
                                     rho_alpha = 0.4, gamma_cols = NULL, 
                                     gamma_option = "C", gamma_alpha = 0.8){
  list(
    rho = list(cols = rho_cols, option = rho_option, alpha = rho_alpha),
    gamma = list(cols = gamma_cols, option = gamma_option, 
                 alpha = gamma_alpha)
  )
}

#' @title Create the summary plot for a TS fit to an LDA model
#'
#' @description Produces a two-panel figure of [1] the change point 
#'   distributions as histograms over time and [2] the time series of the 
#'   fitted topic proportions over time, based on a selected set of 
#'   change point locations. \cr \cr
#'   \code{pred_gamma_TS_plot} produces a time series of the 
#'   fitted topic proportions over time, based on a selected set of change 
#'   point locations. \cr \cr
#'   \code{rho_hist}: make a plot of the change point 
#'   distributions as histograms over time.
#'
#' @param x Object of class \code{TS_fit} produced by \code{\link{TS}}.
#'
#' @param cols \code{list} of elements used to define the colors for the two
#'   panels, as generated simply using \code{\link{set_TS_summary_plot_cols}}. 
#'   Has two elements \code{rho} and \code{gamma}, each corresponding to the
#'   related panel, and each containing default values for entries named
#'   \code{cols}, \code{option}, and \code{alpha}. See
#'   \code{\link{set_gamma_colors}} and \code{\link{set_rho_hist_colors}} for
#'   details on usage.
#'
#' @param bin_width Width of the bins used in the histograms, in units of the
#'   x-axis (the time variable used to fit the model).
#'
#' @param xname Label for the x-axis in the summary time series plot. Defaults
#'   to \code{NULL}, which results in usage of the \code{timename} element
#'   of the control list (held in\code{control$TS_control$timename}). To have
#'   no label printed, set \code{xname = ""}.
#'
#' @param border Border for the histogram, default is \code{NA}.
#'
#' @param selection Indicator of the change points to use. Currently only
#'   defined for "median" and "mode".
#'
#' @param together \code{logical} indicating if the subplots are part of a 
#'   larger LDA plot output. 
#'
#' @param LDATS \code{logical} indicating if the plot is part of a larger 
#'   LDATS plot output.
#' 
#' @return \code{NULL}.
#' 
#' @examples 
#' \donttest{
#'   data(rodents)
#'   document_term_table <- rodents$document_term_table
#'   document_covariate_table <- rodents$document_covariate_table
#'   LDA_models <- LDA_set(document_term_table, topics = 2)[[1]]
#'   data <- document_covariate_table
#'   data$gamma <- LDA_models@gamma
#'   weights <- document_weights(document_term_table)
#'   TSmod <- TS(data, gamma ~ 1, nchangepoints = 1, "newmoon", weights)
#'   TS_summary_plot(TSmod)
#'   pred_gamma_TS_plot(TSmod)
#'   rho_hist(TSmod)
#' }
#'
#' @export
#'
TS_summary_plot <- function(x, cols = set_TS_summary_plot_cols(), 
                            bin_width = 1, xname = NULL, border = NA, 
                            selection = "median", LDATS = FALSE){
  rc <- cols$rho
  rho_cols <- set_rho_hist_colors(x$rhos, rc$cols, rc$option, rc$alpha)
  rho_hist(x, rho_cols, bin_width, xname, border, TRUE, LDATS)
  gc <- cols$gamma
  gamma_cols <- set_gamma_colors(x, gc$cols, gc$option, gc$alpha)
  pred_gamma_TS_plot(x, selection, gamma_cols, xname, TRUE, LDATS)
}

#' @rdname TS_summary_plot
#'
#' @export
#'
pred_gamma_TS_plot <- function(x, selection = "median", 
                               cols = set_gamma_colors(x),
                               xname = NULL, together = FALSE, LDATS = FALSE){
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  if(LDATS){
    par(fig = c(0, 1, 0, 0.3), new = TRUE)
  } else if(together){
    par(fig = c(0, 1, 0, 0.52), new = TRUE)
  } else{
    par(fig = c(0, 1, 0, 1))
  }    
  rhos <- x$rhos
  nrhos <- ncol(rhos)
  if (!is.null(nrhos)){
    if (selection == "median"){
      spec_rhos <- apply(rhos, 2, median)
    } else if (selection == "mode"){
      spec_rhos <- apply(rhos, 2, modalvalue)
    } else {
      stop("selection input not supported")
    }
  } else{
    spec_rhos <- NULL
  }
  x$control$timename <- NULL # to remove from v0.1.0 model fits
  seg_mods <- multinom_TS(x$data, x$formula, spec_rhos,  
                          x$timename, x$weights, x$control)
  nsegs <- length(seg_mods[[1]])
  t1 <- min(x$data[, x$timename])
  t2 <- max(x$data[, x$timename])

  if (is.null(xname)){
    xname <- x$timename
  }
  par(mar = c(4, 5, 1, 1))
  plot(1, 1, type = "n", bty = "L", xlab = "", ylab = "", xaxt = "n", 
       yaxt = "n", ylim = c(0, 1), xlim = c(t1 - 1, t2 + 1))
  yax <- round(seq(0, 1, length.out = 5), 3)
  axis(2, at = yax, las = 1)
  axis(1)
  mtext(side = 2, line = 3.5, cex = 1.25, "Proportion")
  mtext(side = 1, line = 2.5, cex = 1.25, xname)
  ntopics <- ncol(as.matrix(x$data[[x$control$response]]))
  seg1 <- c(0, spec_rhos[-length(rhos)])
  seg2 <- c(spec_rhos, t2)
  time_obs <- rep(NA, nrow(x$data))
  pred_vals <- matrix(NA, nrow(x$data), ntopics)
  sp1 <- 1
  for (i in 1:nsegs){
    mod_i <- seg_mods[[1]][[i]]
    spec_vals <- sp1:(sp1 + nrow(mod_i$fitted.values) - 1)
    pred_vals[spec_vals, ] <- mod_i$fitted.values
    time_obs[spec_vals] <- mod_i$timevals
    sp1 <- sp1 + nrow(mod_i$fitted.values)
  }
  for (i in 1:ntopics){
    points(time_obs, pred_vals[ , i], type = "l", lwd = 3, col = cols[i])
  }
  if(!is.null(spec_rhos)){
    rho_lines(spec_rhos)
  }
}

#' @title Add change point location lines to the time series plot
#'
#' @description Adds vertical lines to the plot of the time series of fitted
#'   proportions associated with the change points of interest.
#'
#' @param spec_rhos \code{numeric} vector indicating the locations along the
#'   x axis where the specific change points being used are located.
#' 
#' @examples 
#' \donttest{
#'   data(rodents)
#'   document_term_table <- rodents$document_term_table
#'   document_covariate_table <- rodents$document_covariate_table
#'   LDA_models <- LDA_set(document_term_table, topics = 2)[[1]]
#'   data <- document_covariate_table
#'   data$gamma <- LDA_models@gamma
#'   weights <- document_weights(document_term_table)
#'   TSmod <- TS(data, gamma ~ 1, nchangepoints = 1, "newmoon", weights)
#'   pred_gamma_TS_plot(TSmod)
#'   rho_lines(200)
#' }
#'
#' @export
#'
rho_lines <- function(spec_rhos) {
  if(is.null(spec_rhos)) {
    return()
  }
  for (spec_rho in spec_rhos) {
    points(rep(spec_rho, 2), c(0, 1), type = "l", lwd = 4, 
           col = rgb(0.3, 0.3, 0.3, 0.6))
  }
}

#' @rdname TS_summary_plot
#'
#' @export
#'
rho_hist <- function(x, cols = set_rho_hist_colors(x$rhos), bin_width = 1, 
                     xname = NULL, border = NA, together = FALSE, 
                     LDATS = FALSE){
  oldpar <- par(no.readonly = TRUE)
  on.exit(par(oldpar))
  if(LDATS){
    par(fig = c(0, 1, 0.3, 0.55), mar = c(1.5, 5, 1, 1), new = TRUE)
  } else if(together){
    par(fig = c(0, 1, 0.54, 1), mar = c(1.5, 5, 1, 1))
  } else{
    par(fig = c(0, 1, 0, 1), mar = c(4, 5, 1, 1))
  } 
  rhos <- x$rhos
  nrhos <- ncol(rhos)
  niter <- nrow(rhos)
  timeobs <- x$data[, x$timename]
  timerange <- range(timeobs)
  timevals <- seq(timerange[1], timerange[2], 1)
  ntimes <- length(timevals) 
  
  binsteps <- seq(timerange[1], timerange[2] + bin_width, bin_width)
  bin1 <- binsteps[1:(length(binsteps)-1)]
  bin2 <- binsteps[2:(length(binsteps))] 
  nbins <- length(bin1)

  maxht <- 1
  if (!is.null(nrhos)){
    hts <- matrix(NA, nbins, nrhos)
    for(i in 1:nbins){
      for(j in 1:nrhos){
        hts[i,j]<- length(which(rhos[,j] >= bin1[i] & rhos[,j] < bin2[i]))
      }
    }
 
    maxht <- max(hts) / niter
    maxht <- ceiling(maxht * 100) / 100
    if (is.null(xname) & !together & !LDATS){
      xname <- x$timename
    }
  }
  plot(1, 1, type = "n", bty = "L", xlab = xname, ylab = "", xaxt = "n", 
       yaxt = "n", ylim = c(0, maxht), xlim = c(bin1[1], bin2[nbins]))
  yax <- round(seq(0, maxht, length.out = 5), 3)
  axis(2, at = yax, las = 1)
  axis(1)
  mtext(side = 2, line = 3.75, cex = 1.25, "Proportion")
  if (!is.null(nrhos)){
    for (i in 1:nbins){
      rho_ord <- order(hts[i,], decreasing = TRUE)
      for(j in 1:nrhos){
        ht_j <- hts[i, rho_ord[j]] / niter
        col_j <- cols[rho_ord[j]]
        rect(bin1[i], 0, bin2[i], ht_j, col = col_j, border = border)      
      }
    }
  }
}

#' @title Prepare the colors to be used in the change point histogram
#'
#' @description Based on the inputs, create the set of colors to be used in
#'   the change point histogram.
#' 
#' @param x \code{matrix} of change point locations (element \code{rhos}) 
#'   from an object of class \code{TS_fit}, fit by \code{\link{TS}}.
#'
#' @param cols Colors to be used to plot the histograms of change points.
#'   Any valid color values (\emph{e.g.}, see \code{\link[grDevices]{colors}},
#'   \code{\link[grDevices]{rgb}}) can be input as with a standard plot. 
#'   The default (\code{rho_cols = NULL}) triggers use of 
#'   \code{\link[viridis]{viridis}} color options (see \code{rho_option}).
#' 
#' @param option A \code{character} string indicating the color option
#'   from \code{\link[viridis]{viridis}} to use if "cols == NULL". Four 
#'   options are available: "magma" (or "A"), "inferno" (or "B"), "plasma" 
#'   (or "C"), "viridis" (or "D", the default option) and "cividis" (or "E").
#'
#' @param alpha Numeric value [0,1] that indicates the transparency of the 
#'   colors used. Supported only on some devices, see 
#'   \code{\link[grDevices]{rgb}}.
#'
#' @return Vector of \code{character} hex codes indicating colors to use.
#'
#' 
#' @examples 
#' \donttest{
#'   data(rodents)
#'   document_term_table <- rodents$document_term_table
#'   document_covariate_table <- rodents$document_covariate_table
#'   LDA_models <- LDA_set(document_term_table, topics = 2)[[1]]
#'   data <- document_covariate_table
#'   data$gamma <- LDA_models@gamma
#'   weights <- document_weights(document_term_table)
#'   TSmod <- TS(data, gamma ~ 1, nchangepoints = 1, "newmoon", weights)
#'   set_rho_hist_colors(TSmod$rhos)
#' }
#'
#' @export 
#'
set_rho_hist_colors <- function(x = NULL, cols = NULL, option = "D", 
                                alpha = 1){
  if(is.null(x)){
    return(NULL)
  }

  nrhos <- ncol(x)
  if (length(cols) == 0){
    cols <- viridis(nrhos, option = option, alpha = alpha, end = 0.9)
  }
  if (length(cols) == 1){
    if (cols == "greys" | cols == "grey" | cols == "grays" | cols == "gray"){
      ggg <- seq(0, 0.8, length.out = nrhos)
      cols <- rep(NA, nrhos)
      for (i in 1:nrhos){
       cols[i] <- rgb(ggg[i], ggg[i], ggg[i], alpha = alpha)
      }
    }
  }
  if (length(cols) > nrhos){
    cols <- cols[1:nrhos]
  }
  if (length(cols) < nrhos){
    nc <- length(cols)
    nt <- nrhos
    msg <- paste0("Fewer colors (", nc, ") provided than change points (", 
                  nt, ")")
    stop(msg)
  }
  cols
}

#' @title Prepare the colors to be used in the gamma time series
#'
#' @description Based on the inputs, create the set of colors to be used in
#'   the time series of the fitted gamma (topic proportion) values.
#' 
#' @param x Object of class \code{TS_fit}, fit by \code{\link{TS}}.
#'
#' @param cols Colors to be used to plot the time series of fitted topic 
#'   proportions.
#' 
#' @param option A \code{character} string indicating the color option
#'   from \code{\link[viridis]{viridis}} to use if "cols == NULL". Four 
#'   options are available: "magma" (or "A"), "inferno" (or "B"), "plasma" 
#'   (or "C"), "viridis" (or "D", the default option) and "cividis" (or "E").
#'
#' @param alpha Numeric value [0,1] that indicates the transparency of the 
#'   colors used. Supported only on some devices, see 
#'   \code{\link[grDevices]{rgb}}.
#'
#' @return Vector of \code{character} hex codes indicating colors to use.
#' 
#' @examples 
#' \donttest{
#'   data(rodents)
#'   document_term_table <- rodents$document_term_table
#'   document_covariate_table <- rodents$document_covariate_table
#'   LDA_models <- LDA_set(document_term_table, topics = 2)[[1]]
#'   data <- document_covariate_table
#'   data$gamma <- LDA_models@gamma
#'   weights <- document_weights(document_term_table)
#'   TSmod <- TS(data, gamma ~ 1, nchangepoints = 1, "newmoon", weights)
#'   set_gamma_colors(TSmod)
#' }
#'
#' @export 
#'
set_gamma_colors <- function(x, cols = NULL, option = "D", alpha = 1){
  if(is.null(x)){
    return(NULL)
  }

  ntopics <- ncol(as.matrix(x$data[x$control$response]))
  if (length(cols) == 0){
    cols <- viridis(ntopics, option = option, alpha = alpha, end = 0.9)
  }
  if (length(cols) == 1){
    if (cols == "greys" | cols == "grey" | cols == "grays" | cols == "gray"){
      ggg <- seq(0, 0.8, length.out = ntopics)
      cols <- rep(NA, ntopics)
      for (i in 1:ntopics){
       cols[i] <- rgb(ggg[i], ggg[i], ggg[i], alpha = alpha)
      }
    }
  }
  if (length(cols) > ntopics){
    cols <- cols[1:ntopics]
  }
  if (length(cols) < ntopics){
    nc <- length(cols)
    nt <- ntopics
    msg <- paste0("Fewer colors (", nc, ") provided than number of topics (", 
                  nt, ")")
    stop(msg)
  }
  cols
}