How to draw a power curve using ggplot2

Question

I want to use ggplot2 to visualize some data that follows a power curve. This has been asked before (Add exp/power trend line to a ggplot) but the answer didn't really help.

One trick is to use stat_function() to create a curve. However, I am unable to get stat_function() and my power curve to work with logarithmic scales.

I illustrate my problems.

Create some sample data and a base plot:

library(ggplot2)

x <- 1:100
pwr <- function(x)x^-2.5
dat <- data.frame(x, y = pwr(x))

p <- ggplot(dat, aes(x = x, y = y)) + 
  geom_point()

p + stat_function(fun = pwr)

Great, let's add some logaritmic scales using a coord_trans(). This works perfectly, except that my straight lines are no longer straight (exactly as the documentation tells me to expect).

p + stat_function(fun = pwr) + coord_trans(x = "log10", y = "log10")

So, try again with coord_x_log10() and coord_y_log10(), but this throws an error:

p + stat_function(fun = pwr) + scale_x_log10() + scale_y_log10()

Error in seq.default(min, max, by = by) : 
  'from' cannot be NA, NaN or infinite

This has most likely to do with the fact that I have to adjust my function to invert the effect of the scales, but I can't quite figure it out.

I can do it with only the logarithmix x-scale:

p + scale_x_log10() + stat_function(fun = function(x)pwr(10^x))

I have no idea how to convert the y-values if I add scale_y_log10().

I'm probably missing something fundamental. Is there an easy way to plot this curve?

Answer 1

Putting my comments into an answer:

The main problem is a bug in stat_function. When used in conjunction with axis transformations it calculates y-values from the transformed and not from the original values. This has been fixed on github very recently.

However, this is not easy to see since the first error occurs when calculating the axis breaks since the bug produces Inf, zero and/or negative y-values. You need to set explicit axis limits to see that the actual problem is stat_function:

p + stat_function(fun = pwr) + 
  scale_x_log10() + scale_y_log10(limits = c(1e-5, 1))
#Warning message:
#Removed 100 rows containing missing values (geom_path).

It's even more obvious if you only transform the x-axis:

p + stat_function(fun = pwr) + 
  scale_x_log10() 

If you can't use the ggplot2 version from github you could use this:

p + geom_line(data = data.frame(x = seq(min(dat$x), max(dat$x), length.out = 100)),
              aes(y = pwr(x))) + 
  scale_x_log10() + scale_y_log10()