library(tidyverse)
library(sf)
library(leaflet)
library(leafpop)
library(formattable)
library(ggiraph)
# install.packages("devtools")
# devtools::install_github("UrbanInstitute/urbnmapr")
library(urbnmapr) # change to include AK/HI on map
library(cowplot)
library(RcppRoll)
library(htmlTable)
library(ggtext)
library(patchwork)
library(Hmisc)
library(httr)
library(rvest)
# require(RSocrata); require(plyr); require(readr); require(here)
# library(rvest)
# library(albersusa)
# Default ggplot theme
theme_set(theme_bw() +
theme(axis.title.x = element_text(size = rel(1.25)),
axis.title.y = element_text(size = rel(1.25)),
plot.caption = element_text(face = "italic", size = 9),
legend.position = "bottom",
legend.background = element_rect(fill = "transparent"),
panel.background = element_rect(fill = "transparent"),
panel.border = element_blank(),
axis.line = element_line(color="black"),
panel.grid.major.x = element_blank(),
panel.grid.minor.x = element_blank(),
panel.grid.minor.y = element_blank(),
panel.grid.major.y = element_line(color = "grey90", linetype = 2)
))
Last Data Update: 2023-12-19 21:31:13.102324. All source documents for this report can be found on github
This document updates weekly each Wednesday after the TN State DOH releases data (usually Wednesday at ~4PM). The CDC releases data each weekday on their site.
For official, weekly TN State Dept of Health Updates, go to the Critical Indicator Report.
For a reliable discussion of COVID-19, vaccines, and emerging variants, listen to Dr. Aronoff’s Medicine Grand Rounds from August 2021.
All ages over 6 months are now eligible for COVID vaccination- Go to one of these sites to sign up:
Update 5/10/2021 New Cases now reported only by the date of detection. Will update with details once it becomes clear how these results will be mixed with the prior reporting method (date of report).
The Tennessee state information is simply a visual summary of the data directly downloaded from TN state Dept of Health, by date reported (rather than date of test). I do not attempt to model, make predictions, or give advice on how to use this data. If you want that sort of thing, go to the VUMC Health Policy site, which has developed some very good local models and Advisory Memos for Tennessee.
For the most up-to-date numbers visit the TN state Dept of Health or local health department sites. In particular, the Metro. Nashville Dashboard has a very nice summary.
Other useful sites for following trends in Tennessee and nationally:
These data are taken directly from the TN state Dept of Health. Previously the DOH reported results by date of report, rather than date of test. They transitioned to reporting results by Date of Test on 5/10/2021, although it is not clear if data reported prior to that will be cleaned up.
Hospitalizations: A current summary of Hospitalizations can be found on the TN DOH site.
# READ IN AND format data ####
tndoh_daily <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-Daily-Case-Info.xlsx"))
names(tndoh_daily) <- tolower(names(tndoh_daily))
# remove the new_test results which our bulk reported or adjusted
tndoh_daily <-
tndoh_daily %>%
mutate(date = lubridate::ymd(date)) %>%
arrange(date) %>%
filter(date > lubridate::ymd("2000-01-01")) # State occasionally releases data with default SAS date of "1960-01-01"for missing values
# min(tndoh_daily$date)
# tndoh_daily <-
# tndoh_daily %>%
# filter(date > lubridate::ymd("2020-01-01")) # data error introduced; filter out erroneous entry for 1960-01-01
tndoh_daily[tndoh_daily$date %in% lubridate::ymd(c("2020-07-13", "2020-07-11", "2020-06-12")), "new_tests"] <- NA_real_
# STATE POPULATION DATA ####
# us_data <- readRDS(here::here("data_rds/us_data.rds"))
us_data <- get_urbn_map("states", sf = TRUE)
us_data <-
us_data %>%
mutate(name = state_name,
stusps = state_abbv)
# state_pop <- readxl::read_excel(here::here("state_populations.xlsx"))
# state_pop$region <- gsub("\\.", "", state_pop$region)
# state_pop <-
# state_pop %>%
# pivot_longer(cols = -c(region, census_2010),
# names_to = "pop_year",
# names_prefix = "pop_",
# values_to = "pop_estimate")
# saveRDS(state_pop, here::here("data_rds/state_pop.rds"))
state_pop <- readRDS( here::here("data_rds/state_pop.rds"))
state_pop_2019 <-
state_pop %>%
filter(pop_year == 2019) %>%
filter(region %in% state.name)
TN_pop2019 <- state_pop_2019$pop_estimate[state_pop_2019$region == "Tennessee"]
# vaccine state info
tndoh_vax_state <- readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_STATE_SUMMARY.XLSX"))
names(tndoh_vax_state) <- tolower(names(tndoh_vax_state))
tndoh_vax_state <-
tndoh_vax_state %>%
mutate(date = lubridate::ymd(date),
vaccine_count_pct = vaccine_count/TN_pop2019*100,
recipient_count_pct = recipient_count/TN_pop2019*100,
recip_fully_vacc_pct = recip_fully_vacc/TN_pop2019*100)
# add trend line calculated by STL (Seasonal and Trend decomposition using Loess)
tndoh_daily$new_cases_ts <- ts(tndoh_daily$new_cases, frequency = 7, start = 1)
tndoh_daily$new_cases_trend.stl <- stl(tndoh_daily$new_cases_ts, s.window = 7)$time.series[, "trend"]
# add Deaths trend line calculated by STL (Seasonal and Trend decomposition using Loess)
tndoh_daily$new_deaths_ts <- ts(tndoh_daily$new_deaths_by_dod, frequency = 7, start = 1)
tndoh_daily$new_deaths_trend.stl <- stl(tndoh_daily$new_deaths_ts, s.window = 7)$time.series[, "trend"]
# create annotations
tndoh_daily <-
left_join(tndoh_daily, tndoh_vax_state,
by = "date") %>%
mutate(pos_tests_new = pos_tests -lag(pos_tests, 1),
pos_pct = pos_tests/total_tests,
pos_pct_today = pos_tests_new/new_tests,
hospitalized_tot = replace_na(total_hosp, 0),
hospitalized_new = total_hosp - lag(total_hosp),
new_cases_100k = new_cases/TN_pop2019*100000,
new_cases_7davg = RcppRoll::roll_meanr(new_cases, 7),
new_cases_7davg_100k = new_cases_7davg/TN_pop2019*100000,
new_deaths_7davg = RcppRoll::roll_meanr(new_deaths_by_dod, 7),
new_hosp_7davg = RcppRoll::roll_meanr(new_hosp, 7),
new_tests_7davg = RcppRoll::roll_meanr(new_tests, 7),
pos_pct_7davg = RcppRoll::roll_meanr(pos_tests_new, 7)/RcppRoll::roll_meanr(new_tests, 7),
tt_daily = sprintf("<tr><td><b>%s</b></td></tr><br>
<tr><td><b>New Cases: %s</b></td></tr><br>
<tr><td><b>New Cases/100k: %s</b></td></tr><br>
<tr><td><b>New Hospitalizations: %s</b></td></tr><br>
<tr><td><b>New Deaths: %s</b></td></tr><br>
<tr><td><b>Test+ %% Today: %1.1f%%</b></td></tr><br>
<tr><td>New Cases 7d Avg/100k: %s</td></tr><br>
<tr><td>New Cases trend: %s</td></tr><br>
<tr><td>Deaths 7d Avg: %1.1f</td></tr><br>
<tr><td>Test+ %% 7d Avg: %1.1f%%</td></tr><br>
<tr><td>Case Fatality Ratio: %1.2f%%</td></tr><br>
<tr><td>New Vaccinations: %s</td></tr><br>
<tr><td>Rec 1st Vaccination %%Pop: %1.2f%%</td></tr>",
paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date), ", ", lubridate::year(date), " (", lubridate::wday(date, label = T), ")"),
comma(new_cases, digits = 0),
comma(new_cases_100k, digits = 0),
comma(new_hosp, 0),
comma(new_deaths_by_dod, digits = 0),
pos_pct_today*100,
comma(new_cases_7davg_100k, 1),
comma(new_cases_trend.stl/TN_pop2019*100000, 1),
comma(new_deaths_7davg, digits = 0),
pos_pct_7davg*100,
total_deaths_by_dod/total_cases*100,
comma(new_vaccine_count, 0),
recipient_count_pct) )
tndoh_daily_limited <-
tndoh_daily %>%
select(date, total_cases, total_hosp, total_deaths_by_dod) %>%
pivot_longer(cols = c(total_cases, total_hosp, total_deaths_by_dod),
names_to = "measure",
values_to = "count") %>%
mutate(var.f = factor(measure, levels = c("total_cases", "total_hosp", "total_deaths_by_dod"),
labels = c("Cases", "Hospitalizations", "Deaths")))
title_tnsdc <- cowplot::ggdraw() +
cowplot::draw_label(paste0("COVID-19 in Tennessee (as of ", max(tndoh_daily$date), ")") , fontface='bold')
tndoh_daily_long <-
tndoh_daily %>%
select(date, new_cases, new_deaths_by_dod, new_hosp) %>%
pivot_longer(cols = c(new_cases, new_deaths_by_dod, new_hosp),
names_to = "var",
values_to = "value") %>%
mutate(var.f = factor(var, levels = c("new_cases", "new_deaths_by_dod", "new_hosp"),
labels = c("Cases", "Hospitalizations", "Deaths")))
label_pos_pct <-
with(tndoh_daily,
sprintf("Overall Positive: %1.1f%%\nToday's Positive: %1.1f%%",
pos_pct[which(date == max(date))]*100,
pos_pct_today[which(date == max(date))]*100))
# AGE DATA ####
tndoh_age <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-Age.xlsx"))
names(tndoh_age) <- tolower(names(tndoh_age))
tndoh_age <-
tndoh_age %>%
filter(age_range != "Pending") %>%
mutate(date = lubridate::ymd(date),
age_range = gsub(" years", "", age_range),
age.f = factor(age_range),
deaths_case_rate = case_when(age.f == "Pending" ~ NA_real_,
TRUE ~ ar_totaldeaths/ar_totalcases*100))
# TN age groups 2019 estimates(taken from https://www.census.gov/data/tables/time-series/demo/popest/2010s-state-detail.html)
tn_age_groups <-
tibble::tibble(age.f = factor(c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"),
levels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")),
populationtn_byage_2019 = c(908244, 858307, 943700, 879850, 854356, 893802, 782810, 481538, 226567))
# c(821984, 855582, 942211, 887898, 850125, 893578, 808164, 516865, 252767)) # OLD ESTIMATES; RECALCULATED
tn_age_groups2 <-
tibble::tibble(age.f2 = factor(c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"), # 16-20 needs to match other age category
levels = c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")),
populationtn_byage_2019_B = c(408605, 585900, 343637, 428409, 943700, 879850, 854356, 893802, 782810, 481538, 226567))
tndoh_age <-
left_join(tndoh_age,
tn_age_groups,
by = "age.f") %>%
group_by(age.f) %>%
mutate(ar_new_cases_7davg = round(RcppRoll::roll_meanr(ar_newcases, 7), 1),
ar_new_cases_100k = ar_newcases/populationtn_byage_2019*100000,
ar_new_cases_100k_7davg = round(RcppRoll::roll_meanr(ar_newcases, 7)/populationtn_byage_2019*100000, 1),
ar_totaldeaths_100k = ar_totaldeaths/populationtn_byage_2019*100000,
ar_new_deaths_7davg = round(RcppRoll::roll_meanr(ar_newdeaths, 7), 1),
ar_new_deaths_100k = ar_newdeaths/populationtn_byage_2019*100000,
ar_new_deaths_100k_7davg = round(RcppRoll::roll_meanr(ar_newdeaths, 7)/populationtn_byage_2019*100000, 1)) %>%
ungroup()
tndoh_age <-
tndoh_age %>%
group_by(date) %>%
mutate(total_new_cases_100k_7davg = sum(ar_new_cases_100k_7davg, na.rm = T),
ar_proportion_100k_7davg = ar_new_cases_100k_7davg/total_new_cases_100k_7davg,
ar_newcases_sum_7davg = sum(ar_new_cases_7davg, na.rm = T),
ar_proportion_7davg = ar_new_cases_7davg/ar_newcases_sum_7davg,
ar_deaths_sum_7davg = sum(ar_new_deaths_7davg, na.rm = T),
ar_newdeaths_proportion_7davg = ar_new_deaths_7davg/ar_deaths_sum_7davg) %>%
ungroup()
# VACCINE by age
tndoh_vax_age <- readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_AGE_GROUPS.XLSX"),
col_types = c("text", "date", rep("numeric", 8)))
names(tndoh_vax_age) <- tolower(names(tndoh_vax_age))
tndoh_vax_age <-
tndoh_vax_age %>%
filter(!grepl("pending", age_group, ignore.case = T )) %>%
mutate(date = lubridate::ymd(date),
age.f2 = factor(age_group,
levels = c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"),
labels = c("<5", "5-11", "12-15", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")))
# tndoh_vax_age <-
# tndoh_vax_age %>%
# filter(!grepl("pending", age_group, ignore.case = T )) %>%
# mutate(date = lubridate::ymd(date),
# age.f = factor(age_group,
# levels = c("0-10", "11-20", "16-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"), # 16-20 needs to match other age category
# labels = c("0-10", "11-20", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")))
tndoh_vax_age <-
left_join(tndoh_vax_age, tn_age_groups2,
by = "age.f2") %>%
rowwise() %>%
mutate(ar_totalvax_pct = vaccine_count/populationtn_byage_2019_B*100,
ar_recipientvax_pct = recipient_count/populationtn_byage_2019_B*100) %>%
ungroup()
tndoh_age_last <- # FOR AGE PROPORTION PLOT
tndoh_age %>%
filter(!is.na(ar_proportion_7davg)) %>%
filter(date == max(date)) %>%
arrange(desc(age_range)) %>%
mutate(sec_axis.ticks_cases = cumsum(ar_proportion_7davg),
sec_axis.label_cases = sprintf("<b>%s:</b> %1.1f%%", age_range, ar_proportion_7davg*100),
sec_axis.ticks_deaths = cumsum(ar_newdeaths_proportion_7davg),
sec_axis.label_deaths = sprintf("<b>%s:</b> %1.1f%%", age_range, ar_newdeaths_proportion_7davg*100)) %>%
select(age_range, ar_proportion_7davg, ar_newdeaths_proportion_7davg,
sec_axis.ticks_cases, sec_axis.label_cases,
sec_axis.ticks_deaths, sec_axis.label_deaths) %>%
mutate(sec_axis.ticks_cases_mid = (sec_axis.ticks_cases + lag(sec_axis.ticks_cases, default = 0))/2,
sec_axis.ticks_deaths_mid = (sec_axis.ticks_deaths + lag(sec_axis.ticks_deaths, default = 0))/2)
# most recent vaccination by age
vax_by_age_mostrecent <-
tndoh_vax_age %>%
filter(!is.na(vaccine_count)) %>%
filter(date == max(date, na.rm = T)) %>%
select(date,age.f2, recipient_count, populationtn_byage_2019_B, ar_recipientvax_pct)
mostrecent_vax_agedate <- max(vax_by_age_mostrecent$date)
vax_by_age_mostrecent_refactor <- select(vax_by_age_mostrecent, -populationtn_byage_2019_B)
vax_by_age_mostrecent_refactor <- vax_by_age_mostrecent_refactor[!vax_by_age_mostrecent_refactor$age.f2 %in% c("5-11", "12-15", "16-20"), ]
vax_by_age_mostrecent_refactor <-
vax_by_age_mostrecent_refactor %>%
mutate(age.f = factor(age.f2,
levels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+"))) %>%
bind_rows(
tibble(age.f = factor(c("0-10", "11-20"),
levels = c("0-10", "11-20", "21-30", "31-40", "41-50", "51-60", "61-70", "71-80", "81+")),
recipient_count = c(vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "5-11"]*5/6, # 0 for <5 not here
vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "5-11"]*1/6 +
vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "12-15"] +
vax_by_age_mostrecent$recipient_count[vax_by_age_mostrecent$age.f2 == "16-20"]))) %>%
left_join(tn_age_groups, by = "age.f") %>%
mutate(ar_recipientvax_pct = recipient_count/populationtn_byage_2019*100)
# Deaths by age table ####
age_death_table <-
tndoh_age %>%
filter(date == max(date, na.rm = T)) %>%
select(date, age.f,populationtn_byage_2019, ar_totalcases, ar_totaldeaths, ar_totaldeaths_100k) %>%
mutate(#ar_totalcases = comma(ar_totalcases, 0),
ar_totaldeaths_pct = sprintf("%1.2f%%", ar_totaldeaths_100k/1000),
ar_totaldeaths_100k = sprintf("%1.2f", ar_totaldeaths_100k)) %>%
left_join({select(vax_by_age_mostrecent_refactor,
age.f, recipient_count, ar_recipientvax_pct )},
by = "age.f")
# COUNTY DATA ####
# tndoh_countynew <- readxl::read_excel(here::here("tn_data_center/TN_COVID19_CountyDaily.xlsx"))
TN_counties <- readRDS(here::here("data_rds/TN_counties.rds")) # see below for source
tndoh_county_coln <- length(names(readxl::read_excel(here::here("data_tndoh/Public-Dataset-County-New.xlsx"), n_max = 2)))
tndoh_countynew <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-County-New.xlsx"),
col_types = c("text","date", rep("numeric", tndoh_county_coln-2)))
names(tndoh_countynew) <- tolower(names(tndoh_countynew))
# create STL new case trend
tndoh_countynew <-
tndoh_countynew %>%
group_by(county) %>%
arrange(date) %>%
nest() %>%
mutate(new_cases_county_ts = map(data, ~ts(.x$new_cases, frequency = 7, start = 1)),
new_cases_county_trend = map(new_cases_county_ts, ~stl(.x, s.window = 7)$time.series[, "trend"])) %>%
select(county, data, new_cases_county_trend) %>%
unnest(cols = c(data, new_cases_county_trend))
tndoh_countynew <-
tndoh_countynew %>%
mutate(date = lubridate::ymd(date),
pct_pos = total_cases/(total_tests)*100,
county_lower = tolower(county)) %>%
group_by(county) %>%
arrange(date) %>%
mutate(new_cases_county_4davg = round(RcppRoll::roll_meanr(new_cases, 4), 1),
new_cases_county_7davg = round(RcppRoll::roll_meanr(new_cases, 7), 1),
max_new_cases_trend = max(new_cases_county_trend, na.rm = T),
new_cases_county_trend_ratio = new_cases_county_trend/max_new_cases_trend,
new_deaths_county_7davg = round(RcppRoll::roll_meanr(new_deaths, 7), 1),
new_testpos_pct = round(new_pos_tests/new_tests*100, 1),
new_testpos_pct_7davg = round(RcppRoll::roll_meanr(new_pos_tests, 7)/RcppRoll::roll_meanr(new_tests, 7)*100, 1),
new_hospitalized_county_7davg = round(RcppRoll::roll_meanr(new_hospitalized, 7), 1),
new_cases_county_7davg_1wkdelta = (new_cases_county_7davg - lag(new_cases_county_7davg, 7))/lag(new_cases_county_7davg, 7)
) %>%
ungroup()
# TN_data <- tigris::counties(state = "TN", cb = TRUE)
# names(TN_data) <- tolower(names(TN_data))
# TN_counties <- sf::st_as_sf(TN_data) %>%
# mutate(county_lower = tolower(as.character(name)))
# TN COUNTY POPULATION FROM US CENSUS 2019 ESTIMATES
# tn_county_pop <- read_csv(here::here("data_csv/TN 2019 county pop estimates.csv"))
# tn_county_pop <-
# tn_county_pop %>%
# mutate(county_lower = tolower(county)) %>%
# select(-county)
# TN_counties <-
# left_join(TN_counties,
# tn_county_pop,
# by = "county_lower")
# saveRDS(TN_counties, here::here("data_rds/TN_counties.rds"))
# TN_counties <- readRDS(here::here("data_rds/TN_counties.rds"))
tndoh_countynew <-
left_join(tndoh_countynew,
TN_counties,
by = "county_lower") %>%
mutate(total_cases_per100kpop = total_cases/population*100000,
new_cases_per100kpop_7davg = new_cases_county_7davg /population*100000,
new_cases_county_trend_per100kpop = new_cases_county_trend/population*100000,
total_tests_pctpop = total_tests/population*100,
new_tests_per100kpop = new_tests/population*100000,
total_deaths_per100kpop = total_deaths/population*100000)
tndoh_countynew_calc <-
tndoh_countynew %>%
group_by(date) %>%
filter(tolower(county) %in% TN_counties$county_lower) %>% # exclude out of state and pending
mutate(#new_testpos_pct = new_pos_tests/new_tests*100,
new_cases_county_7davg_rank = rank(-new_cases_county_7davg),
new_cases_per100kpop_7davg_rank = rank(-new_cases_per100kpop_7davg),
new_cases_rank = rank(-new_cases)) %>%
ungroup() %>%
select(county, date, new_cases_county_7davg_rank,
new_cases_per100kpop_7davg_rank, new_cases_rank)
tndoh_countynew <-
left_join(tndoh_countynew,
tndoh_countynew_calc,
by = c("date", "county"))
# DEMOGRAPHICS
tndoh_demographics <- readxl::read_excel(here::here("data_tndoh/Public-Dataset-RaceEthSex.xlsx"))
names(tndoh_demographics) <- tolower(names(tndoh_demographics))
tndoh_demographics$cat_detail <- gsub("/ ", "/", tndoh_demographics$cat_detail)
tndoh_demographics <-
tndoh_demographics %>%
filter(category == "RACE",
date == max(date)) %>%
dplyr::rename(race = cat_detail,
cases_byrace = cat_totalcases) %>%
mutate(date = lubridate::ymd(date),
race = factor(race,
levels = c("White", "Black or African American", "Asian", "Other/Multiracial", "Pending", "American Indian or Alaska Native", "Native Hawaiian or Other Pacific Islander"),
labels = c("White", "Black", "Asian", "Other", "Pending", "Am.Indian/Alaskan", "Hawaiian/Pacific")))
# US AND STATE DATA ####
# us_data <- tigris::states(cb = TRUE)
# us_data <- sf::st_as_sf(us_data)
# names(us_data) <- tolower(names(us_data))
# us_data <-
# us_data %>%
# filter(stusps %in% state.abb)
# saveRDS(us_data, here::here("data_rds/us_data.rds"))
# us_data <- readRDS(here::here("data_rds/us_data.rds"))
#
#
#
# # state_pop <- readxl::read_excel(here::here("state_populations.xlsx"))
# # state_pop$region <- gsub("\\.", "", state_pop$region)
# # state_pop <-
# # state_pop %>%
# # pivot_longer(cols = -c(region, census_2010),
# # names_to = "pop_year",
# # names_prefix = "pop_",
# # values_to = "pop_estimate")
# # saveRDS(state_pop, here::here("data_rds/state_pop.rds"))
#
# state_pop <- readRDS( here::here("data_rds/state_pop.rds"))
# state_pop_2019 <-
# state_pop %>%
# filter(pop_year == 2019) %>%
# filter(region %in% state.name)
# US DATA FROM THE COVIDTRACKING PROJECT ####
# COVIDTRACKING PROJECT QUIT TRACKING DATA ON MARCH 7, 2021;
# DATA SOURCE CHANGED TO CDD; TEST POSITIVITY RATES NOT AVAILABLE THERE
# covus_currentdata <- readr::read_csv(here::here("data_covidtracking/daily.csv"))
# covus_currentdata <-
# covus_currentdata %>%
# filter(state %in% state.abb,
# !state %in% c("AK", "HI") ) %>%
# mutate(date = lubridate::ymd(date))
# covus_currentdata <- readRDS(here::here("data_covidtracking/covus_current.rds"))
# covus_currentdata <-
# covus_currentdata %>%
# left_join(us_data,
# by = c("state" = "stusps")) %>%
# left_join(state_pop_2019, by = c("name" = "region")) %>%
# mutate(test_density_pct = total/pop_estimate*100,
# case_density_pct = positive/pop_estimate*100)
#
# covus_currentdata <-
# covus_currentdata %>%
# group_by(state) %>%
# arrange(date) %>%
# mutate(new_cases_state_7davg = round(RcppRoll::roll_meanr(positiveIncrease, 7), 1),
# new_deaths_state_7davg = round(RcppRoll::roll_meanr(deathIncrease, 7), 1),
# new_deaths_state_7davg_100k = new_deaths_state_7davg/pop_estimate*100000,
# new_testpos_pct = round(positiveIncrease/totalTestResultsIncrease*100, 1),
# new_testpos_pct_7davg = round(RcppRoll::roll_meanr(positiveIncrease, 7)/RcppRoll::roll_meanr(totalTestResultsIncrease, 7)*100, 1),
# new_cases_state_7davg_100k = new_cases_state_7davg/pop_estimate*100000
# ) %>%
# ungroup()
#
# covus_currentdata <-
# covus_currentdata %>%
# group_by(date) %>%
# mutate(rank_newcases7davg = rank(-new_cases_state_7davg_100k)) %>%
# ungroup()
#
# state_tbl <-
# tibble(state_name = state.name,
# state = state.abb)
#
# covus_currentdata_2163 <-
# covus_currentdata %>%
# mutate(geometry = st_transform(geometry, 2163)) %>%
# left_join(state_tbl, by = "state") %>%
# mutate(tt = sprintf("<tr><td><b><i>%s (%s)</i></b></td></tr>
# <tr><td>New Cases/100k (7d Avg): %s</td></tr>
# <tr><td>Test %% Positive (7d Avg): %1.1f%%</td></tr>
# <tr><td>State Rank (Cases/100k 7d Avg): %s</td></tr>
# <tr><td>Total Cases: %s</td></tr>
# <tr><td>Total Deaths: %s</td></tr>
# <tr><td>New Deaths/d/100k (7d Avg): %1.2f</td></tr>
# <tr><td>CFR: %1.2f%%</td></tr>
# <tr><td>Case Density: %1.1f%% of pop.</td></tr>
# <tr><td>Testing Density: %1.1f%% of pop.</td></tr>",
# toupper(state_name), paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
# comma(new_cases_state_7davg_100k, digits = 1), new_testpos_pct_7davg, rank_newcases7davg, comma(positive, digits = 0),
# comma(death, digits = 0), new_deaths_state_7davg_100k, death/positive*100, case_density_pct, test_density_pct))
# US DATA FROM THE CDC ####
# NOTE: COVIDTRACKING PROJECT QUIT TRACKING DATA ON MARCH 7, 2021;
# TEST POSITIVITY RATES NOT AVAILABLE THERE
# CDC DATA
# cdc data download not available or moved and I can't find it 7/2023
# cdc_currentdata <- readr::read_csv(here::here("data_cdc/cdc_state_info.csv"))
# names(state.name) <- state.abb
# cdc_currentdata <-
# cdc_currentdata %>%
# filter(state %in% state.abb) %>%
# # !state %in% c("AK", "HI") ) %>%
# mutate(date = lubridate::mdy(submission_date),
# region = state.name[state])
#
#
# cdc_currentdata <-
# cdc_currentdata %>%
# # left_join(us_data,
# # by = c("state" = "stusps")) %>%
# left_join(state_pop_2019, by = "region")
#
#
# cdc_currentdata <-
# cdc_currentdata %>%
# group_by(state) %>%
# arrange(date) %>%
# mutate(case_density_pct = tot_cases/pop_estimate*100,
# new_cases_state_7davg = round(RcppRoll::roll_meanr(new_case, 7), 1),
# new_deaths_state_7davg = round(RcppRoll::roll_meanr(new_death, 7), 1),
# new_deaths_state_7davg_100k = new_deaths_state_7davg/pop_estimate*100000,
# new_cases_state_7davg_100k = new_cases_state_7davg/pop_estimate*100000
# ) %>%
# ungroup()
#
# cdc_currentdata <-
# cdc_currentdata %>%
# group_by(date) %>%
# mutate(rank_newcases7davg = rank(-new_cases_state_7davg_100k)) %>%
# ungroup()
#
#
#
# state_tbl <-
# tibble(state_name = state.name,
# state = state.abb)
#
#
# cdc_currentdata_2163 <-
# cdc_currentdata %>%
# filter(date == max(date)) %>%
# left_join(us_data, by = c("state" = "stusps")) %>%
# mutate(geometry = st_transform(geometry, 2163)) %>%
# left_join(state_tbl, by = c("state", "state_name")) %>%
# mutate(tt = sprintf("<tr><td><b><i>%s (%s)</i></b></td></tr><br>
# <tr><td>New Cases/100k (7d Avg): %s</td></tr><br>
# <tr><td>State Rank (Cases/100k 7d Avg): %s</td></tr><br>
# <tr><td>Total Cases: %s</td></tr><br>
# <tr><td>Total Deaths: %s</td></tr><br>
# <tr><td>New Deaths/d/100k (7d Avg): %1.2f</td></tr><br>
# <tr><td>CFR: %1.2f%%</td></tr><br>
# <tr><td>Case Density: %1.1f%% of pop.</td></tr>",
# toupper(state_name), paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
# comma(new_cases_state_7davg_100k, digits = 1), rank_newcases7davg, comma(tot_cases, digits = 0),
# comma(tot_death, digits = 0), new_deaths_state_7davg_100k, tot_death/tot_cases*100, case_density_pct))
# mutate(tt = sprintf("<tr><td><b><i>%s (%s)</i></b></td></tr>
# <tr><td>New Cases/100k (7d Avg): %s</td></tr>
# <tr><td>State Rank (Cases/100k 7d Avg): %s</td></tr>
# <tr><td>Total Cases: %s</td></tr>
# <tr><td>Total Deaths: %s</td></tr>
# <tr><td>New Deaths/d/100k (7d Avg): %1.2f</td></tr>
# <tr><td>CFR: %1.2f%%</td></tr>
# <tr><td>Case Density: %1.1f%% of pop.</td></tr>",
# toupper(state_name), paste0(month.name[lubridate::month(date)], " ", lubridate::mday(date)),
# comma(new_cases_state_7davg_100k, digits = 1), rank_newcases7davg, comma(tot_cases, digits = 0),
# comma(tot_death, digits = 0), new_deaths_state_7davg_100k, tot_death/tot_cases*100, case_density_pct))
# county
tndoh_vax_county <-
readxl::read_excel(here::here("data_tndoh/COVID_VACCINE_COUNTY_SUMMARY.XLSX"),
col_types = c("date", "text", rep("numeric", 8)))
names(tndoh_vax_county) <- tolower(names(tndoh_vax_county))
tndoh_vax_county <-
tndoh_vax_county %>%
mutate(date = lubridate::ymd(date),
county_lower = tolower(county)) %>%
left_join(TN_counties,
by = "county_lower") %>%
rowwise() %>%
mutate(vaccine_count_pct = vaccine_count/population*100,
recipient_count_pct = recipient_count/population*100) %>%
ungroup() %>%
filter(!is.na(statefp))
latest_vax_date <- max(tndoh_vax_county$date, na.rm = T)
tndoh_vax_county_current <-
left_join(
{tndoh_vax_county %>%
filter(date == latest_vax_date,
!is.na(statefp))},
{tndoh_countynew %>%
filter(date == max(date)) %>%
select(county_lower, total_cases_per100kpop, new_cases_per100kpop_7davg)}
) %>%
mutate(immune_pct = recipient_count_pct + total_cases_per100kpop/1000)
# CDC WASTEWATER DATA (added Oct 4, 2022) ----
nwss_dat <- RSocrata::read.socrata(url = "https://data.cdc.gov/resource/g653-rqe2.csv")
# nwss_dat <- content(httr::GET("https://data.cdc.gov/resource/g653-rqe2.csv")) # get data
wwtp_id_details_tn <- readRDS(here::here("data_rds/wwtp_id_details_tn.rds"))
nwss_dat <-
nwss_dat %>%
mutate(original_kpid = key_plot_id) %>%
extract(key_plot_id,
into = c("nwss_source", "state", "wwtp_id", "sample_location", "sample_location_specify", "sample_type"),
regex = "(CDC_BIOBOT|CDC_LUMINULTRA|NWSS)_([a-z]+)_([0-9]+)_(Before treatment plant|Treatment plant)_([0-9]+)*_*(.*)") %>%
rename(key_plot_id = original_kpid) %>%
mutate(wwtp_id = as.numeric(wwtp_id), # required for merging
date = lubridate::ymd(date))
nwss_dat_tn <-
nwss_dat %>%
filter(state == "tn") %>%
left_join(wwtp_id_details_tn, by = "wwtp_id")
# CDC VARIANT DATA
## ├ Get Data ----
# data source: https://covid.cdc.gov/covid-data-tracker/#variant-proportions
# download.file("https://data.cdc.gov/api/views/jr58-6ysp/rows.csv?accessType=DOWNLOAD",
# destfile = here::here("data_cdc/cdc_variant_nowcast.csv")) # contains NOW-CAST results- all US regions
cdc_variants_data <- readr::read_csv(here::here("data_cdc/cdc_variant_nowcast.csv"))
cdc_variants_data$usa_or_hhsregion <- factor(cdc_variants_data$usa_or_hhsregion,
levels = c("USA", 1:10),
labels = c("USA", paste("Region", 1:10)))
cdc_variants_data$variant <- factor(cdc_variants_data$variant)
cdc_variants_data$week_ending <- lubridate::date(lubridate::mdy_hms(cdc_variants_data$week_ending))
cdc_variants_data$published_date <- lubridate::date(lubridate::mdy_hms(cdc_variants_data$published_date))
Interactive Figure: mouse/hover over for details.
pl_new_tnsdc <-
tndoh_daily %>%
filter(date >= max(date) - 120) %>%
pivot_longer(cols = c(new_cases, new_hosp, new_deaths_by_dod),
names_to = "var",
values_to = "value") %>%
# filter(date > lubridate::ymd("2021-01-01")) %>%
ggplot(aes(x = date)) +
geom_rect(aes(xmin = max(tndoh_daily$date)-10, xmax = max(tndoh_daily$date) + 1, ymin = -Inf, ymax = Inf),
fill = "grey", alpha = 0.2) +
ggiraph::geom_col_interactive(aes(y = value, fill = var, tooltip = tt_daily, data_id = date), alpha = 0.5,
position = "identity") +
geom_line(aes(y = new_cases_trend.stl),color = "red", size = 1) +
geom_line(aes(y = new_cases_7davg), color = "red", size = 0.5, linetype=2) +
labs(x=NULL,
y="Results by day reported",
subtitle = "Results by Day, New (Total)",
caption = "Solid line is time series trend. Dotted line is 7d avg\nResults within the past 7-10 days (shaded period) may not be fully reported",
fill = NULL) +
scale_fill_manual(values = c("steelblue",
"red",
"orange"),
labels = c(paste0("Cases: ", comma(tndoh_daily$new_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
" (", comma(tndoh_daily$total_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
paste0("Deaths: ", comma(tndoh_daily$new_deaths_by_dod[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
" (", comma(tndoh_daily$total_deaths[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
paste0("Hospitalizations: ", comma(tndoh_daily$new_hosp[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
" (", comma(tndoh_daily$total_hosp[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")")
)) +
scale_y_continuous(expand = c(0,0),
sec.axis = sec_axis( trans = ~./TN_pop2019*10^5, name="Cases/100k Population")) +
scale_x_date(date_breaks = "1 month", date_labels = "%b" ) +
theme(
legend.position = c(.01, 1),
legend.justification = c(0, 1),
legend.text = element_text(size = 10))
hover_css <- "background-color:orange;stroke:black;fill:orange;"
tooltip_css <- " background-color:rgba(255, 255, 255, 0.9);color:black;stroke:black;opacity:0.1;border:2px solid black;border-radius: 8px;font-family:Arial;font-size: 0.875em;"
ggiraph::girafe(ggobj = pl_new_tnsdc,
width_svg = 8,
options = list(
offx = 40, offy = 40,
opts_hover(css = hover_css),
opts_tooltip(css = tooltip_css)
) )
Deaths plotted by date-of-death. Reporting lags, so deaths within last 1-3 months may not be listed.
Limited to Last 4 months for better detail. New case trend statistic and line estimated by Seasonal and Trend decomposition using Loess (STL) decomposition to account for daily variation- results are similar to a 7-day average but better account for the cyclical variation.
# PLOTS ####
pl_tests_new <-
tndoh_daily %>%
# filter(date - 90) %>%
ggplot(aes(x=date)) +
geom_col(aes(y = new_tests, fill = "Tests"), alpha = 0.5) +
geom_col(aes(y = new_cases, fill = "Cases"), alpha = 0.7) +
labs(x=NULL,
fill = NULL,
y="Tests by day reported",
subtitle = "Tests by Day, New (Total)") +
scale_fill_manual(values = c("Cases" = "steelblue",
"Tests" = "grey"),
labels = c(paste0("Cases:", comma(tndoh_daily$new_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
" (", comma(tndoh_daily$total_cases[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"),
paste0("Tests:", comma(tndoh_daily$new_tests[tndoh_daily$date == max(tndoh_daily$date)], digits = 0),
" (", comma(tndoh_daily$total_tests[tndoh_daily$date == max(tndoh_daily$date)], digits = 0), ")"))) +
scale_y_continuous(expand = c(0,0),
sec.axis = sec_axis( trans=~./TN_pop2019*10^5, name="per 100k Population")) +
scale_x_date(date_breaks = "3 months", date_labels = "%b-%Y" ) +
annotate(geom = "text", x=min(tndoh_daily$date)+1, y = 0.65*max(tndoh_daily$new_tests, na.rm = T), hjust = 0,
label = label_pos_pct, size = 3.5) +
theme(#axis.text.x = element_text(angle = 45, hjust = 1),
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = c(.01, 1),
legend.justification = c(0, 1),
legend.text = element_text(size = 10))
# age ####
pl_age <-
tndoh_age %>%
filter(date == max(date)) %>%
ggplot() +
geom_blank(data = tndoh_age[tndoh_age$ar_totalcases == max(tndoh_age$ar_totalcases), ],
aes(x=age.f, y=ar_totalcases*1.1)) +
geom_bar(aes(age.f, ar_totalcases), alpha = 0.5, fill = "steelblue", stat = "identity", ) +
geom_bar(aes(age.f, ar_totaldeaths), alpha = 0.5, fill = "red", stat = "identity", ) +
geom_text(aes(age.f, ar_totalcases*1.02, vjust =0, label = comma(ar_totalcases, digits = 0)),
color = "steelblue", size = 3) +
geom_text(aes(age.f, ar_totaldeaths*1.01, vjust =0, label = comma(ar_totaldeaths, digits = 0)),
color = "red", size = 3) +
labs(x=NULL,
y="Total Case Count",
subtitle = "Total Cases by Age") +
scale_y_continuous(expand = c(0,0)) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# multi-axis plot
pl_age_death <-
tndoh_age %>%
filter(date == max(date)) %>%
ggplot(aes(age.f, deaths_case_rate, group = date)) +
# pl_age_death <-
# pl_age_death +
geom_point(color = "red") +
geom_line(color = "red")
# Dual axis
ggp <- ggplot_build(pl_age)
ggp2 <- ggplot_build(pl_age_death)
ylim.prim <- ggp$layout$panel_params[[1]]$y.range # in this example, precipitation
ylim.sec <- ggp2$layout$panel_params[[1]]$y.range # in this example, temperatureb <- diff(ylim.prim)/diff(ylim.sec)
b <- diff(ylim.prim)/diff(ylim.sec)
a <- b*(ylim.prim[1] - ylim.sec[1])
pl_age_death_final <-
tndoh_age %>%
filter(date == max(date)) %>%
ggplot(aes(x=age.f)) +
geom_blank(data = tndoh_age[tndoh_age$ar_totalcases == max(tndoh_age$ar_totalcases), ],
aes(x=age.f, y=ar_totalcases*1.1)) +
geom_bar(aes(age.f, ar_totalcases), alpha = 0.5, fill = "steelblue", stat = "identity", ) +
geom_bar(aes(age.f, ar_totaldeaths), fill = "red", stat = "identity", ) +
geom_line(aes(y = a + deaths_case_rate*b, group = date), color = "red", size = 1.5) +
geom_point(aes(y = a + deaths_case_rate*b, group = date), color = "red", size = 2, shape = 21, fill = "white") +
geom_text(aes(age.f, ar_totalcases*1.02, vjust =0, label = comma(ar_totalcases, digits = 0)),
color = "steelblue", size = 3) +
geom_text(aes(age.f, ar_totaldeaths*1.01, vjust =0, label = comma(ar_totaldeaths, digits = 0)),
color = "red", size = 3) +
labs(x=NULL,
y="Total Case Count",
subtitle = "Total Cases by Age") +
scale_y_continuous(expand = c(0,0),
sec.axis = sec_axis(~ (. - a)/b, name = "Crude Case Fatality (%)")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
axis.title.y.right = element_text(color = "red"),
axis.text.y.right = element_text(color = "red"))
race_total <- sum(tndoh_demographics$cases_byrace[tndoh_demographics$race != "Pending"])
# AGE HEATMAP
max_newcaserate_byage <- max(tndoh_age$ar_new_cases_100k_7davg, na.rm = T)
ar_newcase_heatmap_100k <-
tndoh_age %>%
filter(date >= lubridate::ymd("2020-04-01")) %>%
ggplot(aes(x = date,
y = age.f)) +
geom_tile(aes(fill = ar_new_cases_100k_7davg,
color = ar_new_cases_100k_7davg)) +
viridis::scale_fill_viridis(labels = scales::comma,
limits = c(0, max_newcaserate_byage), oob = scales::squish,
option = "plasma", na.value="black") +
viridis::scale_color_viridis(labels = scales::comma,
limits = c(0, max_newcaserate_byage), oob = scales::squish,
option = "plasma", na.value="black") +
labs(fill = "New Cases/100k/d\n(7dAvg)",
color = "New Cases/100k/d\n(7dAvg)",
x = NULL, y = NULL,
subtitle = "New Case Rate, by Age (population adjusted)") +
scale_x_date(date_breaks = "3 months", date_labels = "%b-%Y", expand = c(0, 0)) +
theme_minimal() +
theme(axis.text.y = element_text(size = 12),
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid = element_blank(),
legend.position = "bottom",
legend.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size = 10),
legend.text = element_text(size = 8))
# fig 1 ####
pl_TOP <- cowplot::plot_grid(pl_new_tnsdc, pl_tests_new,
pl_age_death_final, ar_newcase_heatmap_100k,
align = "h", axis = "bt",
nrow = 2)
cowplot::plot_grid(title_tnsdc, pl_TOP, nrow = 2, rel_heights = c(.1, 1))
# ggsave(here::here("figures/pl_tnmap.jpg"), dpi = 600, width = 12, height = 4, units = "in")
SUMMARY BY AGE In Tennessee, 4.22% of residents >80 and 1.76% of those age 71-80 have died of COVID-19, to date.
Age Group | Population (n) | Cases | Cases(% of Total) | Deaths | CFR% | Deaths(% of Total) | Deaths/100k | Deaths (% in Age) | Given 1st Vaccine (% in Age) |
---|---|---|---|---|---|---|---|---|---|
0-10 | 908244 | 241111 | 9% | 18 | 0.007 | 0% | 1.98 | 0.00% | 10.9% |
11-20 | 858307 | 346009 | 13% | 28 | 0.008 | 0% | 3.26 | 0.00% | 47.3% |
21-30 | 943700 | 440444 | 17% | 222 | 0.050 | 1% | 23.52 | 0.02% | 58.2% |
31-40 | 879850 | 401487 | 15% | 599 | 0.149 | 2% | 68.08 | 0.07% | 65.1% |
41-50 | 854356 | 364574 | 14% | 1528 | 0.419 | 5% | 178.85 | 0.18% | 68.9% |
51-60 | 893802 | 346939 | 13% | 3558 | 1.026 | 12% | 398.07 | 0.40% | 75.8% |
61-70 | 782810 | 260806 | 10% | 6207 | 2.380 | 21% | 792.91 | 0.79% | 87.2% |
71-80 | 481538 | 163162 | 6% | 8465 | 5.188 | 28% | 1757.91 | 1.76% | 93.9% |
81+ | 226567 | 80898 | 3% | 9564 | 11.822 | 32% | 4221.27 | 4.22% | 84.7% |
Population within each age group taken from US Census estimates for
TN, 2019CFR=Case Fatality Rate (%) Vaccination data last updated on 2022-08-22. |
NOTES:
A bit of good news. Vaccination data is not updated every day- most recent data provided.
TN now provides detail on number of recipients and number of vaccinations.
Age Group | Given 1st Vaccine (n) | Estimated Pop. (n) | Given 1st Vaccine (% in Age) |
---|---|---|---|
<5 | 10900 | 408605 | 2.67% |
5-11 | 119258 | 585900 | 20.35% |
12-15 | 156446 | 343637 | 45.53% |
16-20 | 230052 | 428409 | 53.70% |
21-30 | 548940 | 943700 | 58.17% |
31-40 | 573213 | 879850 | 65.15% |
41-50 | 589057 | 854356 | 68.95% |
51-60 | 677701 | 893802 | 75.82% |
61-70 | 682573 | 782810 | 87.20% |
71-80 | 452299 | 481538 | 93.93% |
81+ | 191876 | 226567 | 84.69% |
Population within each age group estimated from US Census estimates for TN, 2019Vaccination data last updated on 2022-08-22. |
most recent County Vaccination data update on 2022-08-22
# max(tndoh_vax_county_current$date)
plot_vacc <-
tndoh_vax_county_current %>%
filter(date == latest_vax_date,
!is.na(statefp)) %>%
ggplot() +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = recipient_count_pct),
color = "grey30",
inherit.aes = F) +
geom_sf_text(aes(label=sprintf("%0.1f%%", recipient_count_pct),
# geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos),
geometry = geometry,
color = cut(recipient_count_pct, 2)),
size = 3) +
viridis::scale_fill_viridis(option = "C", direction = 1) +
scale_color_manual(values = c("white", "black"), guide = F) +
labs(fill = "%Pop.:",
title = paste0("Vaccination Coverage (latest update on ", latest_vax_date, ")"),
subtitle = sprintf("First Vaccinations: %s (%1.1f%% of TN Population)\nCompleted Vacc. : %s (%1.1f%% of Pop.)\nTotal Vaccinations: %s",
comma(max(tndoh_vax_state$recipient_count, na.rm = T), 0),
max(tndoh_vax_state$recipient_count_pct, na.rm = T),
comma(max(tndoh_vax_state$recip_fully_vacc, na.rm = T), 0),
max(tndoh_vax_state$recip_fully_vacc_pct, na.rm = T),
comma(max(tndoh_vax_state$vaccine_count, na.rm = T), 0))
) +
theme_void() +
theme(legend.position = "left")
# plot_immune <-
# tndoh_vax_county_current %>%
# filter(date == latest_vax_date,
# !is.na(statefp)) %>%
# ggplot() +
# geom_sf(data = TN_counties) +
# geom_sf(aes(geometry = geometry,
# fill = immune_pct),
# color = "grey30",
# inherit.aes = F) +
# geom_sf_text(aes(label=sprintf("%0.1f%%", immune_pct),
# # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos),
# geometry = geometry,
# color = cut(immune_pct, 2)),
# size = 3) +
# viridis::scale_fill_viridis(option = "C", direction = 1) +
# scale_color_manual(values = c("white", "black"), guide = F) +
# labs(fill = "%Pop.:",
# title = "Vaccination + Documented Prior Infection") +
# theme_void() +
# theme(legend.position = "left")
plot_vacc
# cowplot::plot_grid(plot_vacc, plot_immune,
# nrow = 2)
most recent County data update on 2023-12-09
note: The “New Hospitalizations” in this chart, as reported by the state, may not accurately reflect actual hospitalizations. Right y-axis is population adjusted (per 100k for TN). Ignore for “Test Positivity.”
Counties in the top 10 highest new case rate over the past 7 days; 7-day average trend-line for daily new cases
# add PCT-Pos trend line calculated by STL (Seasonal and Trend decomposition using Loess)
tndoh_daily$pos_pct_today <- replace_na(tndoh_daily$pos_pct_today, 0)
tndoh_daily$pos_pct_today_tl <- ts(tndoh_daily$pos_pct_today*100, frequency = 7, start = 1)
tndoh_daily$pos_pct_today.stl <- stl(tndoh_daily$pos_pct_today_tl, s.window = 7)$time.series[, "trend"]
# add Vaccine trend line calculated by STL (Seasonal and Trend decomposition using Loess)
tndoh_daily$new_vaccine_count <- replace_na(tndoh_daily$new_vaccine_count, 0)
tndoh_daily$new_vaccine_count_tl <- ts(tndoh_daily$new_vaccine_count, frequency = 7, start = 1)
tndoh_daily$new_vaccine_count.stl <- stl(tndoh_daily$new_vaccine_count_tl, s.window = 7)$time.series[, "trend"]
variable.labs <- c( "Cases", "Cases",
"Hospitalizations",
"Deaths", "Deaths", "Deaths",
"Tests", "Test Positivity (%)", "Test Positivity (%)",
"Vaccinations", "Vaccinations")
names(variable.labs) <- c( "new_cases_trend.stl", "new_cases_trend.stl",
"new_hosp",
"new_deaths","new_deaths_by_dod", "new_deaths_trend.stl",
"new_tests", "pos_pct_today", "pos_pct_today.stl",
"new_vaccine_count", "new_vaccine_count.stl")
plot_state_trends <-
tndoh_daily %>%
mutate(pos_pct_today = pos_pct_today*100) %>%
select(date, new_cases_trend.stl, new_deaths_trend.stl, pos_pct_today.stl, new_vaccine_count.stl) %>%
pivot_longer(cols = -date,
names_to = "variable",
values_to = "value") %>%
group_by(variable) %>%
arrange(date) %>%
ggplot() +
geom_line(aes(date, value), size = 1, color = "red") +
labs(x=NULL,
y="Decomposition 7-day Average",
title = "TN State Trends") +
facet_wrap(~variable, scales = "free_y", nrow = 1,
labeller = labeller(variable = variable.labs)) +
theme(strip.text = element_text(size = 14))
# plot_state_trends <-
# tndoh_daily %>%
# mutate(pos_pct_today = pos_pct_today*100) %>%
# # select(date, new_cases_7davg)
# # select(date, new_cases, new_deaths_by_dod, new_tests, pos_pct_today) %>%
# select(date, new_cases_trend.stl, new_deaths_by_dod, pos_pct_today, new_vaccine_count) %>%
# # select(date, new_cases, new_hosp, new_deaths_by_dod, new_tests, pos_pct_today, new_vaccine_count) %>%
# pivot_longer(cols = -date,
# names_to = "variable",
# values_to = "value") %>%
# group_by(variable) %>%
# arrange(date) %>%
# mutate(value_7davg = RcppRoll::roll_meanr(value, 7)) %>%
# ggplot() +
# # plot_state_trends <-
# # plot_state_trends +
# geom_line(aes(date, value), size = 0.6, color = "grey") +
# geom_line(aes(date, value_7davg), size = 1, color = "red") +
# # scale_y_continuous(sec.axis = sec_axis( trans=~./TN_pop2019*10^5, name="per 100k Population")) +
# labs(x=NULL,
# y="Rolling 7-day Average",
# title = "TN State Trends") +
# facet_wrap(~variable, scales = "free_y", nrow = 1,
# labeller = labeller(variable = variable.labs)) +
# theme(strip.text = element_text(size = 14))
# MIDTN counties plot ####
counties_top10 <-
tndoh_countynew %>% #select(new_cases_county_7davg)
ungroup() %>%
filter(date == max(date, na.rm = T),
name %in% TN_counties$name) %>% # exclude the pending and out-of-state
top_n(10, wt = new_cases_county_7davg) %>%
pull(name)
plot_top10_counties <-
tndoh_countynew %>%
filter(name %in% counties_top10) %>%
ggplot() +
# plot_top10_counties <-
# plot_top10_counties +
# geom_bar(stat = "identity") +
geom_area(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
geom_col(aes(date, new_deaths), fill = "red") +
# geom_line(aes(date, new_cases_county_trend), color = "steelblue", size = 1) +
geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
labs(x=NULL,
y="New Results by day",
fill = NULL,
alpha = NULL,
title = "Top 10 New Cases by Day (7d avg.)") +
scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
theme(#axis.text.x = element_text(angle = 90, hjust = 1),
legend.position = "none") +
facet_wrap(~name, nrow = 2,
scales = "free_y")+
theme(strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1))
cowplot::plot_grid(plot_state_trends,
plot_top10_counties,
nrow = 2, rel_heights = c(1.25,2))
deaths plotted by actual date-of-death. reporting typically lags several weeks, so most recent dates will not be fully reported.
The Variant distribution is reported by 10 US Regions. Only data for the US total and Region 4 (SouthEast) is shown here. CDC NowCast data is used for this report.
pl_variants <-
cdc_variants_data %>%
filter(usa_or_hhsregion == "Region 4" | usa_or_hhsregion == "USA",
week_ending > lubridate::ymd(Sys.Date())-365,
published_date == max(published_date)) %>%
ggplot(aes(week_ending, share,
group = variant,
fill = variant)) +
ggiraph::geom_area_interactive(color = "grey50", alpha = 0.6, aes(tooltip = variant), hover_css = "stroke:red;stroke-width:inherit;") +
scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
scale_x_date(expand = c(0,0), date_breaks = "3 month", date_labels = "%b-%Y" ) +
theme(legend.position = "none") +
labs(x = NULL,
y = "Proportion of cases by Variant",
title = "Variant Distribution in Southeastern US (CDC Region 4)",
subtitle = "Hover over with mouse for details",
fill = "Variant:") +
facet_wrap(~usa_or_hhsregion)
ggiraph(code = {print(pl_variants)}, zoom_max = 1, width = .8, width_svg = 7, height_svg = 4)
Sites in TN recently started reporting data to the National Wastewater Surveillance System (NWSS). (Plots added here Oct 2022).
nwss_dat_tn %>%
filter(!is.na(pcr_conc_smoothed)) %>%
ggplot(aes(date, pcr_conc_smoothed/1000000, group = wwtp_id )) +
geom_line() +
facet_wrap(~county_names, nrow = 1, labeller = labeller(county_names = label_wrap_gen(7))) +
labs(title = "TN Wastewater Surveillance",
x= NULL, y = "SARS-COV2 Conc by PCR\n(Flow-Population normalized)",
caption = "Source: CDC, National Wastewater Surveillance System (NWSS)") +
scale_y_continuous(label=comma, ) +
scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
theme(strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1)) +
theme(legend.position = "none")
Counties with the top 10 greatest weekly increase in New Cases, by % Change of 7-day average versus 7 days prior (limited to counties with average of >5 daily cases). Case rates are currently increasing in 44 counties.
plot_top10_increase <-
tndoh_countynew %>%
filter(name %in% counties_top10_increase) %>%
ggplot() +
# plot_top10_increase +
# geom_bar(stat = "identity") +
geom_col(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
geom_col(aes(date, new_deaths), fill = "red") +
geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
labs(x=NULL,
y="New Results by day",
fill = NULL,
alpha = NULL,
subtitle = "New Cases by Day: 10 Counties with the largest 1-week increase") +
scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
theme(#axis.text.x = element_text(angle = 90, hjust = 1),
legend.position = "none") +
facet_wrap(~name, nrow = 2,
scales = "free_y")+
theme(strip.text = element_text(size = 14),
axis.text.x = element_text(angle = 45, hjust = 1))
plot_top10_increase
# print("No increase in any counties")
# *currently disabled due to drop in cases*
counties_near_peak_cases <-
tndoh_countynew %>%
ungroup() %>%
filter(new_cases_county_trend_ratio >= 0.9,
date >= Sys.Date() - 14,
!county %in% c("Pending", "Out of State")) %>%
group_by(county) %>%
filter(new_cases_county_trend_ratio == max(new_cases_county_trend_ratio)) %>%
ungroup()
counties_near_peak_cases_top10 <-
counties_near_peak_cases %>%
slice_max(order_by = new_cases_county_trend, n= 10) %>%
pull(county)
No counties are currently near (at or within >90% of) all-time peak new case rate.
Counties with the top 10 greatest weekly Decrease in New Cases, by % Change of 7-day average versus 7 days prior (limited to counties with average of >5 daily cases).
Only counties with a decline in case averages are shown (may be
<10).
Cases are currently declining in 43 counties
# MIDTN counties plot ####
pl_best10 <-
tndoh_countynew %>%
filter(name %in% counties_top10_decrease) %>%
ggplot() +
geom_col(aes(date, new_cases), fill = "steelblue", alpha = 0.5) +
geom_col(aes(date, new_deaths), fill = "red") +
geom_line(aes(date, new_cases_county_7davg), color = "steelblue", size = 1) +
labs(x=NULL,
y="New Results by day",
fill = NULL,
alpha = NULL,
subtitle = "New Cases by Day: 10 Counties with the largest 1-week Decrease") +
ylim(0,NA) +
scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
theme(#axis.text.x = element_text(angle = 90, hjust = 1),
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text = element_text(size = 14)) +
facet_wrap(~name, nrow = 2,
scales = "free_y")
pl_best10
# cat("There were no counties with a decreased 7d average")
# ifelse(length(counties_top10_decrease) >= 1,
# # IF >=1 with negative trend
# plot(pl_best10),
# # IF No counties with negative trend
# "There are no counties with a decreasing case trend")
County color indicates the New Cases per day per county population (7 day average)
# make county data to list
county_plots <-
tndoh_countynew %>%
filter(!is.na(date),
date >= lubridate::ymd("2020-03-15")) %>%
group_by(name) %>%
nest() %>%
arrange(name) %>%
mutate(pl_new = map2(data, name,
~ggplot(data = .x) +
geom_col(aes(date, new_cases, fill = "Cases"), alpha = 0.5) +
geom_col(aes(date, new_deaths, fill = "Deaths")) +
geom_line(aes(date, new_cases_county_7davg, fill = "Case 7-Day Avg"), color = "steelblue", size = 1) +
# geom_col(aes(date, new_hospitalized, fill = "Hospitalized")) +
labs(x=NULL,
y="New Results by day",
title = glue::glue("{.y} County"),
fill = NULL,
alpha = NULL,
subtitle = "New Results by Date Reported") +
scale_fill_manual(values = c("Cases" = "steelblue",
"Case 7-Day Avg" = "steelblue",
"Deaths" = "red"),
# "Hospitalized" = "orange"),
labels = c(paste0("Total Cases: ", comma(.x$total_cases[.x$date == max(.x$date)], digits = 0)),
paste0("7d Avg: ", comma(.x$new_cases_county_7davg[.x$date == max(.x$date)], digits = 0)),
paste0("Total Deaths: ", comma(.x$total_deaths[.x$date == max(.x$date)], digits = 0))
# paste0("Hospitalized: ", .x$total_hospitalized[.x$date == max(.x$date)])
)
) +
scale_y_continuous(expand = c(0,0), limits = c(0,NA)) +
# ylim(0,NA) +
scale_x_date(date_breaks = "3 month", date_labels = "%b-%Y" ) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "bottom")
))
county_data_byday_now <-
left_join(TN_counties,
{tndoh_countynew %>%
filter(date == max(date, na.rm = T)) %>%
select(name, total_cases, new_cases_county_7davg, new_cases_per100kpop_7davg, new_cases_per100kpop_7davg_rank,
new_testpos_pct, new_testpos_pct_7davg, new_tests, new_pos_tests,
new_cases_rank, new_cases, new_cases_rank, new_cases_county_trend_ratio,
total_deaths, total_hospitalized)}) %>%
mutate(cases = replace_na(total_cases, 0),
deaths = replace_na(total_deaths, 0)) %>%
left_join({county_plots %>% select(name, pl_new)},
by = "name")
lbls <-
sprintf("<strong><i>%s Co.</i></strong><br/>
Total Cases: %s<br/>
New Cases/d/100k, 7dAvg (rank): %1.1f (%s/%s)<br/>
New Cases Today (rank): %s (%s/%s)<br/>
Test +%% Today: %1.1f%% (%s/%s)<br/>
Test +%% (7dAvg): %1.1f%%<br/>
Total Deaths: %s<br/>
Current Case trend vs. Max: %1.1f%%<br/>
<i>Click to view time-course</i>",
toupper(county_data_byday_now$name),
comma(county_data_byday_now$total_cases, digits = 0),
county_data_byday_now$new_cases_per100kpop_7davg, county_data_byday_now$new_cases_per100kpop_7davg_rank, max(county_data_byday_now$new_cases_per100kpop_7davg_rank),
comma(county_data_byday_now$new_cases, digits = 0),county_data_byday_now$new_cases_rank, max(county_data_byday_now$new_cases_rank),
county_data_byday_now$new_testpos_pct, county_data_byday_now$new_pos_tests, county_data_byday_now$new_tests,
county_data_byday_now$new_testpos_pct_7davg,
county_data_byday_now$total_deaths,
county_data_byday_now$new_cases_county_trend_ratio*100
) %>%
lapply(htmltools::HTML)
# pal <- colorQuantile("viridis", domain = county_data_byday_now$new_cases_per100kpop_7davg, n = 10, reverse = FALSE)
pal <- colorNumeric("viridis", domain = county_data_byday_now$new_cases_per100kpop_7davg, reverse = FALSE)
leaflet() %>%
addProviderTiles(providers$Stamen.TonerLite) %>%
addPolygons(data = TN_counties, # borders of all counties
color = "grey",
fill = NA,
weight = 1) %>%
addPolygons(data = county_data_byday_now,
group = "name",
fillColor = ~pal(county_data_byday_now$new_cases_per100kpop_7davg),
fillOpacity = 0.8,
weight = 1,
highlight = highlightOptions(
weight = 5,
color = "red",
bringToFront = TRUE),
label = lbls,
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto")) %>%
addLegend(pal = pal,
values = county_data_byday_now$new_cases_per100kpop_7davg,
title = "New Cases/100k (7d Avg)",
opacity = 0.9,
position = "bottomright") %>%
leafpop::addPopupGraphs(county_data_byday_now$pl_new,
group = "name",
width = 400, height = 200)
Age strata populations are taken from the US Census estimate for TN, year 2019 estimate (the most recent available)
# 2020 vs 2021 ####
# pl_cases_age_2020 <-
# tndoh_age %>%
# filter(date >= lubridate::ymd("2020-06-01"),
# date <= {lubridate::ymd(Sys.Date()) - 365}) %>%
# # date <= lubridate::ymd("2020-08-01")) %>%
# # filter(date >= lubridate::ymd("2021-06-01")) %>%
# # filter(date >= lubridate::ymd("2021-05-01"),
# # age.f %in% c("31-40", "41-50","51-60", "61-70", "71-80", "81+")) %>%
# ggplot(aes(x = date,
# y = ar_new_cases_100k_7davg,
# fill = age.f)) +
# geom_area(alpha = 0.5) +
# theme(legend.position = "none") +
# labs(x = NULL, y = NULL,
# title = "One year ago: Average Daily Cases by Age Group (per 100k)") +
# scale_x_date(date_breaks = "1 month", date_labels = "%b", expand = c(0, 0)) +
# facet_wrap(~age.f, nrow = 1) +
# # facet_wrap(~age.f, nrow = 1, scales = "free_y") +
# theme(axis.text.x.bottom = element_text(angle = 90, size = 7))
# pl_cases_age_peak <-
# tndoh_age %>%
# filter(date >= lubridate::ymd("2020-10-01"),
# date <= lubridate::ymd("2021-02-01")) %>%
# # date <= lubridate::ymd("2020-08-01")) %>%
# # filter(date >= lubridate::ymd("2021-06-01")) %>%
# # filter(date >= lubridate::ymd("2021-05-01"),
# # age.f %in% c("31-40", "41-50","51-60", "61-70", "71-80", "81+")) %>%
# ggplot(aes(x = date,
# y = ar_new_cases_100k_7davg,
# fill = age.f)) +
# geom_area(alpha = 0.5) +
# theme(legend.position = "none") +
# labs(x = NULL, y = NULL,
# title = "2020-2021 Winter peak: Average Daily Cases by Age Group (per 100k)") +
# scale_x_date(date_breaks = "2 months", date_labels = "%b-%Y", expand = c(0, 0)) +
# facet_wrap(~age.f, nrow = 1) +
# # facet_wrap(~age.f, nrow = 1, scales = "free_y") +
# theme(axis.text.x.bottom = element_text(angle = 90, size = 7))
# # axis.text.x = element_text(angle = 45, hjust = 1, size = 7))
pl_cases_age_2021 <-
tndoh_age %>%
filter(date >= lubridate::ymd("2021-06-01")) %>%
# filter(date >= lubridate::ymd("2021-05-01"),
# age.f %in% c("31-40", "41-50","51-60", "61-70", "71-80", "81+")) %>%
ggplot(aes(x = date,
y = ar_new_cases_100k_7davg,
fill = age.f)) +
geom_area(alpha = 0.5) +
theme(legend.position = "none") +
labs(x = NULL, y = NULL,
title = "Current: Average Daily Cases by Age Group (per 100k)") +
scale_x_date(date_breaks = "4 months", date_labels = "%b-%Y", expand = c(0, 0)) +
facet_wrap(~age.f, nrow = 1) +
# facet_wrap(~age.f, nrow = 1, scales = "free_y") +
theme(axis.text.x.bottom = element_text(angle = 90, size = 7))
# axis.text.x = element_text(angle = 45, hjust = 1, size = 7))
max_y <-
max(layer_scales(pl_cases_age_2021)$y$range$range[[2]])
# layer_scales(pl_cases_age_peak)$y$range$range[[2]])
tndoh_schoolage <-
tndoh_age %>%
select(date, age.f, ar_newcases) %>%
filter(age.f %in% c("0-10", "11-20")) %>%
group_by(date) %>%
summarise(ar_newcases_20 = sum(ar_newcases, na.rm = T))
# CASE TREND- STL DECOMPOSITION
tndoh_schoolage$ar_newcases_20_ts <- ts(tndoh_schoolage$ar_newcases_20, frequency = 7, start = 1)
tndoh_schoolage$new_cases_trend.stl <- stl(tndoh_schoolage$ar_newcases_20_ts, s.window = 7)$time.series[, "trend"]
pl_newcases_schoolage <-
tndoh_schoolage %>%
ggplot(aes(x=date)) +
geom_area(aes(y= ar_newcases_20), fill = "steelblue", alpha = 0.4) +
geom_line(aes(y=new_cases_trend.stl), color="red", size=1) +
labs(x = NULL, y = NULL,
title = "Daily New Cases in age groups 0-20") +
scale_x_date(date_breaks = "2 month", date_labels = "%b-%Y", expand = c(0, 0)) +
scale_y_continuous(expand = c(0, 0)) +
theme(axis.text.x.bottom = element_text(angle = 90, size = 7))
cowplot::plot_grid(pl_cases_age_2021 + ylim(c(0, max_y)) +
geom_text(data = age_death_table, size = 3, hjust = 0, vjust=1, fontface = "italic",
aes(x=lubridate::ymd("2021-06-01"), y = max_y,
label = sprintf("Vacc.:%1.1f%%", ar_recipientvax_pct) )),
# pl_cases_age_peak + ylim(c(0, max_y)),
pl_newcases_schoolage,
nrow = 2)
n_ul = 200
# PROPORTION BY AGE ####
# Unadjusted for population size
pl_age_proportion <-
tndoh_age %>%
ggplot(aes(date, ar_proportion_7davg,
group = age_range,
fill = age_range)) +
geom_area(color = "grey50", alpha = 0.6) +
scale_y_continuous(expand = c(0,0), limits = c(0,NA),
sec.axis = dup_axis(
breaks = tndoh_age_last$sec_axis.ticks_cases_mid,
labels = tndoh_age_last$sec_axis.label_cases,
name = NULL
)) +
scale_x_date(expand = c(0,0), date_breaks = "3 month", date_labels = "%b-%Y" ) +
theme(legend.position = "none",
axis.text.y.right = ggtext::element_markdown(size = 12)) +
labs(x = "Date",
y = "Proportion of cases by Age (Avg.)",
fill = "Age Group")
p_case_trend <-
tndoh_daily %>%
ggplot(aes(date, 1)) +
geom_tile(aes(fill = new_cases_7davg_100k )) +
viridis::scale_fill_viridis(limits = c(0, n_ul), oob = scales::squish,
option = "plasma") +
scale_x_date(expand = c(0,0), date_breaks = "3 month") +
scale_y_continuous(expand = c(0,0),
breaks = c(1),
labels = "New Cases\nAvg/d/100k") +
guides(fill = guide_colorbar(label.position = "bottom",
barwidth = 3,
label.theme = element_text(angle = 90, size = 8,
hjust = 0.5, vjust = 0.5))) +
labs(x=NULL, y = NULL, fill = NULL) +
theme(legend.position = "right",
legend.direction = "horizontal",
# legend.
axis.ticks.x = element_blank(),
axis.text.x = element_blank())
p2_legend <- cowplot::get_legend(p_case_trend)
p_case_trend_nolegend <- p_case_trend + theme(legend.position = "none")
pl_age_reference <-
tn_age_groups %>%
mutate(populationtn_byage_2019_prop = populationtn_byage_2019/sum(populationtn_byage_2019)) %>%
arrange(desc(age.f)) %>%
mutate( age_prop_mid = (populationtn_byage_2019_prop + lag(populationtn_byage_2019_prop, default = 0))/2,
age_prop_mid_cumsum = cumsum(age_prop_mid)) %>%
ggplot(aes(x=1, y = populationtn_byage_2019_prop,
fill = age.f)) +
geom_bar(position = "stack", stat = "identity", alpha = 0.6) +
geom_text(aes(x=1, y = age_prop_mid_cumsum,
label = sprintf("%s: %2.1f%%", age.f, populationtn_byage_2019_prop*100) ),
size = 3) +
scale_x_continuous(expand = c(0,0)) +
scale_y_continuous(expand = c(0,0)) +
labs(x="Population%\nby Age", y=NULL, fill = NULL) +
theme(legend.position = "none",
axis.ticks = element_blank(),
axis.text = element_blank(),
axis.line = element_blank())
p_case_trend_nolegend + patchwork::inset_element(p2_legend, 0.90, 0, .95, .90,
align_to = 'full',
clip = F) + patchwork::plot_spacer() +
pl_age_proportion + pl_age_reference +
patchwork::plot_layout(ncol = 2, heights = c(0.075,1), widths = c(1,.125))
Counties are ordered by highest new cases per day (current) based on the most recent 7 day average (per 100k population). Max limit of the color scale set at 200 to increase contrast- some local outbreaks exceed this limit (eg prison outbreaks in Bledsoe, Lake, Wayne Co.)
county_order <-
tndoh_countynew %>%
filter(date == max(date)) %>%
mutate(county_bynewcases = reorder(county, new_cases_per100kpop_7davg)) %>%
select(county, county_bynewcases)
# no longer interactive due to file size issues.
TN_heatmap_100k <-
tndoh_countynew %>%
filter(!county %in% c("Pending", "Out of State")) %>%
left_join(county_order, by = "county") %>%
mutate(new_cases = replace_na(new_cases, 0)) %>%
ggplot(aes(x = date,
y = county_bynewcases,
fill = new_cases_per100kpop_7davg)) +
geom_tile(aes(fill = new_cases_per100kpop_7davg)) +
viridis::scale_fill_viridis(labels = scales::comma,
limits = c(0, n_ul), oob = scales::squish,
option = "plasma") +
viridis::scale_color_viridis(labels = scales::comma,
limits = c(0, n_ul), oob = scales::squish,
option = "plasma") +
labs(fill = "New Cases/100k\n(7dAvg.)",
color = "New Cases/100k\n(7dAvg.)",
x = NULL, y = NULL) +
scale_x_date(date_breaks = "2 month", date_labels = "%b-%Y", expand = c(0, 0)) +
theme_minimal() +
theme(axis.text.y = element_text(size = 8),
panel.grid = element_blank(),
axis.text.x.bottom = element_text(angle = 90),
legend.position = "bottom",
legend.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size = 14),
legend.text = element_text(size = 10))
# x_TN_heatmap_100k <- girafe(ggobj = TN_heatmap_100k,
# options = list(
# opts_hover(css = "stroke-width:0.5")))
TN_heatmap_100k
Notes:
Note: assignment of new cases to county usually lags the total count. The State Health Dept. assigns cases to county of residence, rather than testing site. Individual county health departmens may report numbers differently, such as the number of cases diagnosed within the county, regardless of residence
# fig 2 State ####
pl_tnmap_tnsdc <-
tndoh_countynew %>%
filter(!is.na(date)) %>%
filter(date == max(date)) %>%
mutate(total_cases = case_when(total_cases == 0 ~ NA_real_,
TRUE ~ total_cases)) %>%
ggplot() +
# pl_tnmap_tnsdc <-
# pl_tnmap_tnsdc +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = as.numeric(total_cases)),
inherit.aes = F) +
geom_sf_text(aes(label= comma(total_cases, digits = 0),
geometry = geometry,
color = cut(total_cases, 2)),
size = 2.5) +
viridis::scale_fill_viridis(labels = scales::comma,
# limits = c(0, 100), oob = scales::squish,
option = "plasma") +
# viridis::scale_fill_viridis(option = "D", direction = -1) +
scale_color_manual(values = c( "white", "black"), guide = F) +
labs(fill = "Total Cases:")+
# caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via TN State Data Center") +
theme_void() +
theme(legend.position = "left")
# ggsave(here::here("figures/pl_tnmap.jpg"), dpi = 600, width = 12, height = 4, units = "in")
pl_tnmap_7dDelta <-
tndoh_countynew %>% #select(new_cases_county_7davg_1wkdelta)
filter(!is.na(date)) %>%
filter(date == max(date)) %>%
ggplot() +
# pl_tnmap_7dDelta <-
# pl_tnmap_7dDelta +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = new_cases_county_7davg_1wkdelta*100),
inherit.aes = F) +
scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
geom_sf_text(aes(label=sprintf("%1.1f",new_cases_county_7davg_1wkdelta*100), geometry = geometry),
# color = cut(new_cases_county_7davg_1wkdelta, 2)),
size = 3) +
# viridis::scale_fill_viridis(option = "D", direction = -1) +
# scale_color_manual(values = c( "black", "white"), guide = F) +
labs(fill = "1-Week\nChange %:")+
# caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via TN State Data Center") +
theme_void() +
theme(legend.position = "left")
# pl_tnmap_newhosp <-
# tndoh_countynew %>%
# filter(!is.na(date)) %>%
# filter(date == max(date)) %>%
# ggplot() +
# # pl_tnmap_newhosp <-
# # pl_tnmap_newhosp +
# geom_sf(data = TN_counties) +
# geom_sf(aes(geometry = geometry,
# fill = as.numeric(new_hospitalized)),
# inherit.aes = F) +
# geom_sf_text(aes(label=new_hospitalized, geometry = geometry),
# color = "black",
# # color = cut(new_hospitalized, 2)),
# size = 3) +
# scale_fill_gradient(low = "white", high = "red",
# limits = c(0, max(tndoh_countynew$new_hospitalized [tndoh_countynew$date == max(tndoh_countynew$date, na.rm = T)], na.rm = T))) +
# scale_color_manual(values = c( "white", "black"), guide = F) +
# labs(fill = "New\nHospitalizations:") +
# theme_void() +
# theme(legend.position = "left")
pl_tnmap_new_tnsdc <-
tndoh_countynew %>%
# filter(date == max(date)) %>%
filter(date == max(date)-1) %>%
ggplot() +
# pl_tnmap_new_tnsdc <-
# pl_tnmap_new_tnsdc +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = as.numeric(new_cases)),
inherit.aes = F) +
geom_sf_text(aes(label=new_cases, geometry = geometry),
color = "black",
# color = cut(new_cases, 2)),
size = 3) +
scale_fill_gradient(low = "white", high = "red",
limits = c(0, max(tndoh_countynew$new_cases[tndoh_countynew$date == max(tndoh_countynew$date, na.rm = T)], na.rm = T))) +
labs(fill = "New Cases:") +
# caption = "Visualization by @DrJMLuther\nData Source:@TNDeptofHealth via New York Times") +
theme_void() +
theme(legend.position = "left")
# ggsave(here::here("figures/pl_tnmap_new.jpg"), dpi = 600, width = 12, height = 4, units = "in")
pl_cases_pct <-
tndoh_countynew %>%
# filter(date == max(date, na.rm = T)-1) %>%
filter(date == max(date, na.rm = T)) %>%
ggplot() +
# pl_cases <-
# pl_cases +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = total_cases_per100kpop),
color = "grey30",
inherit.aes = F) +
geom_sf_text(aes(label=sprintf("%0.1f%%", total_cases_per100kpop/1000),
# geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos),
geometry = geometry,
color = cut(total_cases_per100kpop, 2)),
size = 3) +
viridis::scale_fill_viridis(option = "C", direction = 1) +
scale_color_manual(values = c("white", "black"), guide = F) +
# scale_color_manual(values = c( "black", "white"), guide = F) +
labs(fill = "Total Affected\n(Cases %Pop.):") +
theme_void() +
theme(legend.position = "left")
pl_FINAL_tnsdc <-
cowplot::plot_grid(title_tnsdc,
pl_tnmap_new_tnsdc,
pl_tnmap_7dDelta,
# pl_tnmap_newhosp,
# pl_tnmap_newtests_tnsdc,
pl_tnmap_tnsdc,
pl_cases_pct,
nrow = 5,
rel_heights = c(.1,1,1, 1, 1),
align = "y")
unclassified <-
tndoh_countynew %>%
filter(date == max(date),
county_lower %in% c("pending", "out of state"))
# pl_FINAL_tnsdc
pl_FINAL_tnsdc +
annotation_custom(gridExtra::tableGrob(
tibble(
Category = unclassified$county,
Cases = unclassified$new_cases),
rows=NULL, cols = NULL,
theme = gridExtra::ttheme_minimal(base_size = 10)),
xmin=.85, xmax=.95, ymin=0.66, ymax=.9) +
# annotation_custom(gridExtra::tableGrob(
# tibble(
# Category = unclassified$county,
# Cases = unclassified$new_hospitalized),
# rows=NULL, cols = NULL,
# theme = gridExtra::ttheme_minimal(base_size = 10)),
# xmin=.85, xmax=.95, ymin=0.31, ymax=.275) +
annotation_custom(gridExtra::tableGrob(
tibble(
Category = unclassified$county,
Cases = unclassified$total_cases),
rows=NULL, cols = NULL,
theme = gridExtra::ttheme_minimal(base_size = 10)),
xmin=.85, xmax=.95, ymin=0.26, ymax=.375)
The estimate of case fatality is presented as the crude deaths/cases ratio, and is not adjusted for age. This could markedly affect the rate in different counties. This is also affected by the ability to test and identify as many cases as possible, which may differ by county.
pl_deaths_total <-
tndoh_countynew %>%
# filter(date == max(date, na.rm = T)) %>%
filter(date == max(date, na.rm = T)-1) %>%
ggplot() +
# pl_deaths <-
# pl_deaths +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = total_deaths),
color = "grey30",
inherit.aes = F) +
geom_sf_text(aes(label=sprintf("%0.0f", total_deaths),
# geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths, pct_pos),
geometry = geometry,
color = cut(total_deaths, 2)),
size = 3) +
viridis::scale_fill_viridis(option = "C", direction = 1) +
scale_color_manual(values = c("white", "black"), guide = F) +
# scale_color_manual(values = c( "black", "white"), guide = F) +
labs(fill = "Total Deaths:") +
theme_void() +
theme(legend.position = "left")
pl_deaths_percapita <-
tndoh_countynew %>%
# filter(date == max(date, na.rm = T)) %>%
filter(date == max(date, na.rm = T)-1) %>%
ggplot() +
# pl_deaths <-
# pl_deaths +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = total_deaths_per100kpop),
color = "grey30",
inherit.aes = F) +
geom_sf_text(aes(label=sprintf("%0.0f", total_deaths_per100kpop),
# geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_deaths_per100kpop, pct_pos),
geometry = geometry,
color = cut(total_deaths_per100kpop, 2)),
size = 3) +
viridis::scale_fill_viridis(option = "C", direction = 1) +
scale_color_manual(values = c("white", "black"), guide = F) +
# scale_color_manual(values = c( "black", "white"), guide = F) +
labs(fill = "Total Deaths\nper 100k Pop.:") +
theme_void() +
theme(legend.position = "left")
pl_cfr <-
tndoh_countynew %>%
# filter(date == max(date, na.rm = T)) %>%
filter(date == max(date, na.rm = T)-1) %>%
mutate(cfr = total_deaths/total_cases) %>%
ggplot() +
# pl_cfr <-
# pl_cfr +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = cfr*100),
color = "grey30",
inherit.aes = F) +
geom_sf_text(aes(label=sprintf("%0.1f %%", cfr*100),
# geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos),
geometry = geometry,
color = cut(cfr, 2)),
size = 3) +
viridis::scale_fill_viridis(option = "C", direction = 1) +
scale_color_manual(values = c("white", "black"), guide = F) +
# scale_color_manual(values = c( "black", "white"), guide = F) +
labs(fill = "Case Fatality %:") +
theme_void() +
theme(legend.position = "left")
cowplot::plot_grid(title_tnsdc,
pl_deaths_total,
pl_deaths_percapita,
pl_cfr,
nrow = 4, rel_heights = c(.1,1,1,1))
# ggsave(here::here("figures/pl_deaths_map.jpg"), dpi = 600, width = 12, height = 6, units = "in")
The estimated Tennessee county populations are taken from the US Census (2019 estimate). The estimate of case fatality is presented as the crude deaths/cases ratio, and is not adjusted for age. This could markedly affect the rate in different counties
pl_newcases_percapita <-
tndoh_countynew %>% #select(new_cases_per100kpop_7davg)
filter(date == max(date, na.rm = T)) %>%
# filter(date == max(date, na.rm = T)-1) %>%
ggplot() +
# pl_newcases <-
# pl_newcases +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = new_cases_per100kpop_7davg),
color = "grey30",
inherit.aes = F) +
geom_sf_text(aes(label=sprintf("%0.0f", new_cases_per100kpop_7davg),
# geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos),
geometry = geometry,
color = cut(new_cases_per100kpop_7davg, 2)),
size = 3) +
viridis::scale_fill_viridis(option = "C", direction = 1) +
scale_color_manual(values = c("white", "black"), guide = F) +
# scale_color_manual(values = c( "black", "white"), guide = F) +
labs(fill = "Daily Cases per\n100K Pop.\n(7d Avg):") +
theme_void() +
theme(legend.position = "left")
# pl_test_percapita <-
# tndoh_countynew %>%
# mutate(total_tests_pctpop = total_tests/population*100000) %>%
# # select(total_tests_per100kpop)
# filter(date == max(date, na.rm = T)) %>%
# # filter(date == max(date, na.rm = T)-1) %>%
# ggplot() +
# # pl_cfr <-
# # pl_cfr +
# geom_sf(data = TN_counties) +
# geom_sf(aes(geometry = geometry,
# fill = total_tests_pctpop),
# color = "grey30",
# inherit.aes = F) +
# geom_sf_text(aes(label=sprintf("%0.1f ", total_tests_pctpop/100),
# # geom_sf_text(aes(label=sprintf("%0.0f\n%0.0f%%", total_cases, pct_pos),
# geometry = geometry,
# color = cut(total_tests_pctpop, 2)),
# size = 3) +
# viridis::scale_fill_viridis(option = "C", direction = 1) +
# scale_color_manual(values = c("white", "black"), guide = F) +
# # scale_color_manual(values = c( "black", "white"), guide = F) +
# labs(fill = "Tests per\n100k Population:") +
# theme_void() +
# theme(legend.position = "left")
pl_test_positivity <-
tndoh_countynew %>%
filter(date == max(date, na.rm = T)) %>%
ggplot() +
geom_sf(data = TN_counties) +
geom_sf(aes(geometry = geometry,
fill = new_testpos_pct_7davg),
color = "grey30",
inherit.aes = F) +
geom_sf_text(aes(label=sprintf("%0.1f ", new_testpos_pct_7davg),
geometry = geometry,
color = cut(new_testpos_pct_7davg, 2)),
size = 3) +
viridis::scale_fill_viridis(option = "C", direction = 1) +
scale_color_manual(values = c("white", "black"), guide = F) +
labs(fill = "Test+%\n(7d avg):") +
theme_void() +
theme(legend.position = "left")
cowplot::plot_grid(title_tnsdc, pl_newcases_percapita, pl_cases_pct, pl_test_positivity,
nrow = 4, rel_heights = c(.1,1,1,1))
# ggsave(here::here("figures/pl_deaths_map.jpg"), dpi = 600, width = 12, height = 6, units = "in")
Case data for these maps are taken from the CDC. and the 2019 state population estimates are taken from the United States Census Bureau.
Formerly I used the COVID Tracking project, but this site is no longer updating data as of March 7, 2021.
CDC Data Currently unavailable
mouse over state for details
max_newcaserate <- max(cdc_currentdata_2163$new_cases_state_7davg_100k[cdc_currentdata_2163$date == max(cdc_currentdata_2163$date)], na.rm = T)
pl_us_cases_100k <-
cdc_currentdata_2163 %>%
filter(date == max(date)) %>%
ggplot(aes(geometry = geometry,
fill = new_cases_state_7davg_100k)) +
# pl_us_cases_100k <-
# pl_us_cases_100k +
ggiraph::geom_sf_interactive(aes(tooltip = tt,
data_id = state), size = .2) +
viridis::scale_fill_viridis(labels = scales::comma,
limits = c(0, max_newcaserate), oob = scales::squish,
option = "plasma") +
labs(fill = "Cases/100k/d\n(7dAvg)") +
theme_void()
ggiraph::girafe(code = print( pl_us_cases_100k),
options = list(
opts_hover(css = "stroke:orange;stroke-width:2;")))
max_newcaserate_alldates <- max(cdc_currentdata$new_cases_state_7davg_100k, na.rm = T)
state_order <-
cdc_currentdata_2163 %>%
filter(date == max(date)) %>%
mutate(state_bynewcases = reorder(state, new_cases_state_7davg_100k)) %>%
select(state, state_bynewcases)
us_heatmap_100k <-
cdc_currentdata %>%
left_join(state_order, by = "state") %>%
filter(date >= lubridate::ymd("2020-04-01")) %>%
ggplot(aes(x = date,
y = state_bynewcases)) +
geom_tile(aes(fill = new_cases_state_7davg_100k)) +
viridis::scale_fill_viridis(labels = scales::comma,
limits = c(0, max_newcaserate_alldates), oob = scales::squish,
option = "plasma") +
viridis::scale_color_viridis(labels = scales::comma,
limits = c(0, max_newcaserate_alldates), oob = scales::squish,
option = "plasma") +
labs(fill = "Cases/100k/d\n(7dAvg)",
color = "Cases/100k/d\n(7dAvg)",
x = NULL, y = NULL) +
scale_x_date(date_breaks = "1 month", date_labels = "%b", expand = c(0, 0)) +
theme_minimal() +
theme(axis.text.y = element_text(size = 12),
panel.grid = element_blank(),
legend.position = "bottom",
legend.background = element_blank(),
panel.border = element_blank(),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12))
us_heatmap_100k
You made it to the end, so here’s a random XKCD comic (please excuse any profanity- just randomly selected):
Hmisc::markupSpecs$html$session()
R version 4.3.1 (2023-06-16 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 11 x64 (build 22621) Matrix products: default attached base packages: [1] stats graphics grDevices utils datasets methods base other attached packages: [1] rvest_1.0.3 httr_1.4.7 Hmisc_5.1-1 [4] patchwork_1.1.3 ggtext_0.1.2 htmlTable_2.4.1 [7] RcppRoll_0.3.0 cowplot_1.1.1 urbnmapr_0.0.0.9002 [10] ggiraph_0.8.7 formattable_0.2.1 leafpop_0.1.0 [13] leaflet_2.2.0 sf_1.0-14 lubridate_1.9.3 [16] forcats_1.0.0 stringr_1.5.0 dplyr_1.1.3 [19] purrr_1.0.2 readr_2.1.4 tidyr_1.3.0 [22] tibble_3.2.1 ggplot2_3.4.4 tidyverse_2.0.0To cite R in publications use:
R Core Team (2023). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org/.
To cite the Hmisc package in publications use:Harrell Jr F (2023). Hmisc: Harrell Miscellaneous. R package version 5.1-1, https://CRAN.R-project.org/package=Hmisc.
To cite the ggplot2 package in publications use:Wickham H (2016). ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. ISBN 978-3-319-24277-4, https://ggplot2.tidyverse.org.