Блог

Я использую скользящие средние, чтобы сгладить влияние дня недели на распределение вакцины, чтобы увидеть общие тенденции, стратифицированные различными факторами. Я могу создать гистограмму скользящих средних, которая правильно отображает общие данные. Но когда я стратифицирую или создаю фасеты, во вводном периоде появляются «призрачные» полосы убывающей высоты (на которых не должно быть полос). Как я могу этого избежать?

правильный график (без стратификации): g

графики с «призрачными» столбиками в периоде ввода скользящей средней: g facet_grid(race ~., scales=»free_y»)

мой код

 library(tidyverse)
# Make fake data: count of doses per day for 70 days, increasing over the 70 days, with a 50% variance per day-of-week
nPerDay <- floor(sample(5:10, 70, replace=T) * (1   ((1:70)*3/70)) * (.5   (.5*(1:70 %% 7)/6)))
# Use that to create a data frame where one record is the administration of one dose, giving the dose, vaccine brand, 1st or 2nd dose, pt race, amp; pt gender
doses <- data.frame(Admin_date = rep(as.Date("2020-12-31")   1:70, nPerDay)
                    , whichDose = factor(c(rep(1,sum(nPerDay[1:30])), sample(1:2, sum(nPerDay[31:70]), replace=T)))
                    , gender=sample(c("F", "M"), sum(nPerDay), replace=T)
                    , race=sample(LETTERS[1:5], sum(nPerDay), c(.45, .25, .15, .1, .05), replace=T)
                    , brand=sample(c("Pf", "Mo"), sum(nPerDay), replace=T)
)

# plot the doses administered each day, with stacked bars', with bars' color indicating # of 1st or second dose
(ggplot(data=doses, mapping=aes(x=Admin_date))#, fill=whichDose))
    geom_bar(position = "stack")
    geom_line(aes(y=..count.., fill=NULL), stat = "bin", binwidth=1)
)

# Change the bars in the prior plot into rolling 7-day averages, but keep the line as a daily total count.
g <- (
  ggplot(data=doses, mapping=aes(x=Admin_date))#, fill=whichDose)) 
    geom_bar(position = "stack"
             , mapping = aes(y=zoo::rollmean(..count.., 7, align="right", fill=NA))
             , stat="bin", binwidth=1
  )
    geom_line(aes(y=..count.., fill=NULL), stat = "bin", binwidth=1)
    labs(y="doses", fill="Which dose,n7d avg count")
)
g # display this base graph

# explore tha data
g   facet_grid(race~., scales="free_y") # See if the increasing trend and 1st vs 2nd dose porportions or similar across races.
 

Я знаю, что могу избежать этого, создав промежуточный фрейм данных, который предварительно вычислял скользящие средние для стратификации, которую я хочу. Но должен быть способ сделать это на лету в R, за https://yjunechoe .github.io/posts/2020-09-26-demystifying-stat-layers-ggplot2 / и, возможно, функция after_stat() . Но я не могу понять это. Я надеюсь найти простое решение, которым я мог бы поделиться со своей рабочей группой, чтобы другие (у которых меньше опыта в R) могли добавлять функции фасетов в базовый график, чтобы изучить множество факторов, которые у нас есть — у нас гораздо больше, чем пол, раса, бренд и whichDose. Если я смогу избавиться от призрачных полос, они могут просто добавить подобный код, чтобы получить другие стратификации:

 # look at other stratifications
g   facet_grid(gender, scales="free_y")
g   facet_grid(race~brand, scales="free_y")
g   facet_grid(race~gender, scales="free_y")
 

Проблема в том, что после вычисления статистики любое вычисление, которое происходит после статистики, не обязательно учитывает панели. Это создает проблемы zoo::rollmean , потому что он видит только один вектор значений. Следовательно, вам придется перебирать данные по панели.

 library(tidyverse)

nPerDay <- floor(sample(5:10, 70, replace=T) * 
                   (1   ((1:70)*3/70)) * (.5   (.5*(1:70 %% 7)/6)))
doses <- data.frame(
  Admin_date = rep(as.Date("2020-12-31")   1:70, nPerDay),
  whichDose = factor(c(rep(1,sum(nPerDay[1:30])), 
                       sample(1:2, sum(nPerDay[31:70]), replace=T))),
  gender=sample(c("F", "M"), sum(nPerDay), replace=T),
  race=sample(LETTERS[1:5], sum(nPerDay), c(.45, .25, .15, .1, .05), replace=T),
  brand=sample(c("Pf", "Mo"), sum(nPerDay), replace=T)
)


ggplot(data=doses[order(doses$race, doses$Admin_date), ], 
       mapping=aes(x=Admin_date))  
  geom_bar(position = "identity"
           , mapping = aes(y=after_stat(
             unlist(lapply(split(count, PANEL), zoo::rollmean, 
                           7, align = "right", fill = NA))
           ))
           , stat="bin", binwidth=1
  )  
  geom_line(aes(y=..count.., fill=NULL), stat = "bin", binwidth=1)   
  labs(y="doses", fill="Which dose,n7d avg count")  
  facet_grid(race ~ ., scales = "free_y")
#> Warning: Removed 30 rows containing missing values (geom_bar).
 

^{Создано 2021-02-20 пакетом reprex (версия v1.0.0)}

Аналогично подходу @teunbrand (который заслуживает похвалы за его краткое объяснение проблемы, к которой мне нечего добавить), но используя dplyr вспомогательную функцию и вспомогательную функцию, вы могли бы достичь желаемого результата следующим образом:

 library(tidyverse)

set.seed(42)

# Make fake data: count of doses per day for 70 days, increasing over the 70 days, with a 50% variance per day-of-week
nPerDay <- floor(sample(5:10, 70, replace = T) * (1   ((1:70) * 3 / 70)) * (.5   (.5 * (1:70 %% 7) / 6)))
# Use that to create a data frame where one record is the administration of one dose, giving the dose, vaccine brand, 1st or 2nd dose, pt race, amp; pt gender
doses <- data.frame(
  Admin_date = rep(as.Date("2020-12-31")   1:70, nPerDay),
  whichDose = factor(c(rep(1, sum(nPerDay[1:30])), sample(1:2, sum(nPerDay[31:70]), replace = T))),
  gender = sample(c("F", "M"), sum(nPerDay), replace = T),
  race = sample(LETTERS[1:5], sum(nPerDay), c(.45, .25, .15, .1, .05), replace = T),
  brand = sample(c("Pf", "Mo"), sum(nPerDay), replace = T)
)

my_rollmean <- function(count, group) {
  data.frame(group = group, count = count) %>% 
    group_by(group) %>% 
    mutate(roll = zoo::rollmean(count, 7, align = "right", fill = NA)) %>% 
    pull(roll)
}

# Change the bars in the prior plot into rolling 7-day averages, but keep the line as a daily total count.
g <- ggplot(data = doses, mapping = aes(x = Admin_date))  
  geom_bar(
    position = "stack",
    mapping = aes(y = my_rollmean(..count.., ..PANEL..)),
    stat = "bin", binwidth = 1
  )  
  geom_line(aes(y = ..count.., fill = NULL), stat = "bin", binwidth = 1)  
  labs(y = "doses", fill = "Which dose,n7d avg count")

# explore tha data
g   facet_grid(race ~ ., scales = "free_y") # See if the increasing trend and 1st vs 2nd dose porportions or similar across races.
#> Warning: Removed 30 rows containing missing values (position_stack).