2020-06-16-covid-19.en.Rmd

---
title: COVID-19 Trends and Insights
author: D. Kyle Ward
date: 'June 27, 2020'
slug: covid-19-trends
categories:
  - Health
tags: []
subtitle: ''
summary: ''
authors: []
lastmod: '2020-06-16T22:05:27-04:00'
featured: no
image:
  caption: ''
  focal_point: ''
  preview_only: no
projects: []
draft: false
editor_options: 
  chunk_output_type: console
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  echo = FALSE, fig.height = 4, warning = FALSE, out.width = "100%"
)
library(tidyverse)
library(lubridate)
library(zoo)
library(plotly)
library(here)
library(tidycensus)
library(rvest)
library(knitr)
library(kableExtra)
library(readxl)
library(htmltools)
```

```{r}
# Move this css tag outside the chunk to control the width of text
# on the page.

# <style type="text/css">
# .main-container {
#   max-width: 1000px;
#   margin-left: auto;
#   margin-right: auto;
# }
# </style>
```

```{r, message=FALSE}
# Get data from covidtracking.com and store a copy in this repo. To update the
# data, delete the nation and state CSV files in this repo.

nation_csv <- here("static", "data", "2020-06-16-covid-19", "nation.csv")
if (!file.exists(nation_csv)) {
  nation <- read_csv(url("https://covidtracking.com/api/v1/us/daily.csv"))
  nation <- nation %>%
    filter(states == 56) %>%
    mutate(date = ymd(date)) %>%
    arrange(date)
  write_csv(nation, nation_csv)
}

states_csv <- here("static", "data", "2020-06-16-covid-19", "states.csv")
if (!file.exists(states_csv)) {
  states <- read_csv(url("https://covidtracking.com/api/v1/states/daily.csv"))
  states <- states %>%
    filter(!is.na(dateChecked)) %>%
    mutate(
      date = ymd(date),
      region = state.region[match(state, state.abb)]
    ) %>%
    arrange(state, date)
  
  # Grab census population data by state and add to state table
  api_key <- Sys.getenv("CENSUS_API_KEY")
  acs <- get_acs(
    "state",
    variables = "B01003_001E",
    year = 2018
  )
  states <- states %>%
    left_join(
      acs %>% select(GEOID, population = estimate),
      by = c("fips" = "GEOID")
    ) %>%
    filter(!is.na(population))
  
  # Add governor data to the state table
  state_abb <- setNames(state.abb, state.name)
  wikiurl <- "https://en.wikipedia.org/wiki/List_of_current_United_States_governors_by_age"
  data <- wikiurl %>%
    xml2::read_html()
  gov_tbl <- html_table(data)[[1]] %>%
    mutate(state = as.character(state_abb[State])) %>%
    select(state, state_full = State, Party)
  states <- states %>%
    left_join(gov_tbl, by = "state") %>%
    filter(!is.na(Party))
  
  
  # Add lock down dates
  balloturl <- "https://ballotpedia.org/Status_of_lockdown_and_stay-at-home_orders_in_response_to_the_coronavirus_(COVID-19)_pandemic,_2020"
  data <- balloturl %>%
    xml2::read_html()
  data <- html_table(data, fill = TRUE)[[3]]
  data <- data[2:nrow(data), 1:2]
  colnames(data) <- c("state", "date")
  data <- data %>%
    separate(date, c("start", "stop"), sep = "-") %>%
    separate(stop, c(NA, "month", "day", NA), sep = " ") %>%
    separate(day, c("day", NA), sep = "\\[") %>%
    unite(stop, month, day, sep = " ") %>%
    mutate(
      start = paste(start, "2020"),
      start = mdy(start),
      stop = paste(stop, "2020"),
      stop = ifelse(grepl("TBD", stop), "June 16 2020", stop),
      stop = mdy(stop),
      days = stop - start,
      days = ifelse(is.na(days), 0, days)
    ) %>%
    select(state, lockdown_days = days)
  states <- states %>%
    left_join(data, by = c("state_full" = "state"))
  
  write_csv(states, states_csv)
}

nation <- read_csv(nation_csv)
states <- read_csv(states_csv)
```

```{r colors}
# colors used for various charts
default_blue <- I(rgb(31, 119, 180, maxColorValue = 255))
neg_test_tan <- I(rgb(216, 179, 101, maxColorValue = 255))
hospital_red <- I(rgb(222, 45, 38, maxColorValue = 255))
death_gray <- I(rgb(99, 99, 99, maxColorValue = 255))
```


The motivations for this post are simple: 

  1. There is a wealth of data on the COVID-19 pandemic.
  2. That data is misrepresented on a daily basis.
  
The most egregious sin is presenting counts of new cases without accounting for
increased testing. This distorts comparisons between states (and countries) but
also misrepresents temporal trends within them. I present the numbers both ways
in the analysis below and also offer an interactive application for viewing the
data. In addition, I compare the experience of states in general and California
to New York explicitly to offer some data-driven insight as to why states had
such drastically different experiences.

My primary source for COVID-19 data is 
[https://covidtracking.com/](https://covidtracking.com/). All data/code is
linked at the bottom of the article.

## National Data

### New Cases

The chart below shows the number of daily new cases and the 7-day trend. This
chart suggests that we saw a modest reduction in new cases and are now poised to
see an explosion. It's alarming and has prompted renewed wall-to-wall news
coverage of a "second wave". It looks (incorrectly) like an approaching tsunami.

```{r}
temp <- nation %>%
  mutate(sma = rollmean(positiveIncrease, 7, na.pad = TRUE, align = "right"))

plot_ly(data = temp) %>%
  add_trace(
    x = ~date, 
    y = ~round(positiveIncrease, -2),
    type = "bar", name = "New Cases",
    color = default_blue,
    alpha = .15, alpha_stroke = .15
  ) %>%
  add_trace(
    x = ~date, y = ~sma, type = "scatter", mode = "lines", hoverinfo = "none",
    name = "7-day Trend", color = default_blue
  ) %>%
  layout(
    title = "Daily New Cases",
    legend = list(x = .7, y = .95),
    yaxis = list(title = "Cases"),
    xaxis = list(title = "")
  ) %>%
  config(displayModeBar = FALSE)
```

The above graph has been presented in various forms by practically every news
outlet. It is also misleading. The chart below presents the same new-case data
in blue but also includes the number of negative tests shown by tan bars. The
amount of tests that come back negative is growing much faster than those that
come back positive.

```{r}
plot_ly(data = nation) %>%
  add_trace(
    x = ~date,
    y = ~round(positiveIncrease, -2),
    type = "scatter", mode = "line",
    name = "New Cases"
  ) %>%
  add_trace(
    x = ~date,
    y = ~round(negativeIncrease, -2),
    type = "bar",
    name = "Negative Tests",
    color = neg_test_tan
  ) %>%
  layout(
    title = "Daily New Cases and Negative Test Results",
    yaxis = list(title = "Tests"),
    xaxis = list(title = NA)
  ) %>%
  layout(
    legend = list(x = .1, y = .95),
    hovermode = "x unified"
  ) %>%
  config(displayModeBar = FALSE)
```

To accurately present the trend of new cases, they must be shown as a rate per
100 tests performed. With this metric, you can determine if the virus is
spreading or if we are simply doing more tests. The rate of new cases has begun
rising in the past few days, but these charts present a reality that is much
less alarming.

```{r}
temp <- nation %>%
  mutate(
    positive_rate = round(positiveIncrease / totalTestResultsIncrease * 100, 1),
    sma = rollmean(positive_rate, 7, na.pad = TRUE, align = "right")
  )
  
plot_ly(data = temp) %>%
  add_trace(
    x = ~date, y = ~positive_rate, type = "bar", name = "Case Rate",
    color = default_blue, alpha = .15, alpha_stroke = .15
  ) %>%
  add_trace(
    x = ~date, y = ~sma, type = "scatter", mode = "lines", hoverinfo = "none",
    name = "7-day Trend", color = default_blue
  ) %>%
  layout(
    title = "Daily New Cases per 100 Tests",
    yaxis = list(
      title = "New Cases per 100 Tests"
    ),
    xaxis = list(title = NA),
    legend = list(x = .7, y = .95)
  ) %>%
  config(displayModeBar = FALSE)
```

It is also worth noting that the high rates of positive cases in April were in
part the result of limited tests being reserved for people with COVID-like
symptoms. As testing increased, a more accurate measure of the infection rate
has emerged.

### Hospitalizations

Tracking daily hospital usage is a better way to gauge the impact
of COVID-19. This metric is not biased by issues with sample size or frame. In
the chart below, each bar represents the number of people in the hospital on
that day due to the virus. This metric fell from late April to mid June, but
has risen in the past four days.

```{r, warning=FALSE}
temp <- nation %>%
  mutate(
    sma = rollmean(hospitalizedCurrently, 3, na.pad = TRUE, align = "right")
  )

plot_ly(data = temp) %>%
  add_trace(
    x = ~date, y = ~hospitalizedCurrently,
    type = "bar",
    # type = "scatter", mode = "lines",
    name = "Hospitalized",
    color = hospital_red, alpha = .15, alpha_stroke = .15
  ) %>%
  add_trace(
    x = ~date, y = ~sma, type = "scatter", mode = "lines", hoverinfo = "none",
    name = "3-day Trend", color = hospital_red
  ) %>%
  layout(
    title = "Daily Hospitalized",
    yaxis = list(title = "Hospitalized"),
    xaxis = list(title = ""),
    legend = list(x = .7, y = .95)
  ) %>%
  config(displayModeBar = FALSE)
```

### Deaths

Similar to hospital usage, COVID-19 deaths offer a metric that is easier to
interpret. Deaths have fallen since late April; however, deaths are a lagging
indicator. In a second wave scenario, you would expect to see hospitalizations
go up first followed by rising deaths. The bar on June 25th is an odd data point
and likely a data reporting issue. (The [New York Times
data](https://github.com/nytimes/covid-19-data) lists a similar figure of
2,466.)

```{r}
temp <- nation %>%
  mutate(sma = rollmean(deathIncrease, 7, na.pad = TRUE, align = "right"))

plot_ly(data = temp) %>%
  add_trace(
    x = ~date, y = ~deathIncrease, type = "bar", name = "Deaths",
    color = death_gray,
    alpha = .15, alpha_stroke = .15
  ) %>%
  add_trace(
    x = ~date, y = ~sma, type = "scatter", mode = "lines", 
    name = "7-day Trend", color = death_gray
  ) %>%
  layout(
    title = "Daily Deaths",
    # showlegend = FALSE,
    xaxis = list(title = ""),
    yaxis = list(title = "Deaths"),
    legend = list(x = .7, y = .95)
  ) %>%
  config(displayModeBar = FALSE)
```

### National Summary

Perhaps the most important take away from the above review is that daily
hospitalization rates are likely the best statistic to broadcast to the public.
They don't require adjustment based on the amount of testing or lag behind as
much as death statistics. Absolute counts of new cases are the worst metric, but
this is where most news media are focused.

## State Data

While the nation as a whole is trending in the right direction, policy decisions
are made at the state level and individual states could be moving in opposite
directions. Different experiences in each state may also create natural
experiments for us to examine. In addition to the analysis I present here, I've
created an interactive viewer for the state data. Click the image below to open
the app and choose which metric to view, see trends for each state, and compare
states to each other.

<a href="http://us-covid-app.herokuapp.com/" target = "_blank">

![](/img/2020-06-16-covid-19/interactive_app.png)

</a>

### New Case Rate

The table below shows the rate of new cases per day for a handful of states.
Issues with the daily state data create some noise, but the different patterns
are easy to see.

```{r}
temp <- states %>%
  arrange(state, date) %>%
  filter(
    state %in% c("NC", "SC", "GA", "NY"),
    date > ymd("2020-03-20")
  ) %>%
  group_by(state) %>%
  mutate(
    positive_rate = round(positive / totalTestResults * 100, 1)
  ) %>%
  ungroup() %>%
  select(state, date, positive_rate)

plot_ly(data = temp) %>%
  add_trace(
    x = ~date, y = ~positive_rate, type = "scatter", mode = "markers",
    color = ~state
  ) %>%
  layout(
    title = "Daily New Cases per 100 Tests",
    yaxis = list(
      title = "New Cases per 100 Tests"
    ),
    xaxis = list(
      title = NA
    ),
    legend = list(x = .7, y = .95)
  ) %>%
  config(displayModeBar = FALSE)
```

### Death Rate

At the state level, comparing absolute deaths would be misleading. States like
New York and California will tend to have more deaths simply because they have
more people. As I did with new cases, I convert deaths to a rate: deaths per one
million people. Even with this adjustment, the chart below shows how severe the
problem became in NY compared to other states.

```{r}
temp <- states %>%
  arrange(state, date) %>%
  filter(
    state %in% c("NC", "SC", "GA", "NY", "CA", "TX"),
    date > ymd("2020-03-20")
  ) %>%
  mutate(
    death_rate = deathIncrease / population * 1000000
  )

plot_ly(data = temp) %>%
  add_trace(
    x = ~date, y = ~death_rate, type = "scatter", mode = "markers",
    color = ~state
  ) %>%
  layout(
    title = "Daily Deaths by State",
    yaxis = list(
      title = "Deaths per Million"
    ),
    xaxis = list(
      title = NA
    ),
    legend = list(x = .7, y = .95)
  ) %>%
  config(displayModeBar = FALSE)
```

## Regional Differences

The animation below shows the rates of hospitalization and death on a weekly
basis colored by region of the country. A state that experienced any week with a
death rate above 50 per million is labeled.

  * x-axis: deaths per million
  * y-axis: hospitalized per million
  * circle size: state population (Census ACS)
  * color: region

Press "Play" to see two primary patterns emerge:

  1. The separation of the northeastern states around New York City.
  2. Hospitalization and death rates have been falling for nearly all states in
  the past few weeks.

```{r}
# create the weekly stats data frame
weekly_stats_df <- states %>%
  arrange(state, date) %>%
  mutate(
    week = floor_date(date, "weeks"),
    week_id = week(date)
  ) %>%
  group_by(state, week) %>%
  summarize(
    days = n(),
    week_id = first(week_id),
    positiveIncrease = sum(positiveIncrease),
    totalTestResultsIncrease = sum(totalTestResultsIncrease),
    deathIncrease = sum(deathIncrease),
    hospitalizedCurrently = mean(hospitalizedCurrently, na.rm = TRUE),
    hospitalizedCurrently = ifelse(
      is.nan(hospitalizedCurrently), NA, hospitalizedCurrently),
    Party = first(Party),
    population = first(population),
    lockdown_days = first(lockdown_days),
    state_full = first(state_full),
    region = first(region)
  ) %>%
  # filter out early weeks when data doesn't exist for every state.
  # Also filter out any week with <7 days (the latest week)
  filter(week_id > 12, days == 7) %>%
  mutate(
    week_id = week_id - 13, # start sliders at 0
    positive_rate = round(positiveIncrease / totalTestResultsIncrease * 100, 1),
    # some errors in the reported data lead to 1-2 negative numbers
    positive_rate = pmax(positive_rate, 0),
    death_rate = round(deathIncrease / population * 1000000, 4),
    hosp_rate = round(hospitalizedCurrently / population * 1000000, 4),
    lockdown_class = case_when(
      lockdown_days > 60 ~ "1.Long",
      lockdown_days > 30 ~ "2.Medium",
      lockdown_days == 0 ~ "4.None",
      TRUE ~ "3.Short"
    ),
    lockdown_class = factor(
      lockdown_class,
      levels = c("1.Long", "2.Medium", "3.Short", "4.None"),
      ordered = TRUE
    )
  ) %>%
  ungroup()
```


```{r}
assert_zero <- weekly_stats_df %>%
  mutate(hosp_rate = ifelse(is.na(hosp_rate), 0, hosp_rate)) %>%
  # The last week right before publishing features an error/outlier in
  # the NJ data but is otherwise similar to the prior week.
  filter(week_id != 12)

# Create a DF of just those states to be labeled. Label any state that
# exceeded 50 dpm in any week.
labels <- assert_zero %>%
  group_by(state) %>%
  mutate(exceeded_50dpm = ifelse(max(death_rate) > 50, 1, 0)) %>%
  filter(exceeded_50dpm == 1)

num_of_states_over_50dpm <- length(unique(labels$state ))

plot_ly(data = assert_zero) %>%
  add_trace(
    x = ~death_rate, y = ~hosp_rate, 
    type = "scatter", mode = "markers",
    size = ~population,
    color = ~region,
    alpha = .05, alpha_stroke = .05, showlegend = FALSE,
    hoverinfo = "none"
  ) %>%
  add_trace(
    x = ~death_rate, y = ~hosp_rate, 
    type = "scatter", mode = "markers",
    marker = list(
      line = list(
        color = rgb(49, 52, 56, maxColorValue = 255)
      )
    ),
    color = ~region,
    ids = ~state, frame = ~week_id,
    size = ~population,
    hoverinfo = "text",
    text = ~paste0(
      "<br>", state, 
      "<br>", week,
      "<br>DPM: ", round(death_rate, 0),
      "<br>HPM: ", round(hosp_rate, 0)
    )
  ) %>%
  add_text(
    data = labels, text = ~state,
    x = ~death_rate, y = ~hosp_rate,
    ids = ~state, frame = ~week_id,
    textposition = "top", showlegend = FALSE
  ) %>%  
  layout(
    title = "Weekly Metrics by State",
    yaxis = list(
      title = "Hospitalized per Million"
    ),
    xaxis = list(
      title = "Deaths per Million",
      range = c(0, 273)
    )
  ) %>%
  config(displayModeBar = FALSE) %>%
  animation_opts(
    frame = 1500
  ) %>%
  animation_slider(
    currentvalue = list(
      prefix = "Weeks since 3/29: ",
      font = list(size = 12)
    )
  )
```

*Note: The states not reporting hospitalization data run along the x-axis
only*

## New York City

New York City is the [densest city in the
US](https://en.wikipedia.org/wiki/List_of_United_States_cities_by_population),
has the most [air
travel](https://en.wikipedia.org/wiki/List_of_busiest_city_airport_systems_by_passenger_traffic),
and has [transit ridership](/2019/09/03/understanding-transit-ridership-trends/)
that dwarfs other cities. These factors make it a hotbed for viral transmission.
The map below shows the peak death rate experienced in each state, and the
impact of New York City can be seen in the surrounding states. Outside of the
northeast, the peak death rates are much lower. Louisiana stands out in the
south, which the CDC
[attributes](https://www.cdc.gov/mmwr/volumes/69/wr/mm6915e4.htm?s_cid=mm6915e4_w)
to Mardi Gras.

```{r}
peak_deaths <- weekly_stats_df %>%
  group_by(state) %>%
  mutate(peak_death_rate = max(death_rate)) %>%
  filter(death_rate == peak_death_rate) %>%
  slice(1)

plot_geo() %>%
  add_trace(
    z = ~peak_deaths$peak_death_rate,
    hoverinfo = "text",
    # text = state.name,
    text = paste0(
      "<br>", state.name,
      "<br>Peak Rate: ", round(peak_deaths$peak_death_rate, 1)
    ),
    span = I(0),
    colors = "Greys",
    locations = peak_deaths$state, locationmode = 'USA-states'
  ) %>%
  layout(
    geo = list(
      scope = 'usa',
      projection = list(type = 'albers usa'),
      lakecolor = toRGB('white')    
    ),
    dragmode = FALSE
  ) %>%
  colorbar(
    title = paste0(
      "Peak Death Rate",
      "<br>(deaths per million)"
    )
  ) %>%
  config(displayModeBar = FALSE)
```

## Comparison to California

California is the most populous state, but compared to New York, it's experience
with the virus has been far less severe. Comparing the two states reveals
important similarities and differences that could explain the divergence. One
important similarity is that they initiated stay-at-home orders [within a day of
each
other](https://ballotpedia.org/Status_of_lockdown_and_stay-at-home_orders_in_response_to_the_coronavirus_(COVID-19)_pandemic,_2020).

### Density

New York City is the largest, densest city in the country. The table below shows
the top 10 US cities ranked by population density
([source](https://en.wikipedia.org/wiki/List_of_United_States_cities_by_population)).
San Francisco is second on the list, but has only 1/10th the population. In
fact, all the Californian cities in the table combined only hold 1.4 million
people compared to NYCs 8.3. This disparity in population and density is a
significant reason why NYC had a more severe outbreak.

```{r}
city_dens_url <- "https://en.wikipedia.org/wiki/List_of_United_States_cities_by_population"
data <- city_dens_url %>%
    xml2::read_html()
density_tbl <- html_table(data, fill = TRUE)[[5]]
density_tbl <- density_tbl[, c(2:4, 9)] %>%
  select(
    City, State = `State[c]`, Density = `2016 population density`,
    Population = `2019estimate`
  ) %>%
  mutate(
    City = gsub("\\[.*", "", City),
    density_num = gsub("/sq\\smi|,", "", Density),
    density_num = as.numeric(density_num)
  ) %>%
  arrange(desc(density_num))

density_tbl %>%
  mutate(
    State = cell_spec(State, "html", color = case_when(
      State == "California" ~ "red",
      State == "New York" ~ "blue",
      TRUE ~ "gray"
    )),
    City = cell_spec(City, "html", color = case_when(
      .$State == "California" ~ "red",
      .$State == "New York" ~ "blue",
      TRUE ~ "gray"
    ))
  ) %>%
  # filter(density_num >= 10000) %>%
  head(10) %>%
  select(-density_num) %>%
  kable(format = "html", escape = FALSE, row.names = FALSE) %>%
  kable_styling(full_width = FALSE)
```

### Air Travel

Air travel is another important consideration for viral spread. The table below
shows the top six US metro regions in terms of air travel
([source](https://en.wikipedia.org/wiki/List_of_busiest_city_airport_systems_by_passenger_traffic)).
While NYC tops the list, Los Angeles and San Francisco are not far behind. One
potentially important difference not captured below is that NYC attracts more
European tourists, which is [where their virus
originated](https://www.nytimes.com/2020/04/08/science/new-york-coronavirus-cases-europe-genomes.html).
This may or may not prove significant after further study.

```{r}
air_url <- "https://en.wikipedia.org/wiki/List_of_busiest_city_airport_systems_by_passenger_traffic"
data <- air_url %>%
    xml2::read_html()
air_tbl <- html_table(data, fill = TRUE)[[5]] %>%
  filter(`Metropolitan area` %in% c(
    "New York City", "Los Angeles", "Atlanta", "Chicago",
    "Miami", "San Francisco Bay Area"
  )) %>%
  select(
    Metro = `Metropolitan area`,
    `Yearly Passengers` = Totalpassengers,
    `Airport(s)` = `Airport(s) included`
  )

air_tbl %>%
  mutate(
    Metro = cell_spec(Metro, "html", color = case_when(
      Metro %in% c("Los Angeles", "San Francisco Bay Area") ~ "red",
      Metro == "New York City" ~ "blue",
      TRUE ~ "gray"
    ))
  ) %>%
  kable(format = "html", escape = FALSE, row.names = FALSE) %>%
  kable_styling(full_width = FALSE)
```

### Transit

Crowded buses and subways are an excellent place for viruses to spread. In my
post on [transit ridership
trends](/2019/09/03/understanding-transit-ridership-trends/), I used the chart
below to show just how much larger the transit market is in NYC compared to the
rest of the country.

```{r}
xlsx <- here::here(
  "static", "data", "2020-06-16-covid-19", "April 2020 Adjusted Database.xlsx"
)
upt_tbl_raw <- read_excel(xlsx, sheet = "UPT")

upt_tbl <- upt_tbl_raw %>%
  filter(Active == "Active") %>%
  gather(key = "month_year", value = "trips", JAN02:APR20) %>%
  mutate(
    year = substr(
      month_year, start = str_length(month_year) - 1, stop = str_length(month_year)
    ),
    month = substr(month_year, start = 1, stop = 3),
    date = parse_date_time2(month_year, orders = "my")
  ) %>%
  arrange(Agency, date)

# shortens the agency name to the first 4 characters. Improves legend.
shorten_agency <- function(Agency){
  Agency = gsub("[^a-zA-Z0-9 _]", "", Agency)
  word_count = str_count(Agency, "\\W+") + 1
  agency_short = word(Agency, 1, pmin(word_count, 4))
  return(agency_short)
}
```

```{r, out.width="100%", warning=FALSE}
data1 <- upt_tbl %>%
  filter(month_year == "JUN19") %>%
  group_by(Agency) %>%
  summarize(Trips = round(sum(trips, na.rm = TRUE), -4)) %>%
  arrange(desc(Trips)) %>%
  head(10) %>%
  arrange(Trips) %>%
  mutate(Agency = shorten_agency(Agency))

plot_ly(
  data1, y = ~Agency, x = ~Trips, type = "bar", orientation = "h",
  text = ~Agency, textposition = "auto", hoverinfo = 'x'
) %>%
  layout(
    title = "Top 10 Transit Agencies (June 2019 trips)",
    yaxis = list(
      categoryorder = "array",
      categoryarray = ~Agency,
      showticklabels = FALSE
    )
  ) %>%
  config(displayModeBar = FALSE)

transit_multiple <- round(
  data1$Trips[data1$Agency == "MTA New York City"] / 
    data1$Trips[data1$Agency == "Los Angeles County Metropolitan"],
  1
)
```

Unlike the modest difference in air travel, NYC transit is `r transit_multiple`
times larger than the largest market in California during normal operation.
People crammed into enclosed spaces on the NYC city subway were primed to spread
the disease. As the virus spread, San Francisco [closed their
subway](https://sf.curbed.com/2020/3/26/21195566/muni-trains-covid-19-sf-underground-subway)
while Los Angeles cut back transit service after [steep drops in
ridership]((https://www.latimes.com/california/story/2020-04-17/coronavirus-cuts-los-angeles-metro-bus-train-service)).
NYC also saw reduced ridership and ran reduced service, but they only closed the
subway between [1:00 am and 5:00
am](https://abc7ny.com/subway-cleaning-plan-mta-coronavirus-nyc/6155622/).

The chart below updates my original 2019 chart with April 2020 data. During the
peak of the viral outbreak, NYC still recorded more than twice the ridership of
Los Angeles.

```{r, out.width="100%", warning=FALSE}
april20_data <- upt_tbl %>%
  filter(month_year == "APR20") %>%
  group_by(Agency) %>%
  summarize(Trips = round(sum(trips, na.rm = TRUE),-4)) %>%
  arrange(desc(Trips)) %>%
  head(10) %>%
  arrange(Trips) %>%
  mutate(Agency = shorten_agency(Agency))

plot_ly(
  april20_data, y = ~Agency, x = ~Trips, type = "bar", orientation = "h",
  text = ~Agency, textposition = "auto", hoverinfo = 'x'
) %>%
  layout(
    title = "Top 10 Transit Agencies (April 2020 trips)",
    yaxis = list(
      categoryorder = "array",
      categoryarray = ~Agency,
      showticklabels = FALSE
    )
  ) %>%
  config(displayModeBar = FALSE)
```

### Weather

Weather also differentiates California from New York. The chart below compares
average temperatures in San Francisco and Los Angeles with New York City
([source](https://www.ncdc.noaa.gov/cdo-web)). California remained noticeably
warmer than NYC through early April. There is still debate over the influence of
weather on COVID-19, but [some
researchers](https://www.npr.org/sections/goatsandsoda/2020/04/09/830297538/scientists-try-to-figure-out-if-summer-will-slow-the-spread-of-covid-19)
[think](https://www.cebm.net/covid-19/do-weather-conditions-influence-the-transmission-of-the-coronavirus-sars-cov-2/)
it could behave like it's well-known coronavirus cousins.

```{r}
weather_csv <- here(
  "static", "data", "2020-06-16-covid-19", "noaa_weather_2020.csv"
)
weather_tbl <- read_csv(
  weather_csv,
  col_types = cols(
    STATION = col_character(),
    NAME = col_character(),
    LATITUDE = col_double(),
    LONGITUDE = col_double(),
    ELEVATION = col_double(),
    DATE = col_character(),
    TAVG = col_double(),
    TMAX = col_double(),
    TMIN = col_double()
  )
) %>%
  separate(NAME, sep = ", ", into = c("NAME", "REGION")) %>%
  filter(!is.na(TAVG), REGION != "NJ US") %>%
  group_by(REGION, DATE) %>%
  summarize(
    LATITUDE = first(LATITUDE),
    LONGITUDE = first(LONGITUDE),
    TAVG = round(mean(TAVG), 1),
    TMAX = round(mean(TMAX), 1),
    TMIN = round(mean(TMIN), 1)
  ) %>%
  mutate(
    sma = rollmean(TAVG, 10, aligh = "right", na.pad = TRUE)
  )

plot_ly() %>%
  add_trace(
    data = weather_tbl, x = ~DATE, y = ~TAVG,
    color = ~REGION,
    type = "scatter", mode = "lines",
    alpha = .15, alpha_stroke = .15, showlegend = FALSE,
    hoverinfo = "none"
  ) %>%
  add_trace(
    data = weather_tbl, x = ~DATE, y = ~sma,
    color = ~REGION,
    type = "scatter", mode = "lines"
  ) %>%
  layout(
    title = "Average Temperature",
    yaxis = list(title = "Average Temperature"),
    xaxis = list(title = "Date"),
    hovermode = "x unified",
    legend = list(x = .1, y = .95)
  ) %>%
  config(displayModeBar = FALSE)

```

## What about Lockdowns?

Lockdowns are the first factor many consider when COVID-19 outbreaks between
states. Intuition might suggest that states with weak or no lockdown
restrictions would have higher rates cases, hospitalization, and death. Instead,
the comparison of California to New York shows that other factors are more
likely to impact severity.

The animation from before is now color coded by lockdown duration
([source](https://ballotpedia.org/Status_of_lockdown_and_stay-at-home_orders_in_response_to_the_coronavirus_(COVID-19)_pandemic,_2020)).
What it shows is that northeastern states implemented longer lockdowns in
response to a much more severe viral outbreak. Outside of the northeast, states
adopted different approaches and no clear pattern emerged showing that longer or
shorter durations made any meaningful difference.
  
```{r}
plot_ly(data = assert_zero) %>%
  add_trace(
    x = ~death_rate, y = ~hosp_rate, 
    type = "scatter", mode = "markers",
    color = ~lockdown_class,
    colors = c("blue", "gray", "orange", "red"),    
    size = ~population,
    alpha = .05, alpha_stroke = .05, showlegend = FALSE,
    hoverinfo = "none"
  ) %>%
  add_trace(
    x = ~death_rate, y = ~hosp_rate, 
    type = "scatter", mode = "markers",
    marker = list(
      line = list(
        color = rgb(49, 52, 56, maxColorValue = 255)
      )
    ),
    color = ~lockdown_class, 
    colors = c("blue", "gray", "orange", "red"),
    ids = ~state, frame = ~week_id,
    size = ~population,
    hoverinfo = "text",
    text = ~paste0(
      "<br>", state,
      "<br>", week,
      "<br>DPM: ", round(death_rate, 0),
      "<br>HPM: ", round(hosp_rate, 0),
      "<br>Lockdown: ", lockdown_days, " days"
    )
  ) %>%
  add_text(
    data = labels, text = ~state,
    x = ~death_rate, y = ~hosp_rate, 
    ids = ~state, frame = ~week_id,
    textposition = "top", showlegend = FALSE
  ) %>%  
  layout(
    title = "Weekly Metrics by State",
    yaxis = list(title = "Hospitalized per Million"),
    xaxis = list(
      title = "Deaths per Million",
      range = c(0, 273)
    ),
    legend = list(
      title = list(text = "Lockdown<br>Duration")
    )
  ) %>%
  config(displayModeBar = FALSE) %>%
  animation_opts(frame = 1500) %>%
  animation_slider(
    currentvalue = list(
      prefix = "Weeks since 3/29: ",
      font = list(size = 12)
    )
  )
```

# Conclusion

One of the unique things about the modern day is the abundance of information.
This is a great benefit, but it is important that we use and present that
information responsibly. Raw counts must be handled appropriately before making
comparisons across time or geography. Failure to do so will lead to poor
decision making.

The severity of the virus has varied across the country and is likely driven by
underlying differences like density, transit utilization, and weather. Rural
states with shorter or no lockdowns had milder outbreaks than the northeast
in spite of stronger restrictions. Finally, a second wave does appear to be
forming as states re-open. This fact alone does not make it clear what the right
policy response should be. Lockdowns have [negative
consequences](https://www.npr.org/sections/goatsandsoda/2020/06/16/874198026/the-cost-of-thailands-coronavirus-success-despair-and-suicide)
for our lives and livelihoods. Trade offs between these impacts and viral spread
are difficult to balance.

## Reproducibility and Data

All the data and code used to perform this analysis is available on GitHub:

[Data](https://github.com/dkyleward/blog/tree/master/static/data/2020-06-16-covid-19)  
[Code](https://github.com/dkyleward/blog/blob/master/content/post/2020-06-16-covid-19/2020-06-16-covid-19.en.Rmd) 

Other sources of interest:  
[https://covidtracking.com/](https://covidtracking.com/)  
[NY Times Data](https://github.com/nytimes/covid-19-data)