时间戳

首页 > 解决方案 > 如何将负二项式函数拟合到由 R ggplot 中的一个因子包裹的 facet_wrap 的数据？

r - 如何将负二项式函数拟合到由 R ggplot 中的一个因子包裹的 facet_wrap 的数据？

问题描述

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

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)

在此处输入图像描述

标签： rggplot2distributiondata-fitting

解决方案

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

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 包于 2020-04-15 创建(v0.3.0)

还有一些软件包可以为您完成大部分工作。免责声明：我写了 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)

希望有帮助！

推荐阅读