In weinbergerlab/ExcessILI: A package and RShiny app to examine excess ILI from syndromic surveillance data
knitr::opts_chunk$set(
collapse = TRUE,
echo=F,
warning=FALSE,
message=FALSE,
comment = "#>"
)
library(ExcessILI)
library(cdcfluview)
library(reshape2)
library(ggplot2)
library(lubridate)
library(RColorBrewer)
library(plotly)
library(MMWRweek)
library(readr)
Overview
In this example, we will analyze increases of pneumonia and influenza (P&I) above the seasonal baseline, adjusting for year-to-year variations in incidence. We use the weekly data published by the CDC We fit a simple seasonal baseline using harmonics, and adjust for year-to-year variations. Influenza activity is adjusted for as well using data from NREVSS (percent of tests positive for influenza in the previous week). The model is fit through end-of-February 2020 and then extrapolated forward based on the time of year.
## Download the mortality data
pi.data <- pi_mortality(coverage_area='state')
#download the NREVSS data
nrevvs.state <- cdcfluview::who_nrevss(region = c("state"))
clin <- nrevvs.state[["clinical_labs"]]
#data(cdcfluview::hhs_regions)
cw.file <- cdcfluview::hhs_regions
clin2 <- merge(clin, cw.file,
by.x = "region",
by.y = "state_or_territory")
clin2.subsetvars <-
c('region', 'region_number',
'year', 'week', 'wk_date',
'total_a','total_b',
'total_specimens')
clin2 <- clin2[, clin2.subsetvars]
names(clin2)[1:2] <- c("state", "hhs_region")
nrevvs_hhs <- cdcfluview::who_nrevss(region = c("hhs"))
clin.hhs <- nrevvs_hhs[["clinical_labs"]]
clin.hhs.subsetvars <-
c('region',
'wk_date',
"total_a",'total_b',
'total_specimens')
clin.hhs <- clin.hhs[, clin.hhs.subsetvars]
clin.hhs$region <- as.numeric(gsub("Region ", "", clin.hhs$region))
names(clin.hhs) <-
c("hhs_region",
"wk_date",
"hhs_total_a",'hhs_total_b',
'hhs_total_specimens')
clin3 <- merge(clin2, clin.hhs,
by = c("hhs_region", "wk_date"))
clin3$total_a[is.na(clin3$total_a)] <-
clin3$hhs_total_a[is.na(clin3$total_a)]
clin3$total_b[is.na(clin3$total_b)] <-
clin3$hhs_total_b[is.na(clin3$total_b)]
clin3$total_specimens[is.na(clin3$total_specimens)] <-
clin3$hhs_total_specimens[is.na(clin3$total_specimens)]
clin3$state.abb <- state.abb[match(clin3$state, state.name)]
names(clin3) <-
c("hhs_region",
"wk_date",
"state_name",
"MMWRyear", "MMWRweek",
"total_a",'total_b',
'total_specimens',
'total_a_hhs', "total_b_hhs",
'total_specimens_hhs',
"state")
clin3$total_a <- as.numeric(clin3$total_a)
clin3$total_b <- as.numeric(clin3$total_b)
clin3$total_specimens <- as.numeric(clin3$total_specimens)
clin3$flu_pct_adj <- (clin3$total_a + clin3$total_b + 0.5) /
(clin3$total_specimens + 0.5)
clin3$fluN <- clin3$total_a + clin3$total_b + 0.5
clin3$flu.var <- clin3$flu_pct_adj
clin4<-clin3[,c('state','flu_pct_adj', 'wk_date')]
clin4.lag1<-clin4
clin4.lag1$wk_date <- clin4$wk_date + days(7)
names(clin4.lag1) <-c('state','flu_pct_adj_lag1','wk_date')
clin4.lag2<-clin4
clin4.lag2$wk_date <- clin4$wk_date + days(14)
names(clin4.lag2) <-c('state','flu_pct_adj_lag2','wk_date')
clin4.lags <- merge(clin4, clin4.lag1, by=c('state','wk_date'))
clin4.lags <- merge(clin4.lags, clin4.lag2, by=c('state','wk_date'))
#Format and fill mssings with 0s
pi.data$state <- state.abb[match(pi.data$region_name, state.name)]
pi.data$state[pi.data$region_name == 'New York City'] <- 'NYC'
spl1<-split(pi.data, pi.data$state)
min.state <- lapply(spl1, function(x){ x$miss.x<-min(x$total_pni)
return(x)
})
pi.data.clean <- do.call('rbind.data.frame',min.state)
pi.data.clean <- pi.data.clean[!is.na(pi.data.clean$miss.x),]
pi.data.clean2<- merge( pi.data.clean,clin4.lags, by.x=c('week_start', 'state'), by.y=c('wk_date','state'))
pi.data.clean2<- pi.data.clean2[order(pi.data.clean2$state, pi.data.clean2$week_end),]
#What is most appropriate lag to use for NREVSS data?
pi.data.clean2.spl <- split(pi.data.clean2, pi.data.clean2$state)
cor.lags <- sapply(pi.data.clean2.spl, function(x){
cor(x[,c("percent_pni","flu_pct_adj",'flu_pct_adj_lag1','flu_pct_adj_lag2')])
}, simplify='array')
matplot(cor.lags['percent_pni',-1,], type='l')
#Is lag of 0,1,or 2 the best?...shows 1 weeks lag best in 20 states, 2 week lag best in 3 states
table(apply(cor.lags['percent_pni',-1,],2, function(x) which(max(x)==x)))
#Run analysis
excess_deaths1.adjusted <-
excessCases(ds = pi.data.clean2,
datevar = "week_start",
statevar = "state",
denom.var = "all_deaths",
adj.flu = "flu_pct_adj_lag1",
#covs=c("flu_pct_adj", "flu_pct_adj_lag1", "flu_pct_adj_lag2"),
use.syndromes = c("total_pni"),
extrapolation.date = "2020-03-01",
time.res='week')
#dashboardPlot(excess_deaths1.adjusted)
### Extract the quantities of interest
#Which syndrome do you want to plot, and over what time range?
syndrome.select <- 'total_pni' #which syndrome do you want to plot?
n.days<-52 #How many days to plot?
ds <- excess_deaths1.adjusted
#Extract the data needed to plot from the results
dates1 <-
ds[[1]][[1]][[1]]$date
unexplained.cases <-
excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "unexplained.cases")
unexplained.log.rr <-
excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "resid1")
denom <-
excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "denom")
upi <-
excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "upi")
lpi <-
excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "lpi")
obs <-
excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "y")
pred<- excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "pred")
rr <- excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "resid1")
excess_deaths <- excessExtract(ds = ds,
syndrome = syndrome.select,
extract.quantity = "unexplained.cases")
Observed weekly death rate vs seasonal baseline (+/-95% Prediction Interval)
n.days <- 52
select.indices <- (length(dates1)-n.days):length(dates1)
dates<-dates1[select.indices]
states <- dimnames(pred)[[2]]
ages <- dimnames(pred)[[3]]
The black line shows the observed proportion of deaths that were due to Pneumonia & Influenza (P&I) per week. The red line and shaded area represent the 95% Prediction Interval. The latest data is for the week ending r max(dates1)+7.
par(mfrow=c(4,4))
for(i in 1:dim(pred)[2]){
for(j in 1:dim(pred)[3]){
y.range1<-range(c( pred[select.indices,i,j]/denom[select.indices,i,j],obs[select.indices,i,j]/denom[select.indices,i,j], upi[select.indices,i,j]/denom[select.indices,i,j],0))
plot(dates,
pred[select.indices,i,j]/denom[select.indices,i,j],
type='l',
col='red',
ylim=y.range1,
bty='l',
ylab='Proportion',
main=paste(states[i],ages[j]))
points(dates,
obs[select.indices,i,j]/denom[select.indices,i,j],
type='l',
col='black')
polygon(c(dates,
rev(dates)),
c(lpi[select.indices,i,j]/denom[select.indices,i,j],
rev(upi[select.indices,i,j]/denom[select.indices,i,j])),
col = rgb(1, 0, 0, alpha = 0.1),
border = NA)
}
}
Observed deaths/expected deaths by state
rr2<-rr[,,1]
date.mmwrdates <- mmwr_week(dates1)
mmwr.epiyr<- date.mmwrdates$mmwr_year
mmwr.epiyr[date.mmwrdates$mmwr_week<=26] <- mmwr.epiyr[date.mmwrdates$mmwr_week<=26] - 1
mmwr.epiwk <- date.mmwrdates$mmwr_week
mmwr.epiwk[date.mmwrdates$mmwr_week>=27]<-date.mmwrdates$mmwr_week[date.mmwrdates$mmwr_week>=27] - 52
mmwr.epiwk <- mmwr.epiwk +26
check<-cbind.data.frame(date.mmwrdates,mmwr.epiwk, mmwr.epiyr)
These plots show the Observed/Expected number of deaths due to pneumonia and influenza in each week for the 2019-20 year (red) compared to previous years (gray). Values close to 1 indicate that the values for that week are close to what would be expected based on the time of year and influenza activity.
par(mfrow=c(4,4))
for(i in 1:dim(pred)[2]){
y.range1<-c(0,2)
ds2<-cbind.data.frame('epiwk'=mmwr.epiwk,'epiyr'=mmwr.epiyr, rr=rr2[,i])
ds2.c<-dcast(ds2, epiwk~epiyr, value.var='rr', fun.aggregate = mean)
cols1<-c(rep('grey',(ncol(ds2.c)-2) ),'red')
matplot(ds2.c$epiwk ,
exp(ds2.c[,-1]),
type='l',
col=cols1,
ylim=y.range1,
bty='l',
lty=1,
ylab='Observed/Expected',
main=paste(states[i]))
abline(h=1, col='black')
}
Compare Excess P&I mortality vs Excess ILI
Here we compare the observed vs expected number of deaths due to pneumonia and influenza in each week compare to the observed vs expected number of outpatient visits for influenza-like illness (ILI) in each week. we would expect ILI (blue line) to increase earlier than deaths (red line)
ili.data <- ilinet(region = c("state"))
ili.data$state <- state.abb[match(ili.data$region, state.name)]
ili.data <- ili.data[, c("state", "week_start", "ilitotal", "total_patients")]
ili.data <- ili.data[!is.na(ili.data$total_patients),]
ili.data.spl <- split(ili.data, ili.data$state)
min<-sapply(ili.data.spl, function(x) min(x$total_patients))
state.select<-names(min)[which(min>0) ]
ili.data <- ili.data[ili.data$state %in% state.select,]
## Run the main analysis function, adjusting for flu using NREVSS data
excess_cases1 <-
excessCases(ds = ili.data,
datevar = "week_start",
statevar = "state",
denom.var = "total_patients",
adj.flu = "auto",
use.syndromes = c("ilitotal"),
extrapolation.date = "2020-03-01",
time.res='week')
dates.ili <-
excess_cases1[[1]][[1]][[1]]$date
rr.ili <- excessExtract(ds = excess_cases1,
syndrome = "ilitotal",
extract.quantity = "resid1")
rr2.ili <- rr.ili[,,1]
date.mmwrdates.ili <- mmwr_week(dates.ili)
mmwr.epiyr.ili<- date.mmwrdates.ili$mmwr_year
mmwr.epiyr.ili[date.mmwrdates.ili$mmwr_week<=26] <- mmwr.epiyr.ili[date.mmwrdates.ili$mmwr_week<=26] - 1
mmwr.epiwk.ili <- date.mmwrdates.ili$mmwr_week
mmwr.epiwk.ili[date.mmwrdates.ili$mmwr_week>=27]<-date.mmwrdates.ili$mmwr_week[date.mmwrdates.ili$mmwr_week>=27] - 52
mmwr.epiwk.ili <- mmwr.epiwk.ili +26
common.states <- intersect(colnames(rr2), colnames(rr2.ili))
rr2.comp <- rr2[,common.states]
rr2.ili.comp <- rr2.ili[,common.states]
par(mfrow=c(4,4))
for(i in 1:length(common.states)){
y.range1<-c(0,2)
ds2<-cbind.data.frame('epiwk'=mmwr.epiwk,'epiyr'=mmwr.epiyr, rr=rr2.comp[,i])
ds2.c<-dcast(ds2, epiwk~epiyr, value.var='rr', fun.aggregate = mean)
cols1<-c(rep('grey',(ncol(ds2.c)-2) ),'red')
plot(ds2.c$epiwk ,
exp(ds2.c[,'2019']),
type='l',
col='red',
ylim=y.range1,
bty='l',
lty=1,
ylab='Observed/Expected',
main=paste(common.states[i], ' Deaths(red),','ILI(blue)' ))
es2<-cbind.data.frame('epiwk'=mmwr.epiwk.ili,'epiyr'=mmwr.epiyr.ili, rr=rr2.ili.comp[,i])
es2.c<-dcast(es2, epiwk~epiyr, value.var='rr', fun.aggregate = mean)
points(es2.c$epiwk ,
exp(es2.c[,'2019']), type='l', col='blue')
abline(h=1, col='black')
}
## excess deaths
excess_deaths2 <- excess_deaths[dates1 >= as.Date('2020-02-02'),,1]
excess_deaths.state <- apply(excess_deaths2,2,sum)
cumsum_excess_deaths_state <- apply(excess_deaths2,2,cumsum)
matplot(cumsum_excess_deaths_state, type='l', bty='l')
weinbergerlab/ExcessILI documentation built on May 30, 2021, 10:57 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, echo=F, warning=FALSE, message=FALSE, comment = "#>" )
library(ExcessILI) library(cdcfluview) library(reshape2) library(ggplot2) library(lubridate) library(RColorBrewer) library(plotly) library(MMWRweek) library(readr)
Overview
In this example, we will analyze increases of pneumonia and influenza (P&I) above the seasonal baseline, adjusting for year-to-year variations in incidence. We use the weekly data published by the CDC We fit a simple seasonal baseline using harmonics, and adjust for year-to-year variations. Influenza activity is adjusted for as well using data from NREVSS (percent of tests positive for influenza in the previous week). The model is fit through end-of-February 2020 and then extrapolated forward based on the time of year.
## Download the mortality data pi.data <- pi_mortality(coverage_area='state')
#download the NREVSS data nrevvs.state <- cdcfluview::who_nrevss(region = c("state")) clin <- nrevvs.state[["clinical_labs"]] #data(cdcfluview::hhs_regions) cw.file <- cdcfluview::hhs_regions clin2 <- merge(clin, cw.file, by.x = "region", by.y = "state_or_territory") clin2.subsetvars <- c('region', 'region_number', 'year', 'week', 'wk_date', 'total_a','total_b', 'total_specimens') clin2 <- clin2[, clin2.subsetvars] names(clin2)[1:2] <- c("state", "hhs_region") nrevvs_hhs <- cdcfluview::who_nrevss(region = c("hhs")) clin.hhs <- nrevvs_hhs[["clinical_labs"]] clin.hhs.subsetvars <- c('region', 'wk_date', "total_a",'total_b', 'total_specimens') clin.hhs <- clin.hhs[, clin.hhs.subsetvars] clin.hhs$region <- as.numeric(gsub("Region ", "", clin.hhs$region)) names(clin.hhs) <- c("hhs_region", "wk_date", "hhs_total_a",'hhs_total_b', 'hhs_total_specimens') clin3 <- merge(clin2, clin.hhs, by = c("hhs_region", "wk_date")) clin3$total_a[is.na(clin3$total_a)] <- clin3$hhs_total_a[is.na(clin3$total_a)] clin3$total_b[is.na(clin3$total_b)] <- clin3$hhs_total_b[is.na(clin3$total_b)] clin3$total_specimens[is.na(clin3$total_specimens)] <- clin3$hhs_total_specimens[is.na(clin3$total_specimens)] clin3$state.abb <- state.abb[match(clin3$state, state.name)] names(clin3) <- c("hhs_region", "wk_date", "state_name", "MMWRyear", "MMWRweek", "total_a",'total_b', 'total_specimens', 'total_a_hhs', "total_b_hhs", 'total_specimens_hhs', "state") clin3$total_a <- as.numeric(clin3$total_a) clin3$total_b <- as.numeric(clin3$total_b) clin3$total_specimens <- as.numeric(clin3$total_specimens) clin3$flu_pct_adj <- (clin3$total_a + clin3$total_b + 0.5) / (clin3$total_specimens + 0.5) clin3$fluN <- clin3$total_a + clin3$total_b + 0.5 clin3$flu.var <- clin3$flu_pct_adj clin4<-clin3[,c('state','flu_pct_adj', 'wk_date')] clin4.lag1<-clin4 clin4.lag1$wk_date <- clin4$wk_date + days(7) names(clin4.lag1) <-c('state','flu_pct_adj_lag1','wk_date') clin4.lag2<-clin4 clin4.lag2$wk_date <- clin4$wk_date + days(14) names(clin4.lag2) <-c('state','flu_pct_adj_lag2','wk_date') clin4.lags <- merge(clin4, clin4.lag1, by=c('state','wk_date')) clin4.lags <- merge(clin4.lags, clin4.lag2, by=c('state','wk_date'))
#Format and fill mssings with 0s pi.data$state <- state.abb[match(pi.data$region_name, state.name)] pi.data$state[pi.data$region_name == 'New York City'] <- 'NYC' spl1<-split(pi.data, pi.data$state) min.state <- lapply(spl1, function(x){ x$miss.x<-min(x$total_pni) return(x) }) pi.data.clean <- do.call('rbind.data.frame',min.state) pi.data.clean <- pi.data.clean[!is.na(pi.data.clean$miss.x),] pi.data.clean2<- merge( pi.data.clean,clin4.lags, by.x=c('week_start', 'state'), by.y=c('wk_date','state')) pi.data.clean2<- pi.data.clean2[order(pi.data.clean2$state, pi.data.clean2$week_end),]
#What is most appropriate lag to use for NREVSS data? pi.data.clean2.spl <- split(pi.data.clean2, pi.data.clean2$state) cor.lags <- sapply(pi.data.clean2.spl, function(x){ cor(x[,c("percent_pni","flu_pct_adj",'flu_pct_adj_lag1','flu_pct_adj_lag2')]) }, simplify='array') matplot(cor.lags['percent_pni',-1,], type='l') #Is lag of 0,1,or 2 the best?...shows 1 weeks lag best in 20 states, 2 week lag best in 3 states table(apply(cor.lags['percent_pni',-1,],2, function(x) which(max(x)==x)))
#Run analysis excess_deaths1.adjusted <- excessCases(ds = pi.data.clean2, datevar = "week_start", statevar = "state", denom.var = "all_deaths", adj.flu = "flu_pct_adj_lag1", #covs=c("flu_pct_adj", "flu_pct_adj_lag1", "flu_pct_adj_lag2"), use.syndromes = c("total_pni"), extrapolation.date = "2020-03-01", time.res='week')
#dashboardPlot(excess_deaths1.adjusted)
### Extract the quantities of interest #Which syndrome do you want to plot, and over what time range? syndrome.select <- 'total_pni' #which syndrome do you want to plot? n.days<-52 #How many days to plot? ds <- excess_deaths1.adjusted
#Extract the data needed to plot from the results dates1 <- ds[[1]][[1]][[1]]$date unexplained.cases <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "unexplained.cases") unexplained.log.rr <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "resid1") denom <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "denom") upi <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "upi") lpi <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "lpi") obs <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "y") pred<- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "pred") rr <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "resid1") excess_deaths <- excessExtract(ds = ds, syndrome = syndrome.select, extract.quantity = "unexplained.cases")
Observed weekly death rate vs seasonal baseline (+/-95% Prediction Interval)
n.days <- 52 select.indices <- (length(dates1)-n.days):length(dates1) dates<-dates1[select.indices] states <- dimnames(pred)[[2]] ages <- dimnames(pred)[[3]]
The black line shows the observed proportion of deaths that were due to Pneumonia & Influenza (P&I) per week. The red line and shaded area represent the 95% Prediction Interval. The latest data is for the week ending r max(dates1)+7.
par(mfrow=c(4,4)) for(i in 1:dim(pred)[2]){ for(j in 1:dim(pred)[3]){ y.range1<-range(c( pred[select.indices,i,j]/denom[select.indices,i,j],obs[select.indices,i,j]/denom[select.indices,i,j], upi[select.indices,i,j]/denom[select.indices,i,j],0)) plot(dates, pred[select.indices,i,j]/denom[select.indices,i,j], type='l', col='red', ylim=y.range1, bty='l', ylab='Proportion', main=paste(states[i],ages[j])) points(dates, obs[select.indices,i,j]/denom[select.indices,i,j], type='l', col='black') polygon(c(dates, rev(dates)), c(lpi[select.indices,i,j]/denom[select.indices,i,j], rev(upi[select.indices,i,j]/denom[select.indices,i,j])), col = rgb(1, 0, 0, alpha = 0.1), border = NA) } }
Observed deaths/expected deaths by state
rr2<-rr[,,1] date.mmwrdates <- mmwr_week(dates1) mmwr.epiyr<- date.mmwrdates$mmwr_year mmwr.epiyr[date.mmwrdates$mmwr_week<=26] <- mmwr.epiyr[date.mmwrdates$mmwr_week<=26] - 1 mmwr.epiwk <- date.mmwrdates$mmwr_week mmwr.epiwk[date.mmwrdates$mmwr_week>=27]<-date.mmwrdates$mmwr_week[date.mmwrdates$mmwr_week>=27] - 52 mmwr.epiwk <- mmwr.epiwk +26 check<-cbind.data.frame(date.mmwrdates,mmwr.epiwk, mmwr.epiyr)
These plots show the Observed/Expected number of deaths due to pneumonia and influenza in each week for the 2019-20 year (red) compared to previous years (gray). Values close to 1 indicate that the values for that week are close to what would be expected based on the time of year and influenza activity.
par(mfrow=c(4,4)) for(i in 1:dim(pred)[2]){ y.range1<-c(0,2) ds2<-cbind.data.frame('epiwk'=mmwr.epiwk,'epiyr'=mmwr.epiyr, rr=rr2[,i]) ds2.c<-dcast(ds2, epiwk~epiyr, value.var='rr', fun.aggregate = mean) cols1<-c(rep('grey',(ncol(ds2.c)-2) ),'red') matplot(ds2.c$epiwk , exp(ds2.c[,-1]), type='l', col=cols1, ylim=y.range1, bty='l', lty=1, ylab='Observed/Expected', main=paste(states[i])) abline(h=1, col='black') }
Compare Excess P&I mortality vs Excess ILI
Here we compare the observed vs expected number of deaths due to pneumonia and influenza in each week compare to the observed vs expected number of outpatient visits for influenza-like illness (ILI) in each week. we would expect ILI (blue line) to increase earlier than deaths (red line)
ili.data <- ilinet(region = c("state")) ili.data$state <- state.abb[match(ili.data$region, state.name)] ili.data <- ili.data[, c("state", "week_start", "ilitotal", "total_patients")] ili.data <- ili.data[!is.na(ili.data$total_patients),] ili.data.spl <- split(ili.data, ili.data$state) min<-sapply(ili.data.spl, function(x) min(x$total_patients)) state.select<-names(min)[which(min>0) ] ili.data <- ili.data[ili.data$state %in% state.select,] ## Run the main analysis function, adjusting for flu using NREVSS data excess_cases1 <- excessCases(ds = ili.data, datevar = "week_start", statevar = "state", denom.var = "total_patients", adj.flu = "auto", use.syndromes = c("ilitotal"), extrapolation.date = "2020-03-01", time.res='week') dates.ili <- excess_cases1[[1]][[1]][[1]]$date rr.ili <- excessExtract(ds = excess_cases1, syndrome = "ilitotal", extract.quantity = "resid1") rr2.ili <- rr.ili[,,1] date.mmwrdates.ili <- mmwr_week(dates.ili) mmwr.epiyr.ili<- date.mmwrdates.ili$mmwr_year mmwr.epiyr.ili[date.mmwrdates.ili$mmwr_week<=26] <- mmwr.epiyr.ili[date.mmwrdates.ili$mmwr_week<=26] - 1 mmwr.epiwk.ili <- date.mmwrdates.ili$mmwr_week mmwr.epiwk.ili[date.mmwrdates.ili$mmwr_week>=27]<-date.mmwrdates.ili$mmwr_week[date.mmwrdates.ili$mmwr_week>=27] - 52 mmwr.epiwk.ili <- mmwr.epiwk.ili +26
common.states <- intersect(colnames(rr2), colnames(rr2.ili)) rr2.comp <- rr2[,common.states] rr2.ili.comp <- rr2.ili[,common.states] par(mfrow=c(4,4)) for(i in 1:length(common.states)){ y.range1<-c(0,2) ds2<-cbind.data.frame('epiwk'=mmwr.epiwk,'epiyr'=mmwr.epiyr, rr=rr2.comp[,i]) ds2.c<-dcast(ds2, epiwk~epiyr, value.var='rr', fun.aggregate = mean) cols1<-c(rep('grey',(ncol(ds2.c)-2) ),'red') plot(ds2.c$epiwk , exp(ds2.c[,'2019']), type='l', col='red', ylim=y.range1, bty='l', lty=1, ylab='Observed/Expected', main=paste(common.states[i], ' Deaths(red),','ILI(blue)' )) es2<-cbind.data.frame('epiwk'=mmwr.epiwk.ili,'epiyr'=mmwr.epiyr.ili, rr=rr2.ili.comp[,i]) es2.c<-dcast(es2, epiwk~epiyr, value.var='rr', fun.aggregate = mean) points(es2.c$epiwk , exp(es2.c[,'2019']), type='l', col='blue') abline(h=1, col='black') }
## excess deaths excess_deaths2 <- excess_deaths[dates1 >= as.Date('2020-02-02'),,1] excess_deaths.state <- apply(excess_deaths2,2,sum) cumsum_excess_deaths_state <- apply(excess_deaths2,2,cumsum) matplot(cumsum_excess_deaths_state, type='l', bty='l')
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.