In databrew/covid19: COVID-19 App

# To deploy scp -i ~/.ssh/odkkey.pem uci.Rmd [e-mail]:/home/ubuntu/Documents
# Basic knitr options
library(knitr)
opts_chunk$set(comment = NA, 
               # echo = FALSE, 
               warning = FALSE, 
               message = FALSE, 
               error = TRUE, 
               cache = FALSE,
               tidy = TRUE,
              # tidy.opts=list(width.cutoff=60),
              # tidy=TRUE,
               # fig.width = 8.64,
               # fig.height = 4.86,
               fig.path = 'figures/')

## Load libraries
library(covid19) #devtools::install_github('databrew/covid19')
library(ggplot2)
library(lubridate)
library(dplyr)
library(ggplot2)
library(shiny)
library(ggthemes)

# Define parameters
# Basing off: https://www.niusdiario.es/multimedia/nius-te-explica/colpaso-sistema-sanitario-espana-uci-unidad-cuidados-intensivos-hospitales-coronavirus_18_2914020237.html
p <- list(
  # beds_public = 3508,
  # beds_private = 896,
  normal_occupancy = 25,
  # need_hospitalization = 15, # percent which requires hospitalization
  # need_uci = 5, # percent which requires uci
  beds = 5799, # https://www.infolibre.es/noticias/politica/2020/03/21/uci_termometro_coronavirus_covid_19_105153_1012.html
  average_days_in_uci = 10 # total rough estimate
)
options(scipen = '999')


make_prediction <- function(data,
                            n_start = 5,
                            cumulative = FALSE,
                            time_ahead = 7,
                            var = 'uci'){

  sub_data <- data
  # Get the var
  the_var <- paste0(var, ifelse(cumulative, '', '_non_cum'))
  sub_data$var <- as.numeric(unlist(sub_data[,the_var]))
  # narrow down
  sub_data <- sub_data %>%
    dplyr::select(date, var)

  pd <- sub_data %>%
        filter(!is.na(var)) %>%
      mutate(start_date = min(date[var >= n_start])) %>%
      mutate(days_since = date - start_date) %>%
      filter(days_since >= 0) %>%
      mutate(days_since = as.numeric(days_since)) %>%
      mutate(the_weight = 1/(1 + (as.numeric(max(date) - date))))
  fit <- lm(log(var) ~ days_since,
              weights = the_weight,
              data = pd) 


    # Predict days ahead
    day0 <- pd$date[pd$days_since == 0]
    fake <- tibble(days_since = seq(0, max(pd$days_since) + time_ahead, by = 1))
    fake <- fake %>%mutate(date = seq(day0, day0+max(fake$days_since), by = 1))
    fake <- left_join(fake, pd %>% dplyr::select(days_since, var, date))
    fake$predicted <- exp(predict(fit, newdata = fake))
    # fake$predictedlo <- predict(fitlo, newdata = fake)
    ci <- exp(predict(fit, newdata = fake, interval = 'prediction'))
    # cilo <- predict(fitlo, newdata = fake, interval = 'prediction')

    fake$lwr <- ci[,'lwr']
    fake$upr <- ci[,'upr']
    # fake$lwrlo <- ci[,'lwr']
    # fake$uprlo <- ci[,'upr']
    # Doubling time
    dt <- log(2)/fit$coef[2]
    fake %>% mutate(doubling_time = dt)
}

Column {.sidebar}

When will Spain reach its ICU capacity ceiling?

sliderInput(
  "normal_occupancy", label = "Baseline (non COVID-19) occupancy (%)",
  min = 0, max = 100, value = 25
)

sliderInput(
  "beds", label = "UCI beds:",
  min = 4000, max = 10000, value = 5800, step = 100,
)

sliderInput(
  'average_days_in_uci',
  'Number of days in UCI for a COVID-19 patient',
  min = 1, max = 20, step = 1, value = 10
)

sliderInput(
  'time_ahead',
  'Number of days to predict into the future',
  min = 1, max = 30, step = 1, value = 10
)

sliderInput(
  'n_start',
  'The model will be based on the period beginning at which there were more than the following number of COVID-19 UCI admissions',
  min = 10, max = 200, step = 5, value = 10
)

Column

renderPlot({

  n_start = input$n_start
  time_ahead = input$time_ahead

spain_data <- 
  esp_df %>% group_by(date) %>%
                               summarise_at(.vars = vars(uci, deaths, confirmed_cases,
                                            uci_non_cum,
                                            deaths_non_cum,
                                            confirmed_cases_non_cum),
                                            .fun = function(x){sum(x, na.rm = TRUE)})

pd <- make_prediction(data = spain_data,
                      n_start = n_start,
                      cumulative = FALSE,
                      time_ahead = time_ahead,
                      var = 'uci')
# Get the number of daily admissions to "spillover" for the number of days they need to be in UCI
# this function not estimating confidence bounds at this point
spill_over <- function(data,
                       days = p$average_days_in_uci){
  out <- data
  out$predicted_spilled_over <- NA
  out_list <- list()
  for(i in 1:nrow(out)){
    # Get the sub data for up to days days before
    sub_data <- out %>%
      filter(date >= out$date[i] -( days-1),
             date <= out$date[i])
    # Get the sum of ingresado people during that window
    sum_ingresado <- sum(sub_data$var, na.rm = T)
    # Get the predicted sum too
    sum_predicted <- mean(sub_data$predicted, na.rm = TRUE) * nrow(sub_data)
    # Manually replace with observed
    out_predicted <- ifelse(!is.na(sub_data$var[i]),
                            sum_ingresado,
                            sum_predicted)
    message(i, ": ", round(out_predicted))
    # pop back into dataframe
    # out_list[[i]] <- out_predicted
    out$predicted_spilled_over[i] <- out_predicted
  }
  # out <- unlist(out_list)
  return(out)
}
# preds <- spill_over(pd)
average_days_in_uci <- input$average_days_in_uci
normal_occupancy <- input$normal_occupancy
beds <- input$beds

pd <- pd %>% spill_over(days = average_days_in_uci)

already_occupied <- (normal_occupancy/100) * beds
total_beds <- beds

# Shape data for plotting
plot_data <- pd %>%
  dplyr::select(date, var, predicted, predicted_spilled_over) %>%
  mutate(already_occupied = already_occupied) %>%
  tidyr::gather(key, value, var:already_occupied) %>%
  mutate(key = ifelse(key == 'already_occupied',
                      'UCI beds occupied (normal)',
                      ifelse(key == 'predicted_spilled_over',
                             'Predicted COVID-19 UCI beds needed',
                             key)))

basic_data <- plot_data %>%
             filter(key %in% c('UCI beds occupied (normal)',
                               'Predicted COVID-19 UCI beds needed'))

ggplot() +
  geom_bar(data = basic_data,
           aes(x = date,
               y = value,
               fill = key),
           stat = 'identity',
           position = position_stack(),
           alpha = 0.7, width = 1,
           color = 'white',
           lwd = 0.1) +
  geom_text(data = basic_data %>% filter(key == 'Predicted COVID-19 UCI beds needed'),
           aes(x = date,
               y = value + already_occupied * 1.1,
               label = round(value, digits = 0)),
           # stat = 'identity',
           position = position_stack(),
           alpha = 0.7) +
  # scale_y_log10() +
  geom_hline(yintercept = total_beds,
             lty = 2) +
  geom_text(data = tibble(x = min(plot_data$date) +2,
                          y = total_beds + 300,
                          label = 'Total UCI beds'),
            aes(x = x,
                y = y,
                label = label)) +
  geom_text(data = tibble(x = min(plot_data$date) +2,
                          y = already_occupied - 500,
                          label = paste0('Regular UCI occupancy:\n', round(already_occupied), ' (', p$normal_occupancy, '%)')),
            aes(x = x,
                y = y,
                label = label)) +
  theme_simple() +
  theme(legend.position = 'top') +
  scale_fill_manual(name = '',
                     values = c('red', 'black')) +
  labs(x = 'Date',
       y = 'UCI beds',
       title = 'UCI bed capacity, Spain',
       subtitle = paste0('Based on simple log-linear growth model in daily UCI admissions.\nAssumes an average of ', p$average_days_in_uci, ' days in ICU.'),
       caption = paste0('Code at https://github.com/databrew/covid19/blob/master/misc/uci/uci.Rmd\nAssumes: ', total_beds, ' total beds, and a ', p$normal_occupancy, ' non-COVID19 occupancy rate.')) +
  scale_x_date(breaks = sort(unique(plot_data$date)),
               labels = format(sort(unique(plot_data$date)), '%d\n%b'))
},
height = 700)

Model specification pseudo-code

model = log(daily UCI admission) ~ days, weights = 1 / (1 + today - date)

Who made this?

• [e-mail]
• [e-mail]
• [e-mail]

databrew/covid19 documentation built on Aug. 24, 2020, 10:39 a.m.