# Zigzag expanded navigation plots in R: The R package zenplots

#### M. Hofert and R. W. Oldford

#### 2021-09-08

This vignette accompanies the paper “Zigzag expanded navigation plots in R: The R package zenplots”. Note that sections are numbered accordingly (or omitted). Furthermore, it is recommended to read the paper to follow this vignette.

    # attaching required packages
    library(PairViz)
    library(MASS)
    library(zenplots)

## 2 Zenplots

As example data, we use the `olive` data set:

    data(olive, package = "zenplots")

Reproducing the plots of Figure 1:

    zenplot(olive)

    zenplot(olive, plot1d = "layout", plot2d = "layout")

Considering the `str()`ucture of `zenplot()` (here formatted for nicer output):

    str(zenplot)

    function (x, turns = NULL, first1d = TRUE, last1d = TRUE,
              n2dcols = c("letter", "square", "A4", "golden", "legal"),
              n2dplots = NULL,
              plot1d = c("label", "points", "jitter", "density", "boxplot",
                         "hist", "rug", "arrow", "rect", "lines", "layout"),
              plot2d = c("points", "density", "axes", "label", "arrow",
                         "rect", "layout"),
              zargs = c(x = TRUE, turns = TRUE, orientations = TRUE,
                        vars = TRUE, num = TRUE, lim = TRUE, labs = TRUE,
                        width1d = TRUE, width2d = TRUE,
                        ispace = match.arg(pkg) != "graphics"),
              lim = c("individual", "groupwise", "global"),
              labs = list(group = "G", var = "V", sep = ", ", group2d = FALSE),
              pkg = c("graphics", "grid", "loon"),
              method = c("tidy", "double.zigzag", "single.zigzag"),
              width1d = if (is.null(plot1d)) 0.5 else 1,
              width2d = 10,
              ospace = if (pkg == "loon") 0 else 0.02,
              ispace = if (pkg == "graphics") 0 else 0.037, draw = TRUE, ...)

### 2.1 Layout

To investigate the layout options of zenplots a bit more, we need a larger data set. To this end we simply double the olive data here (obviously only for illustration purposes):

    olive2 <- cbind(olive, olive) # just for this illustration

Reproducing the plots of Figure 2:

    zenplot(olive2, n2dcols = 6, plot1d = "layout", plot2d = "layout",
            method = "single.zigzag")

    zenplot(olive2, n2dcols = 6, plot1d = "layout", plot2d = "layout",
            method = "double.zigzag")

    zenplot(olive2, n2dcols = 6, plot1d = "layout", plot2d = "layout",
            method = "tidy")

Note that there is also `method = "rectangular"` (leaving the zigzagging zenplot paradigm but being useful for laying out 2d plots which are not necessarily connected through a variable; note that in this case, we omit the 1d plots as the default (labels) is rather confusing in this example):

    zenplot(olive2, n2dcols = 6, plot1d = "arrow", plot2d = "layout",
            method = "rectangular")

Reproducing the plots of Figure 3:

    zenplot(olive, plot1d = "layout", plot2d = "layout", method = "double.zigzag",
            last1d = FALSE, ispace = 0.1)

    zenplot(olive, plot1d = "layout", plot2d = "layout", n2dcol = 4, n2dplots = 8,
            width1d = 2, width2d = 4)

## 3 Zenpaths

A very basic path (standing for the sequence of pairs (1,2), (2,3), (3,4), (4,5)):

    (path <- 1:5)

    ## [1] 1 2 3 4 5

A zenpath through all pairs of variables (Eulerian):

    (path <- zenpath(5))

    ##  [1] 5 1 2 3 1 4 2 5 3 4 5

If `dataMat` is a five-column matrix, the zenplot of all pairs would then be constructed as follows:

    zenplot(x = dataMat[,path])

The `str()`ucture of `zenpath()` (again formatted for nicer output):

    str(zenpath)

    function (x, pairs = NULL,
              method = c("front.loaded", "back.loaded", "balanced",
                         "eulerian.cross", "greedy.weighted", "strictly.weighted"),
              decreasing = TRUE)

Here are some methods for five variables:

    zenpath(5, method = "front.loaded")

    ##  [1] 5 1 2 3 1 4 2 5 3 4 5

    zenpath(5, method = "back.loaded")

    ##  [1] 1 2 3 1 4 2 5 3 4 5 1

    zenpath(5, method = "balanced")

    ##  [1] 1 2 3 5 4 1 3 4 2 5 1

The following method considers two groups: One of size three, the other of size five. The sequence of pairs is constructed such that the first variable comes from the first group, the second from the second.

    zenpath(c(3,5), method = "eulerian.cross")

    ##  [1] 1 4 2 5 1 6 2 7 1 8 3 4 3 6 7 3 5 8 2

Reproducing Figure 4:

    oliveAcids <- olive[, !names(olive) %in% c("Area", "Region")] # acids only
    zpath <- zenpath(ncol(oliveAcids)) # all pairs
    zenplot(oliveAcids[, zpath], plot1d = "hist", plot2d = "density")

## 4 Build your own zenplots

### 4.3 Custom layout and plots – a spiral of ggplots example

Figure 5 can be reproduced as follows (note that we do not show the plot here due to a CRAN issue when running this vignette):

    path <- c(1,2,3,1,4,2,5,1,6,2,7,1,8,2,3,4,5,3,6,4,7,3,8,4,5,6,7,5,8,6,7,8)
    turns <- c("l",
               "d","d","r","r","d","d","r","r","u","u","r","r","u","u","r","r",
               "u","u","l","l","u","u","l","l","u","u","l","l","d","d","l","l",
               "u","u","l","l","d","d","l","l","d","d","l","l","d","d","r","r",
               "d","d","r","r","d","d","r","r","d","d","r","r","d","d")

    library(ggplot2) # for ggplot2-based 2d plots
    stopifnot(packageVersion("ggplot2") >= "2.2.1") # need 2.2.1 or higher
    ggplot2d <- function(zargs) {
      r <- extract_2d(zargs)
      num2d <- zargs$num/2
      df <- data.frame(x = unlist(r$x), y = unlist(r$y))
      p <- ggplot() +
        geom_point(data = df, aes(x = x, y = y), cex = 0.1) +
        theme(axis.line = element_blank(),
              axis.ticks = element_blank(),
              axis.text.x = element_blank(),
              axis.text.y = element_blank(),
              axis.title.x = element_blank(),
              axis.title.y = element_blank())
      if(num2d == 1) p <- p +
        theme(panel.background = element_rect(fill = 'royalblue3'))
      if(num2d == (length(zargs$turns)-1)/2) p <- p +
        theme(panel.background = element_rect(fill = 'maroon3'))
      ggplot_gtable(ggplot_build(p))
    }

    zenplot(as.matrix(oliveAcids)[,path], turns = turns, pkg = "grid",
            plot2d = function(zargs) ggplot2d(zargs))

### 4.4 Data groups

Split the olive data set into three groups (according to their variable `Area`):

    oliveAcids.by.area <- split(oliveAcids, f = olive$Area)
    # Replace the "." by " " in third group's name
    names(oliveAcids.by.area)[3] <- gsub("\\.", " ", names(oliveAcids.by.area)[3])
    names(oliveAcids.by.area)

    ## [1] "North-Apulia"    "Calabria"        "South-Apulia"    "Sicily"
    ## [5] "Inland-Sardinia" "Coast-Sardinia"  "East-Liguria"    "West-Liguria"
    ## [9] "Umbria"

Reproducing the plots of Figure 6 (note that `lim = "groupwise"` does not make much sense here as a plot):

    zenplot(oliveAcids.by.area, labs = list(group = NULL))

    zenplot(oliveAcids.by.area, lim = "groupwise", labs = list(sep = " - "),
            plot1d = function(zargs) label_1d_graphics(zargs, cex = 0.8),
            plot2d = function(zargs)
                points_2d_graphics(zargs, group... = list(sep = "\n - \n")))

### 4.5 Custom zenpaths

Find the “convexity” scagnostic for each pair of olive acids.

    library(scagnostics)
    Y <- scagnostics(oliveAcids) # compute scagnostics (scatter-plot diagonstics)
    X <- Y["Convex",] # pick out component 'convex'
    d <- ncol(oliveAcids)
    M <- matrix(NA, nrow = d, ncol = d) # matrix with all 'convex' scagnostics
    M[upper.tri(M)] <- X # (i,j)th entry = scagnostic of column pair (i,j) of oliveAcids
    M[lower.tri(M)] <- t(M)[lower.tri(M)] # symmetrize
    round(M, 5)

    ##         [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]
    ## [1,]      NA 0.48952 0.46343 0.45887 0.43914 0.34583 0.31259 0.28413
    ## [2,] 0.48952      NA 0.42276 0.50499 0.44591 0.35855 0.35846 0.31729
    ## [3,] 0.46343 0.42276      NA 0.39700 0.36394 0.31316 0.29534 0.33709
    ## [4,] 0.45887 0.50499 0.39700      NA 0.46454 0.36616 0.29451 0.34888
    ## [5,] 0.43914 0.44591 0.36394 0.46454      NA 0.31977 0.31443 0.36750
    ## [6,] 0.34583 0.35855 0.31316 0.36616 0.31977      NA 0.53726 0.34001
    ## [7,] 0.31259 0.35846 0.29534 0.29451 0.31443 0.53726      NA 0.22231
    ## [8,] 0.28413 0.31729 0.33709 0.34888 0.36750 0.34001 0.22231      NA

Show the six pairs with largest “convexity” scagnostic:

    zpath <- zenpath(M, method = "strictly.weighted") # list of ordered pairs
    head(M[do.call(rbind, zpath)]) # show the largest six 'convexity' measures

    ## [1] 0.5372599 0.5049945 0.4895179 0.4645377 0.4634277 0.4588675

Extract the corresponding pairs:

    (ezpath <- extract_pairs(zpath, n = c(6, 0))) # extract the first six pairs

    ## [[1]]
    ## [1] 7 6
    ##
    ## [[2]]
    ## [1] 4 2
    ##
    ## [[3]]
    ## [1] 2 1
    ##
    ## [[4]]
    ## [1] 5 4
    ##
    ## [[5]]
    ## [1] 3 1
    ##
    ## [[6]]
    ## [1] 4 1

Reproducing Figure 7 (visualizing the pairs):

    library(graph)
    library(Rgraphviz)
    plot(graph_pairs(ezpath)) # depict the six most convex pairs (edge = pair)

Connect them:

    (cezpath <- connect_pairs(ezpath)) # keep the same order but connect the pairs

    ## [[1]]
    ## [1] 7 6
    ##
    ## [[2]]
    ## [1] 4 2 1
    ##
    ## [[3]]
    ## [1] 5 4
    ##
    ## [[4]]
    ## [1] 3 1 4

Build the corresponding list of matrices:

    oliveAcids.grouped <- groupData(oliveAcids, indices = cezpath) # group data for (zen)plotting

Reproducing Figure 8 (zenplot of the six pairs of acids with largest “convexity” scagnostic):

    zenplot(oliveAcids.grouped)

## 5 Advanced features

### 5.1 The structure of a zenplot

Here is the structure of a return object of `zenplot()`:

    res <- zenplot(olive, plot1d = "layout", plot2d = "layout", draw = FALSE)
    str(res)

    ## List of 2
    ##  $ path  :List of 3
    ##   ..$ turns    : chr [1:19] "d" "r" "r" "d" ...
    ##   ..$ positions: num [1:19, 1:2] 1 2 2 2 3 4 4 4 5 6 ...
    ##   .. ..- attr(*, "dimnames")=List of 2
    ##   .. .. ..$ : NULL
    ##   .. .. ..$ : chr [1:2] "x" "y"
    ##   ..$ occupancy: chr [1:8, 1:6] "" "" "" "" ...
    ##  $ layout:List of 6
    ##   ..$ orientations : chr [1:19] "h" "s" "v" "s" ...
    ##   ..$ dimensions   : num [1:19] 1 2 1 2 1 2 1 2 1 2 ...
    ##   ..$ vars         : num [1:19, 1:2] 1 1 2 3 3 3 4 5 5 5 ...
    ##   .. ..- attr(*, "dimnames")=List of 2
    ##   .. .. ..$ : NULL
    ##   .. .. ..$ : chr [1:2] "x" "y"
    ##   ..$ layoutWidth  : num 33
    ##   ..$ layoutHeight : num 44
    ##   ..$ boundingBoxes: num [1:19, 1:4] 0 0 10 11 11 11 21 22 22 22 ...
    ##   .. ..- attr(*, "dimnames")=List of 2
    ##   .. .. ..$ : NULL
    ##   .. .. ..$ : chr [1:4] "left" "right" "bottom" "top"
    ##  - attr(*, "class")= chr [1:3] "zenGraphics" "zenplot" "list"

Let’s have a look at the components. The occupancy matrix encodes the occupied cells in the rectangular layout:

    res[["path"]][["occupancy"]]

    ##      [,1] [,2] [,3] [,4] [,5] [,6]
    ## [1,] ""   "d"  ""   ""   ""   ""
    ## [2,] ""   "r"  "r"  "d"  ""   ""
    ## [3,] ""   ""   ""   "d"  ""   ""
    ## [4,] ""   ""   ""   "r"  "r"  "d"
    ## [5,] ""   ""   ""   ""   ""   "d"
    ## [6,] "l"  "l"  ""   "d"  "l"  "l"
    ## [7,] ""   "u"  ""   "d"  ""   ""
    ## [8,] ""   "u"  "l"  "l"  ""   ""

The two-column matrix `positions` contains in the ith row the row and column index (in the occupancy matrix) of the ith plot:

    head(res[["path"]][["positions"]])

    ##      x y
    ## [1,] 1 2
    ## [2,] 2 2
    ## [3,] 2 3
    ## [4,] 2 4
    ## [5,] 3 4
    ## [6,] 4 4

### 5.2 Tools for writing 1d and 2d plot functions

Example structure of 2d plot based on `graphics`:

    points_2d_graphics

    ## function (zargs, cex = 0.4, box = FALSE, add = FALSE, group... = NULL,
    ##     plot... = NULL, ...)
    ## {
    ##     r <- extract_2d(zargs)
    ##     xlim <- r$xlim
    ##     ylim <- r$ylim
    ##     x <- as.matrix(r$x)
    ##     y <- as.matrix(r$y)
    ##     same.group <- r$same.group
    ##     if (same.group) {
    ##         if (!add)
    ##             plot_region(xlim = xlim, ylim = ylim, plot... = plot...)
    ##         points(x = x, y = y, cex = cex, ...)
    ##         if (box)
    ##             box(...)
    ##     }
    ##     else {
    ##         args <- c(list(zargs = zargs, add = add), group...)
    ##         do.call(group_2d_graphics, args)
    ##     }
    ## }
    ## <bytecode: 0x7fc2480d4db0>
    ## <environment: namespace:zenplots>

For setting up the plot region of plots based on `graphics`:

    plot_region

    ## function (xlim, ylim, plot... = NULL)
    ## {
    ##     if (is.null(plot...)) {
    ##         plot(NA, xlim = xlim, ylim = ylim, type = "n", ann = FALSE,
    ##             axes = FALSE, log = "")
    ##     }
    ##     else {
    ##         fun <- function(...) plot(NA, xlim = xlim, ylim = ylim,
    ##             ...)
    ##         do.call(fun, plot...)
    ##     }
    ## }
    ## <bytecode: 0x7fc27fe8d5d0>
    ## <environment: namespace:zenplots>

Determining the indices of the two variables to be plotted in the current 1d or 2d plot (the same for 1d plots):

    plot_indices

    ## function (zargs)
    ## zargs$vars[zargs$num, ]
    ## <bytecode: 0x7fc23fe8c7b8>
    ## <environment: namespace:zenplots>

Basic check that the return value of `zenplot()` is actually the return value of the underlying `unfold()` (note that, the output of `unfold` and `res` is not identical since `res` has specific class attributes):

    n2dcols <- ncol(olive) - 1 # number of faces of the hypercube
    uf <- unfold(nfaces = n2dcols)

    identical(res, uf) #return FALSE

    ## [1] FALSE

    for(name in names(uf)) {
       stopifnot(identical(res[[name]], uf[[name]]))
    }