library(ggplot2)
library(zoo)

mav=\(x,y)sapply(1:length(x),\(i)mean(x[max(0,i-y):min(length(x),i+y)],na.rm=T))

download.file("https://healthdata.gov/api/views/j8mb-icvb/rows.csv","statespcr.csv") # https://healthdata.gov/dataset/COVID-19-Diagnostic-Laboratory-Testing-PCR-Testing/j8mb-icvb
download.file("https://data.cdc.gov/api/views/xkkf-xrst/rows.csv","statesexcess.csv") # https://www.cdc.gov/nchs/nvss/vsrr/covid19/excess_deaths.htm

pcr=read.csv("statespcr.csv")
excess=read.csv("statesexcess.csv")

for(state in intersect(pcr$state_name,excess$State)){
  pcr2=pcr[pcr$state_name==state,]
  pcr3=pcr2[pcr2$overall_outcome=="Positive",c(6,7)]|>merge(pcr2[pcr2$overall_outcome=="Negative",c(6,7)],by=1,all=T)

  pos=mav(pcr3[,2],10)
  neg=mav(pcr3[,3],10)
  df4=data.frame(date=as.Date(pcr3$date),pcrrate=100*pos/(pos+neg))
  df4$pcrrate[mav(rowSums(pcr3[,-1]),10)<30]=NA

  excess2=excess[excess$State==state&excess$Outcome=="All causes"&excess$Year>=2020&excess$Week.Ending.Date<="2023-05-01",]
  df5=data.frame(date=as.Date(excess2$Week.Ending.Date)-3,excess=100*(excess2$Observed.Number/excess2$Average.Expected.Count-1))

  xy=merge(df5,df4,all=T)
  nonna=1:tail(which(!is.na(xy$excess)),1)
  xy$excess[nonna]=mav(zoo::na.approx(xy$excess[nonna],na.rm=F),3)

  xstart=as.Date("2020-01-01")
  xend=as.Date("2023-05-01")
  xy=xy[xy$date>=xstart&xy$date<=xend,]

  # candidates=c(sapply(c(1,2,5),\(x)x*10^c(-2:5)))
  # ystep=candidates[which.min(abs(candidates-max(xy[,2],na.rm=T)/6))]
  # ystep2=candidates[which.min(abs(candidates-max(xy[,3],na.rm=T)/6))]
  # ystart=min(0,ystep*floor(min(xy[,2],na.rm=T)/ystep))
  # yend=ystep*ceiling(max(xy[,2],na.rm=T)/ystep)
  # yend2=ystep2*ceiling(max(xy[,3],na.rm=T)/ystep2)
  ystart=-20;yend=120;yend2=50;ystep=20;ystep2=10

  secmult=yend/yend2
  seccolor=hcl(135,70,50)

  ggplot(xy,aes(x=date,y=excess))+
  geom_hline(yintercept=c(ystart,0,yend),color="gray75",linewidth=.3,lineend="square")+
  geom_vline(xintercept=c(xstart,xend),color="gray75",linewidth=.3,lineend="square")+
  geom_line(linewidth=.3)+
  geom_line(aes(y=secmult*pcrrate),linewidth=.3,color=seccolor)+
  coord_cartesian(clip="off")+
  scale_x_date(limits=c(xstart,xend),breaks=seq(xstart,xend,"2 months"),expand=expansion(0),date_labels="%b 1 %y")+
  scale_y_continuous(limits=c(ystart,yend),breaks=seq(ystart,yend,ystep),expand=expansion(0),sec.axis=sec_axis(trans=~./secmult,breaks=seq(0,yend2,ystep2),name="Percentage of positive PCR tests (21-day moving average)"))+
  labs(y="Excess mortality percent (21-day moving average)")+
  ggtitle(paste0(state," (r=",sprintf("%.2f",cor(xy[,2],xy[,3],use="complete.obs")),")"))+
  theme(
    axis.text=element_text(size=6,color="black"),
    axis.text.x=element_text(angle=90,vjust=.5,hjust=1),
    axis.ticks=element_line(linewidth=.3,color="gray75"),
    axis.ticks.length=unit(.2,"lines"),
    axis.title=element_text(size=7.5),
    axis.title.y.right=element_text(color=seccolor,margin=margin(0,0,0,5)),
    axis.text.y.right=element_text(color=seccolor),
    axis.title.x=element_blank(),
    panel.background=element_rect(fill="white"),
    panel.grid=element_blank(),
    plot.subtitle=element_text(size=7),
    plot.title=element_text(size=9)
  )

  ggsave(paste0(state,".png"),width=5,height=3.5)
}