In rplzzz/CovMitigation: COVID-19 infection model for Virginia

suppressPackageStartupMessages(library(dplyr))
library(tidyr)
library(ggplot2)
library(CovMitigation)
library(vacovdata)
library(ggthemes)
library(latex2exp)

Assumptions

We have to make a slight approximation here. The growth rate we are interested in is the growth in active cases, which is affected by both the infection rate and the recovery rate. The data from the department of health only has new cases and cumulative cases. For this analysis we went with cumulative cases because the new case data is rather noisy.

Growth observations

#strtdate <- as.Date('2020-03-23')
strtdate <- as.Date('2020-04-05')
ftest_normdate <- as.Date('2020-04-15')
pltdata_new <- 
  select(vacovdata, date, vaNewCases=positiveIncrease, ftest, ntest=totalTestResultsIncrease) %>%
  filter(ntest > 0) %>% 
  arrange(date)

pltdata_old <- 
  select(uvads_covid19, date, vaNewCases, ftest, ntest) %>%
  filter(ntest > 0) %>% 
  arrange(date)

pltdata <- pltdata_new

## Compute the testing fraction correction.
normidx <- which(ftest_normdate == pltdata$date)
pltdata$ftest_relative <- pltdata$ftest / pltdata$ftest[normidx]
#fcorrect <- min(pltdata$ftest) / pltdata$ftest

pltdata$correctedNewCases <- pltdata$vaNewCases / pltdata$ftest_relative 
pltdata$correctedCumCases <- cumsum(pltdata$correctedNewCases)
pltdata$rawCumCases <- cumsum(pltdata$vaNewCases)

pltdata <- filter(pltdata, date >= strtdate)
relidx <- 1
pltdata$correctedNewCasesRel <- pltdata$correctedNewCases / pltdata$correctedNewCases[relidx]
pltdata$correctedCumCasesRel <- pltdata$correctedCumCases / pltdata$correctedCumCases[relidx]
pltdata$rawCumCasesRel <- pltdata$rawCumCases / pltdata$rawCumCases[relidx]

t <- seq_along(pltdata$date) - 1  # start at zero
pltdata$t <- t

## Label positioning
lastidx <- length(t)
lastdate <- pltdata$date[lastidx]

## Create the benchmark data
tds <- c(3, 5, 7, 9, 11, 13)
tdlabs <- paste('td=',tds)
tdlabstex <- sprintf('$t_d = %d$', tds)
gr <- lapply(tds, 
             function(td) {
               r <- log(2)/td
               exp(r*t) / exp(r*t[relidx])
             })
benchmarks <- as.data.frame(gr)
names(benchmarks) <- tdlabs
benchmarks$date <- t + pltdata$date[1]
benchmarks <- pivot_longer(benchmarks, -date, names_to = 'doubling time')
benchplt <- 
  ggplot(data=pltdata, aes(x=date)) + 
  geom_line(data=pltdata, mapping=aes(y=correctedCumCasesRel), size=1.2) +
  geom_line(data=benchmarks, mapping=aes(y=value, group=`doubling time`), color='lightgrey') +
  scale_color_solarized() +
  scale_y_log10() +
  ylab('Corrected cumulative cases (relative)') +
  theme_bw()
for(i in seq_along(gr)) {
  benchplt <- benchplt + geom_label(x=lastdate, y=log10(gr[[i]][lastidx]), label=TeX(tdlabstex[i]))
}
suppressWarnings(print(benchplt))

The rate of testing has been growing rapidly. If you leave out the correction for the number of tests performed, you get this:

benchplt2 <- 
  ggplot(data=pltdata, aes(x=date)) + 
  geom_line(data=pltdata, mapping=aes(y=rawCumCasesRel), size=1.2) +
  geom_line(data=benchmarks, mapping=aes(y=value, group=`doubling time`), color='lightgrey') +
  scale_color_solarized() +
  scale_y_log10() +
  ylab('Cumulative cases (relative)') +
  theme_bw()
for(i in seq_along(gr)) {
  benchplt2 <- benchplt2 + geom_label(x=lastdate, y=log10(gr[[i]][lastidx]), label=TeX(tdlabstex[i]))
}
suppressWarnings(print(benchplt2))

Finally, here is the growth for Charlottesville and Albemarle County. Note, however, that we are getting the testing rate correction from statewide data, since we don't have that at the county level, so there might be some bias in this one.

#strtdate <- as.Date('2020-03-28')


chopltdata <-
  filter(NYTimesCOVID19::cov19county, fips==51540 | fips==51003) %>%
  group_by(date) %>%
  summarise(cases = sum(cases))
chopltdata$newCases <- c(0, diff(chopltdata$cases))
chopltdata <-
  left_join(chopltdata, 
            select(vacovdata, date, vaNewCases=positiveIncrease, ftest, ntest=totalTestResultsIncrease), 
            by='date') %>%
  select(date, newCases, ftest) %>%
  filter(ftest > 0, newCases >0)

## Compute the testing fraction correction.
normidx <- 1
chopltdata$ftest_relative <- chopltdata$ftest / chopltdata$ftest[normidx]
#fcorrect <- min(chopltdata$ftest) / chopltdata$ftest

relidx <- 1
chopltdata$correctedNewCases <- chopltdata$newCases / chopltdata$ftest_relative 
chopltdata$correctedCumCases <- cumsum(chopltdata$correctedNewCases)
chopltdata$rawCumCases <- cumsum(chopltdata$newCases)

chopltdata <- filter(chopltdata, date >= strtdate)
chopltdata$correctedNewCasesRel <- chopltdata$correctedNewCases / chopltdata$correctedNewCases[relidx]
chopltdata$correctedCumCasesRel <- chopltdata$correctedCumCases / chopltdata$correctedCumCases[relidx]
chopltdata$rawCumCasesRel <- chopltdata$rawCumCases / chopltdata$rawCumCases[relidx]

t <- seq_along(chopltdata$date) - 1  # start at zero
chopltdata$t <- t

## Label positioning
lastidx <- length(t)
lastdate <- chopltdata$date[lastidx]

## Create the benchmark data
tds <- c(7,9,11,13,15)
tdlabs <- paste('td=',tds)
tdlabstex <- sprintf('$t_d = %d$', tds)
gr <- lapply(tds, 
             function(td) {
               r <- log(2)/td
               exp(r*t) / exp(r*t[relidx])
             })
benchmarks <- as.data.frame(gr)
names(benchmarks) <- tdlabs
benchmarks$date <- t + chopltdata$date[1]
benchmarks <- pivot_longer(benchmarks, -date, names_to = 'doubling time')
benchplt <- 
  ggplot(data=chopltdata, aes(x=date)) + 
  geom_line(data=chopltdata, mapping=aes(y=correctedCumCasesRel), size=1.2) +
  geom_line(data=benchmarks, mapping=aes(y=value, group=`doubling time`), color='lightgrey') +
  scale_color_solarized() +
  scale_y_log10() +
  ylab('Corrected cumulative cases (relative)') +
  theme_bw()
for(i in seq_along(gr)) {
  benchplt <- benchplt + geom_label(x=lastdate, y=log10(gr[[i]][lastidx]), label=TeX(tdlabstex[i]))
}
suppressWarnings(print(benchplt))

benchpltcho2 <- 
  ggplot(data=chopltdata, aes(x=date)) + 
  geom_line(data=chopltdata, mapping=aes(y=rawCumCasesRel), size=1.2) +
  geom_line(data=benchmarks, mapping=aes(y=value, group=`doubling time`), color='lightgrey') +
  scale_color_solarized() +
  scale_y_log10() +
  ylab('Cumulative cases (relative)') +
  theme_bw()
for(i in seq_along(gr)) {
  benchpltcho2 <- benchpltcho2 + geom_label(x=lastdate, y=log10(gr[[i]][lastidx]), label=TeX(tdlabstex[i]))
}
suppressWarnings(print(benchpltcho2))