Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to get something like Matplotlib's symlog scale in ggplot or lattice?

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))
}
like image 495
Brian B Avatar asked Jan 30 '13 20:01

Brian B


1 Answers

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

like image 119
Julius Vainora Avatar answered Oct 15 '22 15:10

Julius Vainora