8

For very heavy-tailed data of both positive and negative sign, I sometimes like to see all the data on a plot without hiding structure in the unit interval.

When plotting with Matplotlib in Python, I can achieve this by selecting a symlog scale, which uses a logarithmic transform outside some interval, and linear plotting inside it.

Previously in R I have constructed similar behavior by transforming the data with an arcsinh on a one-off basis. However, tick labels and the like are very tricky to do right (see below). enter image description here

Now, I am faced with a bunch of data where the subsetting in lattice or ggplot would be highly convenient. I don't want to use Matplotlib because of the subsetting, but I sure am missing symlog!

Edit:

I see that ggplot uses a package called scales, which solves a lot of this problem (if it works). Automatically choosing tick mark and label placing still looks pretty hard to do nicely though. Some combination of log_breaks and cbreaks perhaps?

Edit 2:

The following code is not too bad

sinh.scaled <- function(x,scale=1){ sinh(x)*scale }
asinh.scaled <- function(x,scale=1) { asinh(x/scale) }



asinh_breaks <- function (n = 5, scale = 1, base=10) 
{
    function(x) {
        log_breaks.callable <- log_breaks(n=n,base=base)
        rng <- rng <- range(x, na.rm = TRUE)
        minx <- floor(rng[1])
        maxx <- ceiling(rng[2])
        if (maxx == minx) 
            return(sinh.scaled(minx, scale=scale))
        big.vals <- 0
        if (minx < (-scale)) {
            big.vals = big.vals + 1
        }
        if (maxx>scale) {
            big.vals = big.vals + 1
        }
        brk <- c()
        if (minx < (-scale)) {
            rbrk <- log_breaks.callable(  c(-min(maxx,-scale), -minx ) )
            rbrk <- -rev(rbrk)
            brk <- c(brk,rbrk)
        }
        if ( !(minx>scale | maxx<(-scale))  ) {
            rng <- c(max(minx,-scale), min(maxx,scale))
            minc <- floor(rng[1])
            maxc <- ceiling(rng[2])
            by <- floor((maxc - minc)/(n-big.vals)) + 1
            cb <- seq(minc, maxc, by = by)
            brk <- c(brk,cb)
        } 
        if (maxx>scale) {
            brk <- c(brk,log_breaks.callable( c(max(minx,scale), maxx )))
        }

        brk

    }
}

asinh_trans <- function(scale = 1) {
    trans <- function(x) asinh.scaled(x, scale)
    inv <- function(x) sinh.scaled(x, scale)
    trans_new(paste0("asinh-", format(scale)), trans, inv, 
              asinh_breaks(scale = scale), 
              domain = c(-Inf, Inf))
}
Community
  • 1
  • 1
Brian B
  • 1,410
  • 1
  • 16
  • 30
  • Usually it ends up being a TON easier to apply your log expansion/contraction to your data before plotting, and plot the resultant data in linear space. Add a few carefully-selected background lines (vertical or horizontal) to highlight the expansion zones, create a sequence of numbers to write to the axes' tickmark locations, and you should be all set. – Carl Witthoft Jan 30 '13 at 21:11
  • 1
    That said (con't what @CarlWitthoft wrote) it shouldn't be that hard to construct a scale transform function that mimics symlog if you look at how the other trans functions work in the **scales** package. – joran Jan 30 '13 at 21:18
  • Well, my thinking was that you'd write your own breaks function that would simply use the built in methods for finding breaks, but applied to each of the three regions separately. – joran Jan 30 '13 at 21:36
  • 2
    See the penultimate post in this thread by Brian Diggs for an example of how to write a function that defines a scale, in this case biexponential, which I think may have some resemblance to your problem: https://groups.google.com/forum/?fromgroups=#!topic/ggplot2/7ddCyXGlKiM – Dennis Jan 31 '13 at 10:11

1 Answers1

15

A solution based on the package scales and inspired by Brian Diggs' post mentioned by @Dennis:

symlog_trans <- function(base = 10, thr = 1, scale = 1){
  trans <- function(x)
    ifelse(abs(x) < thr, x, sign(x) * 
             (thr + scale * suppressWarnings(log(sign(x) * x / thr, base))))

  inv <- function(x)
    ifelse(abs(x) < thr, x, sign(x) * 
             base^((sign(x) * x - thr) / scale) * thr)

  breaks <- function(x){
    sgn <- sign(x[which.max(abs(x))])
    if(all(abs(x) < thr))
      pretty_breaks()(x)
    else if(prod(x) >= 0){
      if(min(abs(x)) < thr)
        sgn * unique(c(pretty_breaks()(c(min(abs(x)), thr)),
                       log_breaks(base)(c(max(abs(x)), thr))))
      else
        sgn * log_breaks(base)(sgn * x)
    } else {
      if(min(abs(x)) < thr)
        unique(c(sgn * log_breaks()(c(max(abs(x)), thr)),
                 pretty_breaks()(c(sgn * thr, x[which.min(abs(x))]))))
      else
        unique(c(-log_breaks(base)(c(thr, -x[1])),
                 pretty_breaks()(c(-thr, thr)),
                 log_breaks(base)(c(thr, x[2]))))
    }
  }
  trans_new(paste("symlog", thr, base, scale, sep = "-"), trans, inv, breaks)
}

I am not sure whether the impact of a parameter scale is the same as in Python, but here are a couple of comparisons (see Python version here):

data <- data.frame(x = seq(-50, 50, 0.01), y = seq(0, 100, 0.01))
data$y2 <- sin(data$x / 3)
# symlogx
ggplot(data, aes(x, y)) + geom_line() + theme_bw() +
  scale_x_continuous(trans = symlog_trans())

enter image description here

# symlogy
ggplot(data, aes(y, x)) + geom_line() + theme_bw()
  scale_y_continuous(trans="symlog")

enter image description here

# symlog both, threshold = 0.015 for y
# not too pretty because of too many breaks in short interval
ggplot(data, aes(x, y2)) + geom_line() + theme_bw()
  scale_y_continuous(trans=symlog_trans(thr = 0.015)) + 
  scale_x_continuous(trans = "symlog")

enter image description here

# Again symlog both, threshold = 0.15 for y
ggplot(data, aes(x, y2)) + geom_line() + theme_bw()
  scale_y_continuous(trans=symlog_trans(thr = 0.15)) + 
  scale_x_continuous(trans = "symlog")

enter image description here

Julius Vainora
  • 47,421
  • 9
  • 90
  • 102