R/drawpostanalysis.R

In NetworkChange: Bayesian Package for Network Changepoint Analysis

#' Plot of latent node cluster
#'
#' Plot latent node cluster using ggplot2.
#' Uses tidyr for data reshaping and viridis for colorblind-friendly palettes.
#'
#' @param mcmcout NetworkChange output
#' @param Y Input raw data
#' @param point.cex node point size.  Default is 3.
#' @param text.cex node label size.  Default is 3.
#' @param segment.size segment size.  Default is 0.1.
#' @param n.cluster number of cluster. Default is 3.
#' @param start start of ts object
#' @param frequency frequency of ts object
#'
#' @references   Jong Hee Park and Yunkyun Sohn. 2020. "Detecting Structural Change
#' in Longitudinal Network Data." \emph{Bayesian Analysis}. Vol.15, No.1, pp.133-157.
#'
#' @return A plot object
#'
#' @importFrom tidyr pivot_longer
#' @importFrom ggrepel geom_text_repel
#' @importFrom ggplot2 scale_color_viridis_d alpha
#' @importFrom rlang .data
#'
#' @export
#'
#' @examples
#'
#'    \dontrun{
#'    set.seed(1973)
#'    ## generate an array with two constant blocks
#'    data(MajorAlly)
#'    Y <- MajorAlly
#'    fit <- NetworkChange(newY, R=2, m=2, mcmc=G, initial.s = initial.s,
#'           burnin=G, verbose=0, v0=v0, v1=v1)
#'    drawPostAnalysis(fit, Y, n.cluster=c(4, 4, 3))
#'    }

drawPostAnalysis <- function(mcmcout, Y, point.cex=3,  text.cex=3,
                             segment.size = 0.1, n.cluster = NULL,
                             start=1, frequency=1){
  m <- attr(mcmcout, "m")
  mcmc <- attr(mcmcout, "mcmc")
  U <- attr(mcmcout, "U")
  V <- attr(mcmcout, "V")
  R <- attr(mcmcout, "R")
  K <- dim(Y)
  T <- K[3]
  if(is.null(dimnames(Y))){
     dimnames(Y)[[1]] <- as.character(1:dim(Y)[1])
     dimnames(Y)[[2]] <- as.character(1:dim(Y)[2])
     dimnames(Y)[[3]] <- as.character(1:dim(Y)[3])
  }
  median.s <- ceiling(apply(attr(mcmcout, "Smat"), 2, median))
  y <- ts(1:T, start=start, frequency=frequency)
  ns <- m + 1
  time <- K[[3]]
  Year <- dimnames(Y)[[3]]

  if(is.null(n.cluster)){
    n.cluster <- rep(3, ns)
  }
  ## regime specific country names
  names <- vector("list", ns)
  for(t in 1:ns){
    names[[t]] <- dimnames(Y)[[1]]
  }

  ## load network data
  end <- c(which(diff(median.s) == 1), time)
  start <- c(1, which(diff(median.s) == 1)+1)
  N <- K[1]

  net <- array(NA, dim=c(N, N, ns))
  D <- t(sapply(1:time, function(t){V[t, order(V[t,], decreasing=TRUE)]}))
  Dmat <- data.frame(Year = Year, `1st` = D[,1], `2nd` = D[,2], check.names = FALSE)
  ## Use tidyr::pivot_longer instead of reshape::melt
  D.long <- tidyr::pivot_longer(Dmat, cols = c("1st", "2nd"),
                                names_to = "variable", values_to = "value")
  D.regime <- matrix(NA, ns, R)

  for(t in 1:ns){
    net[,,t] <- apply(Y[,,start[t]:end[t]], 1:2, sum)
    D.regime[t, ] <- apply(Dmat[start[t]:end[t], 2:3], 2, mean)
  }

  ## multiplot object
  U.list <- .df.list <- title.list <- time.period <- vector("list", ns)
  p.list <- vector("list", ns)
  median.s <- ceiling(apply(attr(mcmcout, "Smat"), 2, median))
  First <- Second <- Size <- Names <- Cluster <- NULL
  for(i in 1:ns){
    time.period[[i]] <- paste0(range(Year[median.s == i])[1], "-", range(Year[median.s == i])[2])
    U.list[[i]] <- U[[i]]
    cls <- kmeans(U.list[[i]], n.cluster[i], nstart = 20)
    .df.list[[i]] <- data.frame(First = U.list[[i]][,1], Second = U.list[[i]][,2],
                               Size = sqrt((U.list[[i]][,1])^2 + (U.list[[i]][,2])^2),
                               Names = names[[i]],
                               Cluster = factor(cls$cluster))
    title.list[[i]] <- paste0("Regime ", i, " (", time.period[[i]], ")")
    p.list[[i]] <- ggplot(.df.list[[i]], aes(x=First, y = Second, label=Names, color=Cluster)) +
      geom_point(aes(colour = Cluster), alpha = 0.5, size = point.cex, show.legend = FALSE) +
      scale_color_viridis_d(option = "D") +
      scale_size_continuous(guide = "none") +
      ggtitle(title.list[[i]]) +
      labs(x = paste("v[1] = ", round(D.regime[i, 1], 2)),
           y = paste("v[2] = ", round(D.regime[i, 2], 2))) +
      ggrepel::geom_text_repel(size = text.cex, segment.size = segment.size,
                               segment.color = alpha("navy", 1/6),
                               colour = "navy", aes(label = Names)) +
      theme_networkchange() +
      theme(axis.title = element_text(size = 8))
  }
  attr(p.list, "title.list") <- title.list
  attr(p.list, "D.long") <- D.long
  attr(p.list, "median.s") <- median.s
  return(p.list)
}