如何通过 R ggplot 中的一个因素将负二项式函数拟合到 facet_wrapped 的数据?

How to fit negative binomial function to data that is facet_wrapped by a factor in R ggplot?

我正在尝试将负二项分布拟合到计数数据,但按比例缩小到 在我的数据中,我必须分离出两个物种的二项式函数绘图。但是,没有简单的方法来在函数中指定它并在两个物种的键中获取带有参数值的线图例。

set.seed(111)
count <- rbinom(500,100,0.1) 
species <- rep(c("A","B"),time = 250)
df <- data.frame(count,species)

#Specifying negative binomial function
negbinom.params <- fitdistr(df$count,"negative binomial", method = "SANN")$estimate
dist.params <- map(list(`Negative Binomial` = negbinom.params),~ map2(names(.),.,~ paste0(.x," = ",round(.y,2))) %>% unlist %>% paste0(.,collapse = ", ")) %>% map2_chr(names(.),., ~ paste(.x,.y,sep=":\n"))

#Plotting
mybinwidth = 2
ggplot(df, aes(x = count, colour = species, fill = species)) + 
  facet_grid(.~species) +  
  geom_histogram(aes(y=..count..),alpha = 0.5, binwidth = mybinwidth) + 
  stat_function(aes(color = "orange"), 
                fun = function(x,size, mu) {
                    mybinwidth * nrow(df) * dnbinom(x,size = size, mu = mu)
                },
                args=fitdistr(df$count, "negative binomial", method="SANN")$estimate, 
                xlim=c(0,50),n=20)

你说得对,要做到这一点有点痛苦。我对您的示例进行了一些调整,以更清楚地显示两种不同的分布。这是我尝试使您的方法起作用的尝试:

library(ggplot2)
library(MASS)
#> Warning: package 'MASS' was built under R version 3.6.2

set.seed(111)

df <- data.frame(
  count = rnbinom(500, rep(c(5, 10), each  = 250), 0.5),
  species = rep(c("A", 'B'), each = 250)
)

# Not the prettiest formatting, but it'll show the point
ests <- sapply(split(df$count, df$species), function(x) {
  est <- fitdistr(x, "negative binomial", method = "SANN")$estimate
  formatted <- paste0(names(est)[1], " = ", format(est, digits = 2)[1], ",",
                      names(est)[2], " = ", format(est, digits = 2)[2])
  formatted
})

mybinwidth <- 1

spec_A = df[df$species == "A",]
spec_B = df[df$species == "B",]

ggplot(df, aes(count)) +
  geom_histogram(binwidth = mybinwidth,
                 aes(fill = species), alpha = 0.5,
                 position = "identity") +
  stat_function(aes(color = "A"), 
                data = data.frame(species = "A"),
                fun = function(x, size, mu) {
                  mybinwidth * nrow(spec_A) * dnbinom(x,size = size, mu = mu)
                },
                args = fitdistr(spec_A$count, "negative binomial", method="SANN")$estimate, 
                xlim = c(0, max(df$count)), n= max(df$count) + 1, inherit.aes = FALSE) +
  stat_function(aes(color = "B"), 
                data = data.frame(species = "B"),
                fun = function(x, size, mu) {
                  mybinwidth * nrow(spec_B) * dnbinom(x,size = size, mu = mu)
                },
                args = fitdistr(spec_B$count, "negative binomial", method="SANN")$estimate, 
                xlim = c(0, max(df$count)), n= max(df$count) + 1, inherit.aes = FALSE) +
  scale_colour_discrete(labels = unname(ests), name = "fit") +
  facet_wrap(~ species)
#> Warning: `mapping` is not used by stat_function()
#> Warning: `data` is not used by stat_function()
#> Warning: `mapping` is not used by stat_function()
#> Warning: `data` is not used by stat_function()

reprex package (v0.3.0)

于 2020-04-15 创建

还有一些包可以为您完成大部分工作。免责声明:我写了 ggh4x,所以我不是没有偏见的。您也可以将ggplot代码替换为以下内容(假设预处理类似)

library(ggh4x)
ggplot(df, aes(count)) +
  geom_histogram(binwidth = mybinwidth,
                 aes(fill = species), alpha = 0.5,
                 position = "identity") +
  stat_theodensity(aes(colour = species,
                       y = after_stat(count * mybinwidth)),
                   distri = "nbinom") +
  scale_colour_discrete(labels = unname(ests), name = "fit") +
  facet_wrap(~ species)

希望对您有所帮助!