stat_theodensity with scale_x_log10

[FelixTheStudent]

Cool package, love it!

I found that scale_x_log10 does not work together with stat_theodensity, but this would be a nice feature for log-normal data.

Here's a minimal example:

library(tidyverse)
library(ggh4x)

x <- data.frame(x=rlnorm(1000, 1,1))
p<- x %>% ggplot(aes(x))+geom_histogram(bins=100)+
  stat_theodensity(aes(y = stat(count)),
                   distri = "lnorm")
p  # works
p + scale_x_log10() # error
p + coord_trans(x="log10") # error

The error produced is:

Error in computing default starting values. Warning message: Computation failed in stat_theo_density(): Error in start.arg.default(obs, distname) : values must be positive to fit a lognormal distribution

This is with R version 3.6.3. Not sure if this is a bug or if what I'm asking is unreasonable (scales and coords are a tricky concept).

Cheers,

Felix

[teunbrand]

Hi, thank you for bringing this to my attention.

I can indeed confirm that the scale_x_log10() produces an error. This happens for understandable reasons; scale transformations occur prior to the calculation of any statistic, so the statistic is calculated on the transformed data. The log10 transformation makes the lognormal distribution a normal distribution with negative values, which the fitdistrplus::fitdist() function (that estimates the parameterisation under the hood) can't handle for lognormal distributions.

For this I can see 2 possible solutions:

• Simply use the normal distribution if you use the log10 scale; but this work only for normal-lognormal scale/stat pairs.
• I should make an option to untransform the data prior to fitting the data and then re-transform the result.

With regards to the coordinate transformation, I don't think it is stat_theodensity() that is causing the problems, but the histogram. If we compare the following 3 plots, the histogram fails when the binwidth is wide.

df <- data.frame(x = rlnorm(1000, 1, 1))

ggplot(df, aes(x)) +
  geom_histogram(bins = 100) +
  coord_trans(x = "log10") # Error

ggplot(df, aes(x)) +
  geom_histogram(binwidth = 0.1) +
  coord_trans(x = "log10") # Works

ggplot(df, aes(x)) +
  stat_theodensity(distri = "lnorm") +
  coord_trans(x = "log10") # Works

Because histograms are parameterised under the hood as rectangles, if the xmin is below zero, that would produce an error. Maybe this is an issue that might be raised at ggplot2 itself? The solution I'd put forward is to choose a smaller binwidth for now.

As a small aside, the overlay of a histogram and a (scaled) density curve works best if the stat(count) is multiplied by the width of a histogram's bin (which is variable when you set bins = 100). I recommend the following for your example:

ggplot(df, aes(x)) +
  geom_histogram(binwidth = 0.1) +
  stat_theodensity(aes(y = stat(count * 0.1)), 
                   distri = "lnorm") +
  coord_trans(x = "log10")

Do you think the untransform option would solve the issue?

Best wishes

[FelixTheStudent]

Amazing answer, thank you. I like the coord_trans option a lot!I found two more interesting tweaks that I'll leave here as a note to others and probably my future self.

x <- data.frame(x=rlnorm(1000, 1,1))

# stat(density) achieves the same as stat(count*.1) for obvious reasons
ggplot(x, aes(x,stat(density)))+
  geom_histogram(binwidth = 0.1 )+
  stat_theodensity(distri = "lnorm") +
  coord_trans(x = "log10")

# manually supplying breaks solves the xmin issue. 
# If desired, one can make them logarithmic so they appear to have
# equal width in final plot:
logseq <- function(x, length = 100, abs_min=1e-5){
  # exponentially spaces sequence covering the range of x (+/- ten percent)
  pmax(abs_min, exp(seq(log(.9*min(x)), log(1.1*max(x)), length.out = length)))
}
ggplot(x, aes(x,stat(density)))+
  geom_histogram(breaks = logseq(x) )+
  stat_theodensity(distri = "lnorm") +
  coord_trans(x = "log10")

This solves the issue, as far as I'm concerned. Thanks again for the nice work and helpful package!