COVID-19 in Montreal

September 19, 2020 - 4 minutes
Statistics R Shiny

This post is an analysis of the COVID-19 situation in Montreal demonstrated in Shiny Apps.

Historical data are downloaded from Frederic Harnois’s Github repository.

The new data can either downloaded from his Github repo or scraping from Santé Montréal. An example code is as follows,

library(rvest)
library(dplyr)
library(stringr)

url <- 'https://santemontreal.qc.ca/population/coronavirus-covid-19/situation-du-coronavirus-covid-19-a-montreal/'

file_name <- paste0('data/', Sys.Date() - 1, '.csv')

tbl <- read_html(url) %>% 
       html_nodes("table") %>%
       html_table() %>%
       .[[4]] %>%
       select(1:2) %>%
       rename('Arrondissements' = 1, 'Cas Confirmés' = 2) %>%
       mutate('Cas Confirmés' = str_replace_all(.[[2]], ' ', ''))

tbl[34, 1] <- 'Territoire à confirmer'

readr::write_csv(tbl, path = file_name)

The data contains the name of the sub-regions (Arrondissements) and confirmed cases, the file name is the date of the data. To show the weekly cases, we first have to manipulated the data. We first merge the data by date then subtract daily confirmed cases to the previous day, which will turn into daily increase cases. And finally, we aggregated them by weeks.

To obtained the the weekly adjusted cases to the population, we obtain the population from the report on Ville de Montréal about demography.

One thing may be notice that, after making subtractions, for some date, the increase is negative, I assume that this is some mis-reporting of the previous day. To simplify this, I just make negative cases into 0. But the real problem of negative cases should go further investigate.

library(dplyr)
library(tibble)
library(tidyr)


data_dir <- "mtl_data"      # path to folder that holds multiple .csv files
csv_files <- fs::dir_ls(data_dir, regexp = "\\.csv$")

data <- lapply(csv_files, readr::read_csv)

name <- tibble(data[[1]][,1]) %>% slice(-(n()-1), -n())

# extract data and convert to numeric
data <- purrr::map(data, 2) %>% 
        bind_cols() %>% 
        mutate_all( ~ stringr::str_replace(., "<", "")) %>%
        mutate_all(as.numeric) %>% 
        slice(-(n()-1), -n())

# convert cumulated data to daily increase
daily_inc <- t(apply(data, 1, 'diff'))
daily_inc[daily_inc < 0] <- 0
daily_inc <- cbind(data[,1], daily_inc)

date <- c(stringr::str_remove_all(csv_files, pattern = c('mtl_data/|.csv')))
week <- format(as.Date(date, format = '%Y-%m-%d'), format="%W")

# merge data
cvmtl_inc <- bind_cols(name, daily_inc)
colnames(cvmtl_inc) <- c('name', date)

# view by week
date_week <- tibble(date, week)

# map_idx
NUM <- c(24, 9, 71, 7, 27, 72, 11, 1, 10, 3,
         6, 17, 18, 23, 2, 74, 16, 75, 5, 13,
         22, 8, 19, 25, 15, 14, 76, 77, 21, 12,
         20, 26, 4)

pop <- c(134245, 42796, 3823, 19324, 166520, 32448, 48899,
         18980, 6973, 20151, 18413, 44489, 76853, 136024, 20276,
         3850, 84234, 5050, 23954, 69297, 104000, 31380, 106743, 
         139590, 98828, 78305, 4958, 921, 78151, 69229, 89170, 143853, 20312)

cvmtl_wk <- cvmtl_inc %>% 
  pivot_longer(starts_with('2020'), names_to = 'date', values_to = 'cases') %>%
  left_join(date_week) %>% 
  group_by(name, week) %>% 
  summarise(total = sum(cases)) %>%
  pivot_wider(names_from = 'week', values_from = 'total') %>% 
  add_column('30' = NA, .after = '29') %>% # no data for week 30
  add_column(NUM, .after = 'name') %>% 
  add_column(pop, .after = 'NUM') %>% 
  pivot_longer(4:ncol(.), names_to = 'week', values_to = 'case') %>% 
  mutate(adj = round(case/pop*10^5, 2))

In terms of the demonstration in Shiny Apps, there are three components, User Interface (UI), Server and shinyApp, for simplicity I put all three in one file. This may be improve by reactive function in the future.

library(tmap)
library(sf)
library(dplyr)

load('cvmtl_wk.rda')
mtl <- read_sf("mtl.shp")
week_range <- range(cvmtl_wk$week)
cvmtl <- mtl %>% select(NUM) %>% left_join(cvmtl_wk) %>% select(-NUM)

ui <- fluidPage(
  
  titlePanel("COVID-19 in Montreal by Weeks"),
  
  sidebarLayout(
    
    sidebarPanel(
      
      selectInput("select", label = h3("Legend Style"), 
                  choices = list("Fixed" = 'fixed', "Flexible" = 'pretty'), 
                  selected = 'fixed'),
      
      sliderInput(inputId = 'week', 
                  label = h3("Weeks"), 
                  min = week_range[1], 
                  max = week_range[2],
                  step = 1,
                  value = week_range[2], 
                  animate = TRUE),
      
      checkboxInput("checkbox", label = "Population Adjusted", value = FALSE),

      hr(),
      helpText(paste('Starting from lastest data availabe: 2020-10-01 (Week 39)'))
    ),
    
    mainPanel(
      tmapOutput("cvmtl_map")
    )
  )
  
)

server<-function(input, output) {
  
  output$cvmtl_map <- renderTmap({
    
    cvmtl <- filter(cvmtl, week == input$week)
    var_name <- ifelse(input$checkbox, 'adj', 'case')

      tm_shape(cvmtl) +  tm_polygons(var_name, title = "Cases", 
                                     style = input$select,
                                     breaks = c(0, 10, 20, 50, 100, 200, 400, Inf),
                                     legend.reverse = TRUE,
      ) +
        tm_layout(aes.palette = list(seq = "YlGnBu"))
  }) 
}

shinyApp(ui = ui, server = server)

I use a free ShinyApps account, active hour may be limited to 25 hours per month. There may be updated in the future in terms of functionality and data but the code is just a working example at the time of published date.

Linear Mixed Models via rstan

August 30, 2020 - 6 minutes
Statistics R

Saving Simple Features (sf) as RDS

April 10, 2020 - 1 minutes
Miscellaneous R

COVID-19 in Quebec

March 31, 2020 - 7 minutes
Statistics R ggplot2