Italy/UK death slope comparison

covid_ItUK.R

## adapted from https://rviews.rstudio.com/2020/03/05/covid-19-epidemiology-with-r/

## CC-BY-SA 3.0 https://creativecommons.org/licenses/by-sa/3.0/

library(tidyverse)
library(lubridate)
library(directlabels)
library(colorspace)
library(emmeans)
library(broom)
library(robust)

jhu_url <- paste("https://raw.githubusercontent.com/CSSEGISandData/", 
  "COVID-19/master/csse_covid_19_data/", "csse_covid_19_time_series/", 
  "time_series_19-covid-Deaths.csv", sep = "")

dd <- read_csv(jhu_url)

dd2 <- (
    dd
    %>% rename(province = "Province/State"
             , country_region = "Country/Region") %>%
        pivot_longer(-c(province, country_region, Lat, Long),
                     names_to = "Date", values_to = "cumulative_deaths")
    %>% group_by(country_region)
    ## adjust JHU dates back one day to reflect US time, more or less
    %>%     mutate(Date = mdy(Date) - days(1),
                   daily_deaths=c(NA,diff(cumulative_deaths)))
    %>% filter(Date>=as.Date("2020-02-23"), ## match DS data
               country_region %in% c("Italy","United Kingdom"),
               ## omit Gibraltar/Channel Islands/etc.
               is.na(province) | province=="United Kingdom")
    %>% mutate(diffDate=as.numeric(Date-min(Date)))
    %>% ungroup()
)

theme_set(theme_classic())
mm <- lm(log(daily_deaths)~diffDate*country_region,
         data=dd2,
         subset=daily_deaths>0)
summary(mm)


## try robust lm
cc <- mm$call
cc[[1]] <- quote(lmRob)
mmr <- eval(cc)

summary(mmr)

aa <- (augment(mm)
        ##       newdata=expand.grid(diffDate=unique(dd2$diffDate),
        ##                           country_region=unique(dd2$country_region)))
    %>% mutate(lwr=exp(.fitted-qnorm(0.975)*.se.fit),
               upr=exp(.fitted+qnorm(0.975)*.se.fit),
               daily_deaths=exp(.fitted),
               Date=diffDate+min(dd2$Date))
)

gg0 <- (ggplot(dd2,aes(Date,daily_deaths,colour=country_region))
    + scale_y_log10() 
    + geom_point()
    ## + geom_smooth(method="lm")
    + geom_line(data=aa)
    + geom_ribbon(data=aa,aes(ymin=lwr,ymax=upr,fill=country_region),
                  alpha=0.2,colour=NA)
    + scale_colour_discrete_qualitative()
    + geom_dl(aes(label=country_region),method="smart.grid")
    + theme(legend.position="none")
)

ggsave("ituk.png")
png("ituk2.png",width=6,height=3,units="in",res=150)
plot(emtrends(mm,var="diffDate","country_region",transform="log"))
dev.off()