如何使用ggplot2绘制功率曲线

Question

I希望使用ggplot2来可视化遵循幂函数曲线的一些数据。这个问题以前已经被问过了（Add exp/power trend line to a ggplot），但是答案并没有真正帮助到我。

一个技巧是使用stat_function()来创建曲线。然而，我无法使stat_function()和我的幂函数曲线在对数刻度下工作。

我说明我的问题。

创建一些示例数据和基本图：

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)

很好，让我们使用coord_trans()添加一些对数刻度。这完美地起作用，只是我的直线不再是直的了（正如文档告诉我所期望的那样）。

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

因此，请再次尝试使用coord_x_log10()和coord_y_log10()，但是这会引发错误：

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

这可能与我需要调整函数以反转比例尺的效果有关，但我还无法完全理解。

只有对数x轴比例尺时，我才能做到：

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

如果我添加了scale_y_log10()，我不知道如何转换y值。可能是我缺少一些基本的东西。有没有简单的方法来绘制这条曲线？

- Andrie

---

如果ggplot2这个部分不是决定性因素，那么有一个叫做poweRlaw的包可以实现这个功能：https://cran.r-project.org/web/packages/poweRlaw/index.html。 - ulfelder

两个问题：您需要设置轴限制以避免零/负值。例如，scale_x_log10(limits = c(1, 200)) + scale_y_log10(limits = c(1e-5, 1))。然而，stat_function与轴变换不兼容：https://dev59.com/CWDVa4cB1Zd3GeqPfrHm - Roland

@Roland 我也尝试手动设置限制，但是如果你不知道数据的确切情况，那么这样做通常有多大用处呢... - joran

@BrandonBertelsen 我是故意这样做的。geom_point() 使用数据框绘制点。stat_function() 不会触及底层数据，而只使用函数。我认为它们不会相互干扰。 - Andrie

ggplot(dat, aes(x = x, y = y)) + geom_point() + coord_trans(ytrans = "log10", xtrans = "log10")...这不是你想要的吗？ - Brandon Bertelsen

Answer 1

将我的评论转化为答案：

主要问题在于stat_function中存在一个bug。当与轴变换一起使用时，它会从经过转换的而不是原始值计算y值。这个问题最近已经在github上得到了修复。

然而，这并不容易看出，因为第一个错误发生在计算轴断点时，由于该bug产生了Inf、零和/或负的y值。您需要设置明确的轴限制才能看到实际的问题是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).

如果只转换x轴，则更加明显：

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

如果您无法使用来自github的ggplot2版本，您可以使用以下内容：

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()

- Roland