Как подогнать отрицательную биномиальную функцию к данным, которые facet_wrapped фактор в R ggplot? - PullRequest
Новая
       69

Как подогнать отрицательную биномиальную функцию к данным, которые facet_wrapped фактор в R ggplot?

1 голос
/ 15 апреля 2020

Я пытаюсь подогнать отрицательное биномиальное распределение для подсчета данных, но сократил до подсчета, как в в этом примере В моих данных я должен разделить графики биномиальной функции для двух видов. Тем не менее, это не простой способ указать это в функции и получить легенды строки со значениями параметров в ключе для обоих видов. 

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)

enter image description here

1 Ответ

2 голосов
/ 15 апреля 2020

Вы правы, это немного больно, чтобы получить права. Я немного адаптировал ваш пример, чтобы более четко показать два разных дистрибутива. Вот моя попытка заставить ваш подход работать: 

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

Создано в 2020-04-15 пакетом Представить (v0 .3.0)

Существуют также пакеты, которые выполняют большую часть этой работы за вас. Отказ от ответственности: я написал ggh4x, поэтому я не беспристрастный. Вы также можете заменить код ggplot следующим (при условии аналогичной предварительной обработки) 

library(ggh4x) # devtools::install_github("teunbrand/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)

Надеюсь, что это помогло!

...