under_reporting.Rmd

---
title: "Data completeness state"
author: "Dan Weinberger"
date: "4/22/2020"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ExcessILI)
library(cdcfluview)
library(reshape2)
library(ggplot2)
library(lubridate)
library(RColorBrewer)
library(plotly)
library(MMWRweek)
library(readr)
library(rjson)
library(htmlTable)
library(RSocrata)
library(pdftools)
library(readr)
library(gsubfn)
library(INLA)
library (RCurl)
library(rjags)
library(HDInterval)
library(pbapply)
library(parallel)

```
## Backfilling
#NOTE THE DATA ARE MISSING FOR JULY-OCT 2019

```{r}
#Data from Andrew: NOTE: it looks like July-oct 2019 is missing from these data
wk11 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data11.csv')
wk11$max.date.report <- as.Date('2020-03-14')
wk11$report.date <- as.Date('2020-03-27')

wk12 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data12.csv')
wk12$max.date.report <- as.Date('2020-03-21')
wk12$report.date <- as.Date('2020-04-03')

wk13 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data13.csv')
wk13$max.date.report <- as.Date('2020-03-28')
wk13$report.date <- as.Date('2020-04-10')

wk15 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data15.csv')
wk15$max.date.report <- as.Date('2020-04-11')
wk15$report.date <- as.Date('2020-04-17')

wk16 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data16.csv')
wk16$max.date.report <- as.Date('2020-04-18')
wk16$report.date <- as.Date('2020-04-24')

wk17 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data17.csv')
wk17$max.date.report <- as.Date('2020-04-25')
wk17$report.date <- as.Date('2020-05-01')

wk18 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data18.csv')
wk18$max.date.report <- as.Date('2020-05-02')
wk18$report.date <- as.Date('2020-05-08')

wk19 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data19.csv')
wk19$max.date.report <- as.Date('2020-05-09')
wk19$report.date <- as.Date('2020-05-15')


wk20 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data20.csv')
wk20$max.date.report <- as.Date('2020-05-16')
wk20$report.date <- as.Date('2020-05-22')

wk21 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data21.csv')
wk21$max.date.report <- as.Date('2020-05-23')
wk21$report.date <- as.Date('2020-05-29')

wk22 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data22.csv')
wk22$max.date.report <- as.Date('2020-05-30')
wk22$report.date <- as.Date('2020-06-05')

wk23 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data23.csv')
wk23$max.date.report <- as.Date('2020-06-06')
wk23$report.date <- as.Date('2020-06-12')

wk24 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data24.csv')
wk24$max.date.report <- as.Date('2020-06-13')
wk24$report.date <- as.Date('2020-06-19')

wk25 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data25.csv')
wk25$max.date.report <- as.Date('2020-06-20')
wk25$report.date <- as.Date('2020-06-26')

wk26 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data26.csv')
wk26$max.date.report <- as.Date('2020-06-27')
wk26$report.date <- as.Date('2020-07-03')

wk27 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data27.csv')
wk27$max.date.report <- as.Date('2020-07-04')
wk27$report.date <- as.Date('2020-07-10')

wk28 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data28.csv')
wk28$max.date.report <- as.Date('2020-07-11')
wk28$report.date <- as.Date('2020-07-17')

wk29 <- 
  read.csv('./Archives_FluView_state_data/State_Custom_Data29.csv')
wk29$max.date.report <- as.Date('2020-07-18')
wk29$report.date <- as.Date('2020-07-24')


all.reports <- rbind.data.frame(wk11,wk12,wk13,wk15, wk16, wk17, wk18, wk19, wk20, wk21, wk22, wk23,wk24, wk25)


all.reports$epiyr <-
  as.numeric(as.character(substr(all.reports$SEASON,1,4)))

all.reports$year <- all.reports$epiyr
all.reports$year[all.reports$WEEK<=26] <-
  all.reports$epiyr[all.reports$WEEK<=26] +1
all.reports$week.death <-
  mmwr_week_to_date(all.reports$year, all.reports$WEEK)+6

#Fix formatting for the count variables
all.reports$NUM.INFLUENZA.DEATHS <-
  gsub(',','',all.reports$NUM.INFLUENZA.DEATHS) 
all.reports$NUM.INFLUENZA.DEATHS <-  as.numeric(as.character(all.reports$NUM.INFLUENZA.DEATHS)) 

all.reports$NUM.PNEUMONIA.DEATHS <-
  gsub(',','',all.reports$NUM.PNEUMONIA.DEATHS) 
all.reports$NUM.PNEUMONIA.DEATHS <-  as.numeric(as.character(all.reports$NUM.PNEUMONIA.DEATHS)) 

all.reports$TOTAL.DEATHS <-
  gsub(',','',all.reports$TOTAL.DEATHS) 
all.reports$TOTAL.DEATHS <-  as.numeric(as.character(all.reports$TOTAL.DEATHS))

all.reports2 <- all.reports[,c("week.death","SUB.AREA" ,"report.date","TOTAL.DEATHS") ]
```

```{r}
compare.m <- melt(all.reports2[,c("week.death","SUB.AREA" ,"report.date","TOTAL.DEATHS" )], id.vars  =c("week.death",'SUB.AREA',"report.date"))

names(compare.m) <- 
  c('death_date','state','report_date','variable','N_deaths')

compare.m$report_date <-as.Date(compare.m$report_date)
compare.m$death_date <-as.Date(compare.m$death_date)

compare.m$complete.weeks <- round(as.vector(difftime(compare.m$report_date,
          compare.m$death_date, units='weeks')))

compare.c <- dcast( compare.m[,c('death_date','state','report_date','N_deaths')], state+ death_date ~ report_date, value.var='N_deaths', fun.aggregate = mean)


```


## Set up reporting triangle
```{r}
#Ignore reports from first week
#Filter if death date hasn't happened yet or happened a week ago
#Ignore reports from first week
compare.m.alt <-   compare.m
compare.m.alt <-
  compare.m.alt[compare.m.alt$complete.weeks<=26,]
compare.m.alt <-
  compare.m.alt[as.numeric((compare.m.alt$report_date - compare.m.alt$death_date)) > 6,]

#Only use death data from Apr 1 onwards to account for changes in reporting during pandemic
compare.m.alt <- 
  compare.m.alt[compare.m.alt$death_date >=as.Date('2020-03-29'),]

#Split by state and death date
compare.m.alt.spl <- 
  split(compare.m.alt, paste0(compare.m.alt$state,compare.m.alt$death_date)) 

#Sort by report date
compare.m.alt.spl <- 
  lapply(compare.m.alt.spl, function(x){
    x[order(x$report_date),]
    })

#Subtract previous value to calculate new value
compare.m.alt.spl <- lapply(compare.m.alt.spl, function(x){
  
  #Remove redundant rows if deaths was >6 months ago
    #If no observations before 26 weeks, just take last obs
  if(max(x$complete.weeks)==26 & min(x$complete.weeks)==26){
    x <-   x[nrow(x),] #take most complete observation
    x$complete.weeks <-999
  }
  
    #If some observations before and after 26 weeks

  if(max(x$complete.weeks)==26 & min(x$complete.weeks)<=26){
    x$wk26 <-0
    x$wk26[x$complete.weeks>=26] <-1
    x.early<-x[x$wk26==0,]  
    x.late <-x[x$wk26==1,]  
    x.late <-   x.late[nrow(x.late),, drop=F] #take most complete observation
    #x.late$complete.weeks <- 999
    x <- rbind.data.frame(x.early, x.late)
    x$wk26 <- NULL
  }
  
#Combine first and second weeks of observation
  x$complete.weeks[x$complete.weeks==1] <-2
  
  x$new.reports <- x$N_deaths
  x$first.week.measured<- x$complete.weeks[1]

  if(nrow(x)>1){
    
    #If first observation we have for a week s after week 2, put in the 'most complete' column 999
    if(x$complete.weeks[1]>2){
      x$complete.weeks[1] <-999
      }
    
    for(i in 2:nrow(x)){
       x$new.reports[i]<-x$N_deaths[i]-x$N_deaths[(i-1)]
    }
  }
  return(x)
}) 

compare.m.alt.spl2 <-
  do.call(rbind.data.frame,compare.m.alt.spl)
#replace negative values with 0s
compare.m.alt.spl2$new.reports[compare.m.alt.spl2$new.reports<0] <-0

us.test <- compare.m.alt.spl2[compare.m.alt.spl2$state=='US',]
unique(us.test$first.week.measured)
unique(us.test$complete.weeks)

```


## read in model
```{r}
source('./functions/jags_negbin4.R')


```


```{r}
exclude.states <- c('Connecticut','North Carolina','US')
compare.m.alt.spl2 <- compare.m.alt.spl2[!(compare.m.alt.spl2$state %in% exclude.states),]
```

```{r}
jags.func <- function(state.select){
ds.select <- compare.m.alt.spl2[compare.m.alt.spl2$state==state.select,]

ds.select$complete.weeks[ds.select$complete.weeks>26] <-999

ds.select$first.week.measured[ds.select$first.week.measured>26] <-26

ds7 <- acast( ds.select[,c('death_date','complete.weeks',"new.reports", 'first.week.measured' )],  death_date ~ complete.weeks, value.var="new.reports" , fun.aggregate = sum)

ds7.length <- acast( ds.select[,c('death_date','complete.weeks',"new.reports", 'first.week.measured' )],  death_date ~ complete.weeks, value.var="new.reports" , fun.aggregate = length)

#If we don't have observation, set to NA
ds7[ds7.length==0] <-NA

#what is first week for which we have a report for each state? For dates where we don't have >26 weeks of data and don't observed the first week, we need to sum th eprobabilities (betas)

first.measured.date <- as.data.frame(unique(ds.select[,c('death_date','first.week.measured')]))

max.measured <- max(first.measured.date$first.week.measured[first.measured.date$first.week.measured!=999])

first.measured.date$first.week.measured[first.measured.date$first.week.measured==999] <- max.measured+1

#get rid of columns where we don't have a meausurement
#ds7 <- ds7[,c(1:max.measured, dim(ds7)[2])]

st1 <- ds7
date.sum <- apply(st1,1,sum, na.rm=T) #total observations for 
death_date <- as.Date(dimnames(st1)[[1]])
death_yr <- year(death_date)
death_week <- week(death_date)
death_epiyr <- death_yr

death_epiyr[death_week<=26] <- 
  death_yr[death_week<=26]-1
death_epiyr.index <-as.numeric(as.factor(death_epiyr))

st2 <- st1[, -ncol(st1)] #remove column '99'

beta.priors <- rep(0.1, times=(ncol(st2)))

reporting.triangle <- st1

reporting.triangle[reporting.triangle[,1]>0, ncol(reporting.triangle)] <-NA

max_D <- ncol(reporting.triangle)-1

##############################################################
#Model Fitting
##############################################################
inits1=list(".RNG.seed"=c(123), ".RNG.name"='base::Wichmann-Hill')
inits2=list(".RNG.seed"=c(456), ".RNG.name"='base::Wichmann-Hill')
inits3=list(".RNG.seed"=c(789), ".RNG.name"='base::Wichmann-Hill')


##############################################
#Model Organization
##############################################
model_spec<-textConnection(model_string_negbin4)
model_jags<-jags.model(model_spec, 
                       inits=list(inits1,inits2, inits3),
                       data=list('n.dates' =
                                   nrow(reporting.triangle),
                                 'n' = reporting.triangle,
                                 'D' = ncol(reporting.triangle)-1,
                                                    
                                 
                                 alphat.shape.prior=0.001,
                                alphat.rate.prior=0.001,
                'N.first.obs'=(first.measured.date$first.week.measured-1),
                           'beta.priors'=beta.priors
                                 
                                 ),
                       n.adapt=5000, 
                       n.chains=3)

params<-c('sum.n','sum.lambda',
         'beta.logged', 'alpha','sum.beta')

##############################################
#Posterior Sampling
##############################################
posterior_samples<-coda.samples(model_jags, 
                                params, 
                                n.iter=5000)
posterior_samples.all<-do.call(rbind,posterior_samples)
#post1.summary<-summary(posterior_samples)
#post_means<-colMeans(posterior_samples.all)
out.list=list('posterior_samples.all'=posterior_samples.all,'date.sum'=date.sum,'death_date'=death_date)
return(out.list)
}
```

```{r, eval=T}
#states.test <-'US'
#mod1 <-jags.func('US')

n_cores<- detectCores() -1
states.test <- unique(compare.m.alt.spl2$state)
#states.test <- c('New York', 'California', 'Louisiana', 'Georgia', 'Florida','Kentucky')

  cl <- makeCluster(n_cores)
  clusterEvalQ(cl, {
    library(lubridate, quietly = TRUE)
    library(reshape2, quietly = TRUE)
    library(rjags, quietly = TRUE)
})
  clusterExport(cl, c('jags.func','compare.m.alt.spl2','model_string_negbin4'), environment())

mod1<-pblapply(cl = cl,X=states.test,FUN=jags.func)

stopCluster(cl)
names(mod1) <- states.test
saveRDS(mod1,'./jags_results/mod1.rds')


```


```{r}
mod1<- readRDS('./jags_results/mod1.rds')

jags_extract <- function(ds){

  posterior_samples.all <- ds$posterior_samples.all
  
  death_date <- ds$death_date
  
  date.sum <- ds$date.sum
  
  post_means<-apply(posterior_samples.all, 2, median)
  sample.labs<-names(post_means)
  
  ci<-t(hdi(posterior_samples.all, credMass = 0.95))
  ci<-matrix(sprintf("%.1f",round(ci,1)), ncol=2)
  row.names(ci)<-sample.labs
  post_means<-sprintf("%.1f",round(post_means,1))
  names(post_means)<-sample.labs
  
  #pred.index <- grep('sum.lambda',sample.labs)
  pred.index2 <- grep('sum.n',sample.labs)
  
  beta.index <- grep('beta.logged',sample.labs)
  alpha.index <- grep('alpha',sample.labs)
  baseline.index <- grep('baseline.n',sample.labs)
  sum.lambda.index <- grep('sum.lambda',sample.labs)
  
  pred.means <- post_means[pred.index2]
  lambda.means <- post_means[sum.lambda.index]
  
  #lambda.samps <-
  #  t(posterior_samples.all[,sum.lambda.index])
   pred.samps <-
    t(posterior_samples.all[,pred.index2])
   
  #divide obs and predicted
  complete.prop.est.iter <- 
    apply(pred.samps,2,function(x) date.sum/x) 
  
  complete.prop.est.iter <- 
    cbind.data.frame('n.weeks.since.death'=(nrow(complete.prop.est.iter)+1):2, complete.prop.est.iter)
  
  beta.log.means <- as.numeric(post_means[beta.index])
  probs <- exp(beta.log.means)
  probs.samps <- posterior_samples.all[,beta.index]
  alpha <- as.numeric(post_means[alpha.index])
  
  pred.ci <- ci[pred.index2,]
  all.preds<- as.numeric(as.character(cbind(pred.means, pred.ci)))
  all.preds<- as.data.frame(matrix(all.preds,ncol=3))
  names(all.preds) <- c('pred.med','pred.lcl','preds.ucl')
  
  all.preds <- cbind.data.frame(all.preds, 'obs'=date.sum, 'death_date'=death_date)
  
  out.list=list('preds'=all.preds,'probs'=probs,'probs.samps'=probs.samps,'complete.prop.est.iter'=complete.prop.est.iter)
  return(out.list)
}

res1 <- lapply(mod1, jags_extract)
all.preds1 <- lapply(res1,'[[', 'preds')
complete.prop.est.iter <- lapply(res1,'[[', 'complete.prop.est.iter')


for(i in 1:length(all.preds1)){
  all.preds1[[i]]$state <- names(all.preds1)[i]
  complete.prop.est.iter[[i]]$state <- names(all.preds1)[i]
}

all.preds.df <- do.call('rbind.data.frame', all.preds1)
complete.prop.est.iter.df <- do.call('rbind.data.frame', complete.prop.est.iter)

write.csv(all.preds.df,'./outputs/NobBs.preds.csv')
saveRDS(complete.prop.est.iter.df,'./outputs/NobBs.complete.prop.est.iter.df.rds')

prop.report.wk <- sapply(res1,'[[', 'probs')
complete.wk <- apply(prop.report.wk, 2,cumsum)
complete.wk <- cbind.data.frame('week'=2:(nrow(complete.wk)+1), complete.wk)
write.csv(complete.wk,'./outputs/NobBs.complete.csv')

matplot(complete.wk$week, complete.wk[,-1], type='l')

prop.report.wk.iter <- sapply(res1,'[[', 'probs.samps', simplify='array')
complete.wk.iter <- apply(prop.report.wk.iter, c(1,3),function(x) cumsum(exp(x)))
dimnames(complete.wk.iter)[[1]] <- 2:(dim(complete.wk.iter)[1]+1)
saveRDS(complete.wk.iter,'./outputs/NobBs.complete.iters.rds')

#this is based on obs/exp from NobBS
# complete.prop.est.median.df <-
#   cbind.data.frame(complete.prop.est.iter.df[,1], apply(complete.prop.est.iter.df[,-1],1, median))
# 
# complete.wk.iter.median <- apply(complete.wk.iter, c(1,3), median)


```


```{r}

par(mfrow=c(4,4), mar=c(2,2,1,1))
for(i in 1:length(all.preds1)){
ds.plot <- all.preds1[[i]]
ds.plot <- ds.plot[ds.plot$death_date>=as.Date('2020-01-01'),]
state <- names(all.preds1)[i]

yrange <- range(c(ds.plot$pred.med,ds.plot$obs))
plot(ds.plot$death_date,ds.plot$pred.med, type='l', lty=c(1,2,2), col=c('black','gray','gray'), ylim=yrange, main=state)
points(ds.plot$death_date,ds.plot$obs, type='l',pch=16 ,col='red', lty=2)
polygon(c(ds.plot$death_date, rev(ds.plot$death_date)), c(ds.plot$pred.lcl, rev(ds.plot$preds.ucl)), col=rgb(0,0,1,alpha=0.1), border=F )
abline(v=as.Date('2020-04-18'), lty=2, col='gray')

}  


```

## As a check, let's look at the backfilling patterns for certain states
```{r, fig.width=6, fig.height=8}

check1.m <- melt(compare.m.alt[,c('state','death_date','report_date', 'N_deaths')], id.vars=c('state','death_date','report_date'))
check1.m$time.report <- as.numeric(check1.m$report_date - check1.m$death_date)
check1.m <- check1.m[check1.m$time.report>7,]
check1.c <- acast(check1.m, state~death_date~report_date)
dates1 <- as.Date(dimnames(check1.c)[[2]])
check2 <- check1.c[,dates1 >=as.Date('2020-01-01'),]
dates2 <- dates1[dates1 >=as.Date('2020-01-01')]

#plot.states <- c('Delaware', 'Louisiana', 'New York', 'Florida')
plot.states <- c('Delaware', 'Louisiana', 'New York', 'Florida', 'South Carolina','Washington' )
par(mfrow=c(3,2))
for(i in plot.states){
  check3 <- check2[i,,]
  dates1 <- as.Date(dimnames(check3)[[1]])
check3 <- check3[dates1 >=as.Date('2020-01-01'),]
dates2 <- dates1[dates1 >=as.Date('2020-01-01')]
  dates2 <- dates1[dates1 >=as.Date('2020-01-01')]
matplot(dates2 ,check3, type='l', bty='l', xaxt='n', xlab='', ylab='Provisional Count', main=i, ylim=c(0, max(check3, na.rm=T)))
axis(side=1, at=seq.Date(from=min(dates2), to=max(dates2), length.out=length(dates2)), label=dates2)
abline(v=as.Date('2020-04-01'))
}
```