stat_function的ggplot比例转换不准确

Question

我有一些右偏数据，我想使用 ggplot 直观地比较分布拟合与常规比例和对数比例的数据。但是，当我使用 scale_x_continuous() 或 scale_x_log10() 转换分布曲线时，转换无法正确转换。

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

注意对数正态曲线的平均值与直方图的平均值不匹配。为什么不？有没有办法让它们匹配？

在另一篇不同但相关的帖子中，建议的答案是包含参数inherit.aes = FALSE，但这在这里没有帮助。

我使用的是 R 版本 3.4.3 和 ggplot2 版本 2.2.1。

Answer 1

首先，当使用对数正态分布时，请记住默认是使用 自然 对数，不是 以 10 为底的对数.上图的部分问题是由于对数底数的混合。

让我们首先生成对数正态随机变量 的示例观测值，meanlog -4 和 sdlog 1，即

library(ggplot2)
library(gridExtra)

set.seed(42)

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

我们将从在标准 x 轴上绘制密度开始。我将使用geom_histogram 和stat = "density" 以便缩放条形并且无需使用美学y = ..density.. 这与您的原始图非常相似，只是没有尝试缩放x 轴。

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

现在，回想一下，如果

然后

log 是自然对数。

在对数刻度上绘制生成的数据示例的一种方法如下。请注意，日志转换在 geom_historgram 的映射中是显式的，并且 stat_function 使用的是 dnorm 不是 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")

现在，要转换 x 轴，您需要使用 ggplot2::scale_x_continuous 和 trans = "log" 参数。当将此变换应用于图形时，x 轴的比例会被修改，stat_function 的评估将发生在变换后的 x 值上，而不是原始值上。因此，您需要定义函数以使用dnorm(log(x))，如下所示：

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), ")"))

值得注意的是，第二张图中 x 轴刻度的标签是整数值，x 轴标签是 log(x)，而在第三张图中，x 轴刻度是表达式，标签是计划“x”。确保您使用的是描述性刻度和轴标签。

Answer 2

目标（如果最初没有这样说明）仍然是查看 log10 尺度上的对数正态数据和分布。为了达到这个目标，密度（pdf）需要为 log10 尺度。 （感谢分享以下代码的同事！）

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

查看 log10 尺度的数据时，它也会有不同的 pdf。下面的函数和代码将普通 pdf 转换为 log10-scale 的 pdf。

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