mixed model for 'cursed' Morgan Stanley data

morganstanley_cursed_mixed.R

## setup from  https://www.tjmahr.com/morgan-stanley-cursed-covid-plot/

library(tidyverse)
theme_set(theme_bw())
library(lme4)
library(ggeffects)
library(colorspace)
library(directlabels)

# a helper function to download the data from github
# in case you want to play along
path_blog_data <- function(x) {
  file.path(
    "https://raw.githubusercontent.com",
    "tjmahr/tjmahr.github.io/master/_R/data",
    x
  )
}

data <- readr::read_csv(
  path_blog_data("2020-05-12-states_daily_4pm_et.csv"),
  col_types = cols(
    date = col_date("%Y%m%d"),
    state = col_character(),
    inIcuCurrently = col_number(),
    .default = col_skip()
  ),
  progress = FALSE
)

data <- data %>%
  filter(
    as.Date("2020-04-28") <= date,
    date <= as.Date("2020-05-11")
  )


## new stuff starts here

## convenient to fit to an integer date starting from zero
ddm <- data %>% mutate(ndate = as.numeric(date-min(date)))

## fit random-slopes model
m1 <- lmer(inIcuCurrently ~ ndate + (ndate|state), ddm,
           na.action = na.exclude)

## add state-by-state predicted values
ddm <- (mutate(ddm, pred = predict(m1))
  %>% group_by(state)
  %>% mutate(nstate = ifelse(
                 all(is.na(inIcuCurrently)) |
                 min(inIcuCurrently, na.rm = TRUE) < 350,
             "", state))
  %>% ungroup()
)


## population-level predicted values + CIs
## (ignores uncertainty in RE)
d0 <- ggpredict(m1)$ndate %>% full_join(select(ddm, ndate, date),
                                        by = c("x" = "ndate"))

gg0 <- (ggplot(ddm, aes(date, inIcuCurrently, colour = state))
  + geom_point()
  + geom_line(aes(y=pred), alpha = 0.2)
  + geom_line(data=d0, aes(y=predicted), lwd=2, colour = "black",
              alpha = 0.7)
  + geom_ribbon(data=d0, aes(y = predicted,
                             ymin=conf.low, ymax = conf.high),
                colour = NA, alpha = 0.2)
  + scale_color_discrete_qualitative(guide=guide_none())
)

gg0L <- gg0 + geom_dl(aes(label = nstate), method =
                                    list(dl.trans(x = x+0.2),
                                         "last.bumpup"))

print(gg0L)
ggsave(gg0L, file = "morganstanley.png")