Vertical Histogram

In the process of munging data for my current project I came across the need to compare (visually) the difference between two modes within the same dataset. I was using a simple scatterplot and setting the alpha in the hopes that the over-plotting would indicate which was the major mode. Unfortunately, the size of the data overwhelmed this approach.

I only wanted to use a single image and it was important that I keep the scatterplot to show other features of the data. I started looking for a way to combine a histogram (rotated 90 degrees) with the scatterplot to help describe the density within the plot. A quick search for how to do this in R turned up empty so I decided to implement my own version of such a plot.

Certainly, there are other ways to describe the features that I am trying to present here but in this particular case the following code worked out nicely. Hopefully it proves useful to others as well.

plot.vertical.hist <- function(data,breaks=500) {

    agg <- aggregate(data$Y, by=list(xs=data$X), FUN=mean)
    hs <- hist(agg$x / 10000, breaks=breaks, plot=FALSE)

    old.par <- par(no.readonly=TRUE)
    mar.default <- par('mar')
    mar.left <- mar.default
    mar.right <- mar.default
    mar.left[4] <- 0
    mar.right[2] <- 0

    # Main plot 
    par (fig=c(0,0.8,0,1.0), mar=mar.left)
    plot (agg$xs, agg$x / 10000,
          xlab="X", ylab="Y",
          main="Vertical Histogram Side Plot",
          pch=19, col=rgb(0.5,0.5,0.5,alpha=0.5))
    grid ()

    # Vertical histogram of the same data
    par (fig=c(0.8,1.0,0.0,1.0), mar=mar.right, new=TRUE)
    plot (NA, type='n', axes=FALSE, yaxt='n',
          xlab='Frequency', ylab=NA, main=NA,
          xlim=c(0,max(hs$counts)),
          ylim=c(1,length(hs$counts)))
    axis (1)
    arrows(rep(0,length(hs$counts)),1:length(hs$counts),
           hs$counts,1:length(hs$counts),
           length=0,angle=0)

    par(old.par)
    invisible ()
}

Results look similar to the following:

Initially, I experimented with rug or barplot(..., horiz=TRUE). Unfortunately, rug isn't available on the left or right side and would suffer from the same problem that the alpha settings did and I was unable to get the alignment worked out when using barplot.

Column-Major Confusion

As a programmer coming from a language like C I am used to understanding multidimensional arrays in the following way:

int matrix[3][3] = {
    { 1, 2, 3 },
    { 4, 5, 6 },
    { 7, 8, 9 }
};

Understandably, this is a little bit of a setback when trying to grok how a language such as R handles matrices. The example above is 'row-major' but R uses 'column-major' by default. Note that I'm not describing memory layout here - which coincidentally is the same as my description - I referring to how the matrices are presented to the programmer. To complete the example, here is how I would create the same matrix as above in R:

> matrix(1:9,ncol=3)
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

There is no problem with this per se. In fact, I'd imagine that R programmers that come to something like C may feel similar unease in the paradigm shift. I'm finding that this fact is glazed over in some texts offering an introduction to R. This happens in non-obvious ways. In particular, I found this example today:

> matrix(c(1,0,0,0,1,0,0,0,1),ncol=3)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1

The identity matrix is a particularly bad choice in this regard as it gives no indication of the true layout being used. It is probably good for any tutorial using matrices to cover an obvious simple case first to set the stage before moving directly to something like the identity matrix (or any other symmetric matrix for that matter).

Go vector or go home

My programming experience progressed mostly along the lines of: C, C++, shell, Java, Java, Ruby, Python, Java, Java. Only recently have I started exploring the likes of Haskell, Erlang and R. Well that evolution bit me a little while back when I tried using an iterative approach in R where a vectorized solution would have been better.

I was dealing with a vector of timestamps that were formatted as 'seconds since the epoch' and what I wanted was to limit that vector to weekend timestamps only.

My naive approach was to construct a simple loop over the values and apply a function to each element. I was only dealing with about 20,000 elements but the time to do this was painfully slow - roughly 20 seconds - so I started investigating an apply-like approach. R provides several ways to do this depending on the input/output requirements: lapply, sapply, and vapply. All three resulted in behavior similar to the simple loop.

The function to test for weekend-ness is as follows:

is.weekend <- function(x) {
   tm <- as.POSIXlt(x,origin="1970/01/01")
   format(tm,"%a") %in% c("Sat","Sun")
}

I don't know the specific details of date/time conversion in R but I was pretty sure that this was not the bottleneck. After a little searching I came upon a different approach. Instead of looping over each element I should have been passing the entire vector around to the functions. I believe that the apply functions take the vector as an argument but do the manual loop internally. However, R supports a more native approach to handling vectors: vectorized operations.

Looping:

use.sapply <- function(data) {
   data[sapply(data$TIME,FUN=is.weekend),]
}

system.time(v <- use.sapply(csv.data))
  user  system elapsed
19.456   6.492  25.951

Vectorized:

use.vector <- function(data) {
   data[is.weekend(data$TIME),]
}

system.time(v <- use.vector(csv.data))
 user  system elapsed
0.032   0.020   0.052

Duly noted.

Hotness

We have an internal image that floated around work several years ago that details network utilization of TCP over a wide variety of configurations. It is a heatmap created in matlab that is just sweet, sweet eye candy. We actually hung it on the outside of a cube for a short while and people couldn't help but stop and look at it.

It is entirely dysfunctional, mind you. The designer tried to combine eight parameters - with all variations - into a individual 2D plot (3D if you consider color a dimension). It was definitely an internal tool - there were only two or three of us who could decipher the layout enough to say anything about the data. That was fine by us; we basically made up the entire population of people who cared.

Fast forward a few years. I'm currently working on a technical report that could use the data we used to create that plot and, as luck would have it, I'm also the only one from the original group still at the company. In order to be able to include this data in the report there needs to occur a certain amount of reformatting - first my brain, then that plot. I wasn't the original designer and, although I have access to the code, I don't know matlab so I'm pretty much stuck. I decided to rework the data in R.

The thing about that original plot was that it had a certain je ne sais quoi: it made you look. I wanted to keep that so I immediately investigated heatmap functionality available in R.

Really? Ouch. Not much available there. I came up with two resources that were helpful: A Wikipedia entry about a Mandelbrot set animation in R; A stackoverflow answer that mentioned rasterImage in a comment. The first site lead me to the color set used in our original plot and the second gave me the pointer I needed to get the job done. I'll leave what follows as a reminder for myself and a helpful nudge for those who face a similar problem in the future.

---

hmap.example <- function () {

    code <- c("colfun <- colorRampPalette(c(...))",
        "my.colors <- colfun(10000)","xs <- 1:100",
        "X  <- outer(xs,rev(xs))",
        "C1 <- matrix(my.colors[X],100,100)", "X  <- outer(xs,xs)",
        "C2 <- matrix(my.colors[X],100,100)", "X  <- outer(rev(xs),xs)",
        "C3 <- matrix(my.colors[X],100,100)",
        "plot(c(-100,100),c(-100,100),type='n')",
        "rasterImage(C1,1,1,100,100)",
        "rasterImage(C2,-100,1,1,100)", "rasterImage(C3,-100,-100,1,1)",
        "abline(v=0,col='black',lwd=5)", "abline(h=0,col='black',lwd=5)")

    colfun <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
                    "#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))

    my.colors <- colfun(10000)
    xs <- 1:100
    X  <- outer(xs,rev(xs))
    C1 <- matrix(my.colors[X],100,100)
    X  <- outer(xs,xs)
    C2 <- matrix(my.colors[X],100,100)
    X  <- outer(rev(xs),xs)
    C3 <- matrix(my.colors[X],100,100)
    plot(c(-100,100),c(-100,100),type='n',axes=FALSE,xlab='',ylab='')
    rasterImage(C1,1,1,100,100)
    rasterImage(C2,-100,1,1,100)
    rasterImage(C3,-100,-100,1,1)
    abline(v=0,col='black',lwd=5)
    abline(h=0,col='black',lwd=5)
    text(1,1:length(code)*-6,labels=code,cex=0.8,pos=4,family="mono")
}

---

And the result: