ggplot scale transformation inaccurate for stat_function

Question

I have some right-skewed data, and I'd like to visually compare distribution fits to the data on the regular scale and the log-scale using ggplot. However, when I use the scale_x_continuous() or scale_x_log10() to convert the distribution curves, the transformation does not translate correctly.

x <- rlnorm(1000, meanlog = -4, sdlog = 1)
ggplot(data.frame(x)) +
  geom_histogram(aes(x, y = ..density.. * 25)) +
  scale_x_log10() +
  stat_function(fun = "dlnorm",
                args = list(meanlog = -4,
                            sdlog = 1))

Notice how the mean of the lognormal curve does not match the mean of the histogram. Why not? Is there a way to get them to match?

In a different but related post, a suggested answer was to include the argument inherit.aes = FALSE, but that does not help here.

I am using R version 3.4.3 and ggplot2 version 2.2.1.

Answer 1

First, when working with the log normal distribution recall that the default is to work with the natural logarithm, not the base 10 logarithm. Part of the issue with the graphic above is due to the mixing of logarithm bases.

Let's work through this by first generating the example observations of a log normal random variable with meanlog -4 and sdlog 1, that is,

library(ggplot2)
library(gridExtra)

set.seed(42)

dat <- data.frame(x = rlnorm(1000, meanlog = -4, sdlog = 1))

We will start by plotting the density on the standard x-axis. I'll use geom_histogram with stat = "density" so that the bars are scaled and there is no need to use the aesthetic y = ..density.. This is very similar to your original plot, just no attempt to scale the x axis.

ggplot(dat) +
  geom_histogram(mapping = aes(x = x), stat = "density")  +
  stat_function(fun = "dlnorm",
                args = list(meanlog = -4, sdlog = 1),
                n = 501,
                color = "red")

Now, recall that if

then

where the log is the natural logarithm.

One way to plot the generated data example on the log scale is as follows. Note that the log transform is explicit in the mapping for geom_historgram and that the stat_function is using dnorm not dlnorm.

ggplot(dat) +
  geom_histogram(mapping = aes(x = log(x)), stat = "density")  +
  stat_function(fun = "dnorm",
                args = list(mean = -4, sd = 1),
                n = 501,
                color = "red")

Now, to transform the x-axis you will want to use ggplot2::scale_x_continuous with the trans = "log" argument. When this transform is applied to the graphic, the scale of the x-axis is modified and the evaluation of the stat_function will occur on the transformed x values, not the original values. Thus, you'll need to define the function to use dnorm(log(x)) as shown below:

ggplot(dat) +
  geom_histogram(mapping = aes(x = x), stat = "density")  +
  stat_function(fun = function(x, ...) {dnorm(log(x), ...) },
                args = list(mean = -4, sd = 1),
                n = 501,
                color = "red") +
  scale_x_continuous(trans = "log",
                     breaks = exp(seq(-6, 0, by = 2)),
                     labels = paste("exp(", seq(-6, 0, by = 2), ")"))

It is worth noting that that labels for the x-axis ticks in the second plot are integer values and the x-axis label is log(x) whereas in the third graphic the x-axis ticks are expressions and the label is plan "x." Make sure you are using descriptive tick and axis labels.

Answer 2

The goal (if not stated so originally) is still to view the lognormal data and distribution on the log10 scale. In order to reach this goal, the density (the pdf) needs to be for the log10-scale. (Thank you to a colleague who shared the code below!)

## generate data:
x <- rlnorm(1000, meanlog = -4, sdlog = 1)

## generate sequence of x values for the curve.
xx <- seq(min(x), max(x), length = 1000)
## Calculated the density for each xx value.
## Here, density is based on the lognormal distribution.
pdf <- dlnorm(xx, -4, 1)

## Repeat for log(xx).
xx_ln <- log(xx)
## This density is based on the normal distribution.
pdf_norm <- dnorm(xx_ln, -4, 1)

## As a reminder, the pdf's for the distributions are different:
head(cbind(pdf, pdf_norm))

When looking at the data on a log10-scale, it will also have a different pdf. The function and code below transforms the normal pdf into a pdf for the log10-scale.

## Function: numerical integration stuff for log10 distribution plots
## essentially transforms pdf_norm to log10 base.
## step_size = Riemann sum-- step size to integrate over.
## x_10 = x values after a log10-transformation
## pdf_norm == pdf values for normal distribution (see above)
num_int <- function(df){
  df$step_size <- c(diff(df$xx_10), NA)
  int <- sum(df$step_size * df$pdf_norm, na.rm = T)
  return(data.frame(int))
}

## to complete the numerical integration, need log10(values)
xx_10 <- log10(xx)
curve_df <- data.frame(xx, xx_10, pdf, pdf_norm)
int <- num_int(curve_df) 
curve_df$pdf_10 <- curve_df$pdf_norm / as.numeric(int)

## replace Inf rows with NA
## (not necessary with the example code)
curve_df %<>%
  mutate(pdf = replace(pdf, pdf == Inf, NA),
         pdf_norm = replace(pdf_norm, pdf_norm == Inf, NA),
         pdf_10 = replace(pdf_10, pdf_10 == Inf, NA))


ggplot() +
  geom_histogram(data = data.frame(x), aes(x = x, y = ..density..)) + 
  geom_line(data = curve_df,
            aes(xx, pdf_10), col="blue", size = I(1.2), linetype = 1) +
  scale_x_log10()