COVID Tracking Project Enhancements to Johns Hopkins Case/Fatality Data

.libPaths()

options(warn=-1) 
options(scipen = 20)
options(datatable.print.topn=100)
options(datatable.showProgress=FALSE)
usualsuspects <- c(
    'tidyverse', 'data.table', 'pryr', 'plyr','readxl', 'bit', 'grid', 'gridExtra',     
    'rvest', 'magrittr','lubridate','rlist', 'tictoc', 'skimr','pracma','TTR',
    'fst','feather','rio', 'Rcpp','tidyquant', 'Matrix', 'arrow',
    'knitr', 'kableExtra', 'microbenchmark', "Metrics",
    'ggplot2','RColorBrewer',
    'magick',"webshot","IRdisplay",    
    'incidence','outbreaks','EpiEstim', 'modeltime', 
    'forecast', 'timetk', 'h2o','xgboost', 'catboost', 'lightgbm'
) 
suppressMessages(invisible(lapply(usualsuspects, library, character.only = TRUE)))
funcsdir <- "/steve/r/functions"
funcsfile <- "rfunctions.r"
setwd(funcsdir)
source(funcsfile)
cppdir <- "c:/steve/c++/rstudio"
setwd(cppdir)
blanks(1)
lsf.str()
set.seed(543)
blanks(2)

diffna <- function (var) c(NA,diff(var))
prodna <- function (var) prod(var,na.rm=TRUE)
sumna <- function (var) sum(var,na.rm=TRUE)
meanna <- function (var) round(mean(var,na.rm=TRUE))
freadfac <- function (var) data.table::fread(var,stringsAsFactors=TRUE)
fstrddt <- function(var) fst::read_fst(var,as.data.table=TRUE)
movavg <- function (var, lags) round(data.table::frollmean(var,lags))
meancov <- function(var) round(mean(var))
wkday <- function(dt) lubridate::wday(dt,label=TRUE)
blanks(2)

url <- "https://covidtracking.com/data/national/hospitalization"
css <- "td"
html <- (read_html(url) %>% html_nodes(css) %>% html_text())
head(html,20)
blanks(2)

html <- html %>% gsub("Currently on ventilator/Now on ventilator","",.) %>% gsub("Currently in ICU/Now in ICU","",.) %>%  
            gsub("Currently hospitalized/Now hospitalized","",.)  %>% gsub("Date","",.)
dt <- data.table(matrix(html,ncol=4,byrow=TRUE)) %>% setnames(.,c("date","hospitalized","ICU",'ventilator'))
meta(dt)
blanks(1)
mykab(head(dt))
blanks(2)

invisible(dt[,c("date","hospitalized","ICU","ventilator"):=.(mdy(date),gsub(",","",hospitalized),gsub(",","",ICU),gsub(",","",ventilator))][,
            hospitalized:=ifelse(hospitalized=="N/A",NA,as.integer(hospitalized))][,ICU:=ifelse(ICU=="N/A",NA,as.integer(ICU))][,
            ventilator:=ifelse(ventilator=="N/A",NA,as.integer(ventilator))])
dt <- dt[order(date)]
setkey(dt,date)
meta(dt)
blanks(1)
prheadtail(dt,howmany=10)
blanks(2)

wdir <- "/data/covidtracking"
setwd(wdir)
fname <- "hospitalizations.fst"
write_fst(dt,fname)
blanks(2)

prheadtail(dt[date>"2020-03-25"])
blanks(2)

wd <- "c:/steve/covid/magrsave"
setwd(wd)
fname <- tail(list.files( pattern = "mfinal*"),1)
print(fname)
mfinalus <- import(fname,setclass="data.table")
meta(mfinalus)
print(max(mfinalus$date))
blanks(2)

mykab(mfinalus[,.(daycases=sumna(daycases),daydeaths=sumna(daydeaths)),.(date)][
    ,.(.N,daycases=sumna(daycases),daydeaths=sumna(daydeaths)),.(wkday(date))][
    ,`:=`(pctcases=round(100*daycases/sum(daycases),2),pctdeaths=round(100*daydeaths/sum(daydeaths),2))])
blanks(2)

mykab(mfinalus[,.(daycases=sumna(daycases),daydeaths=sumna(daydeaths)),.(aftermay=date>"2020-05-31")][
    ,.(aftermay,fatalrate=round(100*(daydeaths/daycases),2))])
blanks(2)

covidjhu <- mfinalus[,.(daycases=sumna(daycases),daydeaths=sumna(daydeaths)),.(date)]
covidfinal <- dt[covidjhu][,dayofweek:=wkday(date)][,aftermay:=date>"2020-05-31"][,idx:=.N:1]
meta(covidfinal)
blanks(1)
prheadtail(covidfinal) 
blanks(2)

wdir <- "/data/covidtracking"
setwd(wdir)

fname <- "coviddaily.fst"
write_fst(covidfinal,fname)
fname <- "coviddaily.feather"
write_feather(covidfinal,fname)
fname <- "coviddaily.parquet"
write_parquet(covidfinal,fname)
blanks(2)

Build the data.table to be used for forecasting.

In [15]:

ncovidfinal <- covidfinal[date>"2020-03-25"][,c("daycases_lag21","hospitalized_lag14","ICU_lag7")
                          :=.(shift(daycases, 21L, fill = NA, type = "lag"),
                             shift(hospitalized, 14L, fill = NA, type = "lag"),
                             shift(ICU, 7L, fill = NA, type = "lag"))]
nr <- nrow(ncovidfinal)
meta(ncovidfinal)
blanks(2)

ggplot(ncovidfinal,aes(x=daycases,y=ICU, color=aftermay)) + geom_point() + geom_smooth(method='loess', se=FALSE)

ggplot(ncovidfinal[22:nr],aes(x=daycases_lag21,y=daydeaths,color=aftermay)) + geom_point() + geom_smooth(method='loess', se=FALSE)

ggplot(ncovidfinal[15:nr],aes(x=hospitalized_lag14,y=daydeaths,color=aftermay)) + geom_point() + geom_smooth(method='loess', se=FALSE)

ggplot(ncovidfinal,aes(x=ICU_lag7,y=daydeaths,color=aftermay)) + geom_point() + geom_smooth(method='loess', se=FALSE)

depvar <- c("daydeaths")
indvars <- c("aftermay","dayofweek","daycases_lag21","hospitalized_lag14","ICU_lag7")
vars <- c(depvar,indvars)
blanks(2)

lst <- lm(as.formula(paste(depvar,"~",paste(indvars,collapse="+"),sep="")),data=ncovidfinal[22:nr,..vars])
summary(lst)
blanks(2)

v <- var(ncovidfinal[22:nr]$daydeaths)
m <- mse(ncovidfinal[22:nr]$daydeaths,predict(lst))
r2 <- 1-m/v
v
m
r2
blanks(2)

trainspm <- xgb.DMatrix(sparse.model.matrix(as.formula(paste0(depvar,"~.")), data = ncovidfinal[22:nr,..vars]), label =  ncovidfinal[22:nr][[depvar]])
str(trainspm)
bst <- xgb.train(data = trainspm, max_depth = 6, eta = .3, nthread = 2, nrounds = 10, subsample = .3, tree_method="exact",
              gamma=5,verbose=1, subsample=.5,lambda=2,alpha=2,objective = "reg:squarederror",eval_metric = "rmse")
blanks(2)

v <- var(ncovidfinal[22:nr]$daydeaths)
m <- mse(ncovidfinal[22:nr]$daydeaths,predict(bst,newdata=trainspm))
r2 <- 1-m/v
v
m
r2
blanks(2)

importance <- xgb.importance(feature_names = colnames(trainspm), model = bst)
xgb.plot.importance(importance, rel_to_first = TRUE, xlab = "Relative Importance")
blanks(2)

newdt <- data.table(date=ncovidfinal[idx==1]$date+1:7,
                    aftermay=rep(TRUE,7),
                    dayofweek=ncovidfinal[between(idx,1,7),dayofweek],
                    daycases_lag21=ncovidfinal[between(idx,15,21),daycases],
                    hospitalized_lag14=ncovidfinal[between(idx,8,14),hospitalized], 
                    ICU_lag7=ncovidfinal[between(idx,1,7),ICU]                                       
                   )
meta(newdt)
blanks(1)
mykab(newdt)
blanks(2)

newspm <- xgb.DMatrix(sparse.model.matrix(as.formula(paste0("~.")), data = newdt[,..indvars]))
dimnames(newspm)
blanks(2)

round(predict(lst,newdata=newdt))
sum(round(predict(lst,newdata=newdt)))
blanks(1)
round(predict(bst,newdata=newspm))
sum(round(predict(bst,newdata=newspm)))
blanks(1)
round((predict(lst,newdata=newdt)+predict(bst,newdata=newspm))/2)
sum(round((predict(lst,newdata=newdt)+predict(bst,newdata=newspm))/2))
blanks(2)

1: 1753 2: 2390 3: 2504 4: 2405 5: 2416 6: 2180 7: 1702

vars = c("aftermay","dayofweek","daycases_lag21","hospitalized_lag14","ICU_lag7","daydeaths")
vars <- c("date",vars)
blanks(2)

newdtt <- rbind(tail(ncovidfinal[,..vars],7),
    cbind(newdt,predicted=round((predict(lst,newdata=newdt)+predict(bst,newdata=newspm))/2)) 
              ,use.names=F)
mykab(newdtt)
blanks(2)

wd <- "c:/steve/covid/modeling"
setwd(wd)
ofile <- paste0("newdtt",newdtt[7]$date,".png")
setnames(newdtt,c("date","aftermay","weekday","cases_l21","hosp_l14","ICU_l7","deaths"))
kbl(newdtt[,-c("aftermay")]) %>%
    kable_paper("striped", full_width = FALSE) %>%  	 
    column_spec(6, color = "white",
    background = spec_color(c(rep(1,7),rep(3,7)), begin = 0.4, end = 0.7)) %>%
    as_image(width=10,file=ofile)

setnames(newdtt,c("date","aftermay","weekday","cases_l21","hosp_l14","ICU_l7","deaths"))
kbl(newdtt[,-c("aftermay")]) %>%
    kable_paper("striped", full_width = FALSE) %>%  	 
    column_spec(6, color = "white",
    background = spec_color(c(rep(1,7),rep(3,7)), begin = 0.4, end = 0.7)) %>%
    as.character() %>%
    IRdisplay::display_html()