R/mosaic.plotly.R

In lessR: Less Code with More Comprehensive Results

# Mosaic plot and association plot, Plotly only
# Called from: Chart.R (form="mosaic", form="assoc")


# Contingency table of x and the by variables, the counts summed from y
#   when given, as frequencies of a summary table. Cases with a missing
#   value of any variable are dropped.
.mosaic_table <- function(x, by, y, x_name, by_names) {

  vars <- list(x)
  if (!is.null(by))
    vars <- c(vars, if (is.data.frame(by)) as.list(by) else list(by))
  names(vars) <- c(x_name, by_names)
  vars <- lapply(vars, function(v) if (is.factor(v)) v else factor(v))

  ok <- complete.cases(as.data.frame(vars))
  w <- if (is.null(y)) rep(1, length(vars[[1L]])) else as.numeric(y)
  ok <- ok & !is.na(w)

  tab <- tapply(w[ok], lapply(vars, function(v) v[ok]), sum)
  tab[is.na(tab)] <- 0
  tab
}


# Expected counts under mutual independence of all the variables, which
#   for two variables is the usual independence model, with the Pearson
#   residuals and the chi-square test of that model
.mosaic_resid <- function(tab) {

  n <- sum(tab)
  dm <- dim(tab)
  marg <- lapply(seq_along(dm), function(k) apply(tab, k, sum) / n)
  e <- n * Reduce(outer, marg)
  dim(e) <- dm
  dimnames(e) <- dimnames(tab)

  r <- (tab - e) / sqrt(e)
  r[e == 0] <- 0

  df <- prod(dm) - 1 - sum(dm - 1)
  chisq <- sum(r^2)
  p <- if (df > 0) pchisq(chisq, df, lower.tail=FALSE) else NA

  list(e=e, r=r, chisq=chisq, df=df, p=p)
}


# Fill for each residual, the shading_hcl() scheme of the vcd package:
#   positive residuals in the hue of the theme, negative in its
#   complement, more saturated and darker past |r| of 2 and of 4
# for the default theme the pair is the first two colors of "hues",
#   blue and brown, 180 degrees apart; hue NA, the gray and white
#   themes, has no chroma
.mosaic_shade <- function(r, hue) {
  i <- ifelse(abs(r) > 4, 1, ifelse(abs(r) > 2, 0.5, 0))
  chroma <- if (is.na(hue)) 0 else 20 + 80*i
  h <- if (is.na(hue)) 0 else ifelse(r >= 0, hue, (hue + 180) %% 360)
  hcl(h, chroma, 90 - 40*i)
}


# HCL hue of the first fill color of the theme
.mosaic_hue <- function() {
  clr <- substr(.to_hex(getOption("bar_fill_discrete")[1L]), 1L, 7L)
  luv <- convertColor(t(col2rgb(clr)) / 255, from="sRGB", to="Luv")
  (atan2(luv[, 3L], luv[, 2L]) * 180 / pi) %% 360
}


# the residual key as legend entries, one per shading level, each drawn
#   outside the plotted range so it appears only in the legend
.mosaic_key <- function(p, hue) {

  r.lvl <- c(5, 3, 0, -3, -5)
  lbl <- c("> 4", "2 to 4", "-2 to 2", "-4 to -2", "< -4")
  if (is.na(hue)) {  # the sign shows only in the border, so one entry a level
    r.lvl <- r.lvl[1:3]
    lbl <- c("|r| > 4", "|r| 2 to 4", "|r| < 2")
  }
  clr <- .mosaic_shade(r.lvl, hue)
  for (i in seq_along(r.lvl))
    p <- plotly::add_trace(p, x=-10, y=-10, type="scatter",
           mode="markers", name=lbl[i], showlegend=TRUE,
           hoverinfo="skip", inherit=FALSE,
           marker=list(symbol="square", size=14, color=clr[i],
                       line=list(color="#4D4D4D", width=1)))
  p
}


.mosaic.plotly <- function(tab, form="mosaic", fill=NULL,
                           residuals=FALSE, theme="colors",
                           main=NULL, digits_d=2,
                           style_opts=.plotly_style()) {

  rs <- .mosaic_resid(tab)
  nd <- length(dim(tab))
  var.nm <- names(dimnames(tab))
  lvl <- dimnames(tab)
  gray <- theme %in% c("gray", "white")
  hue <- if (gray) NA else .mosaic_hue()
  # the variable names a step above the level labels, not two
  lab.size <- round(13 * getOption("lab_size", 1))
  ax.size <- round(13 * getOption("axis_size", 0.9))
  # the side variable names sit beyond the widest level label, in data
  #   units of a panel of about 300 pixels, the smaller end of a typical
  #   display: a character about 0.8 of the axis font size wide, enough
  #   for capitals, a name line about 1.4 of the label size
  chr.w <- 0.8 * ax.size / 300
  lab.h <- 1.4 * lab.size / 300
  side.x <- function(lv, left)
    if (left) -0.015 - max(nchar(lv)) * chr.w - 0.8 * lab.h
    else 1.015 + max(nchar(lv)) * chr.w + 0.8 * lab.h
  top.y <- 1.015 + 1.5 * ax.size / 300
  edge <- "#4D4D4D"  # gray30

  # fill mode: residuals=TRUE shades each tile by its Pearson residual;
  #   otherwise, with two or more variables, one fill for each level of
  #   the last variable, recycled as for a bar chart
  shade <- residuals
  if (shade && nd < 2) {
    cat("\n"); stop(call.=FALSE, "\n------\n",
      "residuals=TRUE  shades each tile by its departure from\n",
      "  independence, which needs at least two variables: specify\n",
      "  by  as well as  x.\n\n")
  }
  highlight <- nd > 1 && !shade && !is.null(fill)
  if (highlight) fill <- .to_hex(rep_len(fill, length(lvl[[nd]])))
  else if (!shade && !is.null(fill)) fill <- .to_hex(fill)

  # each cell of the table as a row: its levels, observed, expected, r
  cells <- expand.grid(lvl, stringsAsFactors=FALSE)
  cells$o <- as.vector(tab)
  cells$e <- as.vector(rs$e)
  cells$r <- as.vector(rs$r)
  cell.fill <- function(k) {  # k indexes cells
    if (highlight) fill[match(cells[k, nd], lvl[[nd]])]
    else if (shade) .mosaic_shade(cells$r[k], hue)
    else (fill %||% getOption("bar_fill_discrete", "gray80"))[1L]
  }
  hover <- function(k) {
    txt <- paste0(var.nm, ": ", unlist(cells[k, seq_len(nd)]),
                  collapse="<br>")
    paste0(txt, "<br>Count: ", cells$o[k],
           if (nd > 1) paste0("<br>Expected: ",
                              .fmt(cells$e[k], digits_d),
                              "<br>Residual: ",
                              .fmt(cells$r[k], digits_d)))
  }
  cell.of <- function(idx)  # row of cells for a vector of level indices
    1L + sum((idx - 1L) * cumprod(c(1L, dim(tab)))[seq_len(nd)])

  shapes <- list()
  ann <- list()
  hx <- hy <- ht <- NULL
  rect <- function(x0, x1, y0, y1, k) {
    dash <- if (shade && cells$r[k] < 0) "dash" else "solid"
    shapes[[length(shapes) + 1L]] <<- list(type="rect",
      x0=x0, x1=x1, y0=y0, y1=y1, fillcolor=cell.fill(k),
      line=list(color=edge, width=1, dash=dash), layer="below")
    hx <<- c(hx, (x0 + x1) / 2)
    hy <<- c(hy, (y0 + y1) / 2)
    ht <<- c(ht, hover(k))
  }
  lab <- function(x, y, text, xa="center", ya="middle", angle=0,
                  size=ax.size)
    ann[[length(ann) + 1L]] <<- list(x=x, y=y, text=text, xref="x",
      yref="y", xanchor=xa, yanchor=ya, textangle=angle,
      showarrow=FALSE, font=list(size=size))


  if (form == "mosaic") {

    gap <- c(0.03, 0.02, 0.012, rep(0.008, max(0, nd - 3)))
    margins <- lapply(seq_len(nd), function(k) apply(tab, seq_len(k), sum))
    pos <- vector("list", nd)  # piece positions per level, for labels

    split <- function(x0, x1, y0, y1, k, idx) {
      if (k > nd) {
        rect(x0, x1, y0, y1, cell.of(idx))
        return(invisible())
      }
      m <- margins[[k]]
      sub <- if (k == 1L) as.vector(m)
             else as.vector(do.call(`[`, c(list(m), as.list(idx), TRUE)))
      L <- length(sub)
      p <- if (sum(sub) > 0) sub / sum(sub) else rep(0, L)
      rows <- k %% 2L == 1L  # odd levels stack rows, even place columns
      len <- if (rows) y1 - y0 else x1 - x0
      g <- min(gap[k], len / (2 * L))
      sz <- p * (len - g * (L - 1))
      at <- if (rows) y1 else x0
      for (i in seq_len(L)) {
        a0 <- at
        a1 <- if (rows) at - sz[i] else at + sz[i]
        pos[[k]][[length(pos[[k]]) + 1L]] <<-
          list(i=i, idx=idx, a=(a0 + a1) / 2, b0=if (rows) x0 else y0,
               b1=if (rows) x1 else y1)
        if (rows) split(x0, x1, a1, a0, k + 1L, c(idx, i))
        else split(a0, a1, y0, y1, k + 1L, c(idx, i))
        at <- if (rows) a1 - g else a1 + g
      }
    }
    split(0, 1, 0, 1, 1L, integer(0))

    # level labels: the first variable down the left, the second across
    #   the top over the first row, the third down the right beside the
    #   last column, each with its variable name
    for (q in pos[[1L]]) lab(-0.015, q$a, lvl[[1L]][q$i], xa="right")
    x.lo <- side.x(lvl[[1L]], TRUE)
    lab(x.lo, 0.5, var.nm[1L], angle=-90, size=lab.size)
    x.hi <- 1.02
    if (nd >= 2L) {
      for (q in pos[[2L]]) if (q$idx[1L] == 1L)
        lab(q$a, 1.015, lvl[[2L]][q$i], ya="bottom")
      lab(0.5, top.y, var.nm[2L], ya="bottom", size=lab.size)
    }
    if (nd >= 3L) {
      n2 <- length(lvl[[2L]])
      for (q in pos[[3L]]) if (q$idx[2L] == n2)
        lab(1.015, q$a, lvl[[3L]][q$i], xa="left")
      x.hi <- side.x(lvl[[3L]], FALSE)
      lab(x.hi, 0.5, var.nm[3L], angle=90, size=lab.size)
      x.hi <- x.hi + lab.h
    }
    x.rng <- c(x.lo - lab.h, x.hi)
    y.rng <- c(-0.02, if (nd >= 2L) top.y + lab.h else 1.02)
  }


  else {  # association plot: rows for the combinations of the leading
          #   variables, a column for each level of the last

    if (nd < 2L) {
      cat("\n"); stop(call.=FALSE, "\n------\n",
        "An association plot displays the departure of each cell from\n",
        "  independence, so it needs at least two variables: specify\n",
        "  by  as well as  x.\n\n")
    }
    row.vars <- seq_len(nd - 1L)
    # rows in the order of the first variable, then the next within it
    ord <- do.call(order, lapply(row.vars,
                                 function(v) match(cells[[v]], lvl[[v]])))
    row.key <- do.call(paste, c(cells[row.vars], sep="\r"))
    rows <- unique(row.key[ord])
    n.col <- length(lvl[[nd]])
    col.w <- 1 / n.col
    w.max <- max(sqrt(cells$e))

    hp <- vapply(rows, function(k) max(0, cells$r[row.key == k]), 0)
    hn <- vapply(rows, function(k) max(0, -cells$r[row.key == k]), 0)
    band <- pmax(hp + hn, 1e-6)
    g <- 0.15 * sum(band) / length(rows)
    s <- 1 / (sum(band) + g * (length(rows) - 1))
    top <- 1
    for (j in seq_along(rows)) {
      base <- top - hp[j] * s
      shapes[[length(shapes) + 1L]] <- list(type="line", x0=0, x1=1,
        y0=base, y1=base, line=list(color="gray50", width=1, dash="dot"),
        layer="below")
      for (k in which(row.key == rows[j])) {
        cc <- match(cells[k, nd], lvl[[nd]])
        w <- 0.9 * col.w * sqrt(cells$e[k]) / w.max
        xm <- (cc - 0.5) * col.w
        y.tip <- base + cells$r[k] * s
        rect(xm - w/2, xm + w/2, min(base, y.tip), max(base, y.tip), k)
      }
      k1 <- which(row.key == rows[j])[1L]
      # the last row variable labels each row on the right; the first
      #   labels its group of rows on the left
      if (nd >= 3L)
        lab(1.015, base, cells[k1, nd - 1L], xa="left")
      top <- top - band[j] * s - g * s
    }
    first <- vapply(rows, function(k) strsplit(k, "\r")[[1L]][1L], "")
    for (lv in unique(first)) {
      js <- which(first == lv)
      ys <- vapply(js, function(j) 1 - (sum(band[seq_len(j)]) -
                   band[j]/2 + g*(j - 1)) * s, 0)
      lab(-0.015, mean(ys), lv, xa="right")
    }
    x.lo <- side.x(lvl[[1L]], TRUE)
    lab(x.lo, 0.5, var.nm[1L], angle=-90, size=lab.size)
    x.hi <- 1.02
    if (nd >= 3L) {
      x.hi <- side.x(lvl[[nd - 1L]], FALSE)
      lab(x.hi, 0.5, var.nm[nd - 1L], angle=90, size=lab.size)
      x.hi <- x.hi + lab.h
    }
    for (cc in seq_len(n.col))
      lab((cc - 0.5) * col.w, 1.015, lvl[[nd]][cc], ya="bottom")
    lab(0.5, top.y, var.nm[nd], ya="bottom", size=lab.size)
    x.rng <- c(x.lo - lab.h, x.hi)
    y.rng <- c(-0.02, top.y + lab.h)
  }


  p <- plotly::plot_ly(x=hx, y=hy, text=ht, type="scatter", mode="markers",
         hoverinfo="text", showlegend=FALSE,
         marker=list(size=12, opacity=0))

  if (shade) p <- .mosaic_key(p, hue)
  if (highlight)
    for (i in seq_along(fill))
      p <- plotly::add_trace(p, x=-10, y=-10, type="scatter",
             mode="markers", name=lvl[[nd]][i], showlegend=TRUE,
             hoverinfo="skip", inherit=FALSE,
             marker=list(symbol="square", size=14, color=fill[i],
                         line=list(color=edge, width=1)))

  if (shade) {  # the test of independence beneath the key
    p.txt <- if (is.na(rs$p)) "" else
             if (rs$p < 0.001) "p-value < 0.001"
             else paste("p-value =", .fmt(rs$p, 3))
    ann[[length(ann) + 1L]] <- list(x=1.02, y=0, xref="paper",
      yref="paper", xanchor="left", yanchor="bottom", showarrow=FALSE,
      align="left", font=list(size=ax.size),
      text=paste0(if (gray) "Dashed border:<br>negative<br><br>",
                  p.txt))
  }

  ax <- list(visible=FALSE, fixedrange=TRUE)
  title.size <- round(16 * getOption("main_size", 1))
  p <- plotly::layout(p,
    title=if (!is.null(main) && nzchar(main))
            list(text=main, font=list(size=title.size)) else NULL,
    xaxis=c(ax, list(range=x.rng)),
    yaxis=c(ax, list(range=y.rng, scaleanchor="x")),
    shapes=shapes, annotations=ann,
    showlegend=shade || highlight,
    legend=list(x=1.02, xanchor="left", y=1, yanchor="top",
      title=list(text=if (shade) "Pearson<br>residuals" else var.nm[nd],
                 font=list(size=lab.size)),
      itemclick=FALSE, itemdoubleclick=FALSE,
      bordercolor=.to_hex(style_opts$panel_border), borderwidth=1,
      font=list(size=ax.size)),
    plot_bgcolor=.to_hex(style_opts$window_fill),
    paper_bgcolor=.to_hex(style_opts$window_fill),
    margin=list(t=if (is.null(main)) 20 else round(title.size * 2.2),
                l=10, r=10, b=10))

  p <- .plotly_font(p)
  list(plot=p, resid=rs, tab=tab)
}