Comments to Ethical Skeptic (part 2) - sars2.net

library(data.table);library(ggplot2)

agecut=\(x,y)cut(x,c(y,Inf),paste0(y,c(paste0("-",y[-1]-1),"+")),T,F)

kim=\(x)ifelse(x>=1e3,ifelse(x>=1e6,paste0(x/1e6,"M"),paste0(x/1e3,"k")),x)

new=fread("https://www2.census.gov/programs-surveys/popest/datasets/2020-2023/national/asrh/nc-est2023-agesex-res.csv")
old=fread("https://www2.census.gov/programs-surveys/popest/datasets/2010-2020/national/asrh/nc-est2020-agesex-res.csv")

old=old[SEX==0&AGE!=999,.(age=AGE,pop=unlist(.SD[,-(1:4)]),year=rep(2010:2020,each=.N))]
new=new[SEX==0&AGE!=999,.(age=AGE,pop=unlist(.SD[,-(1:3)]),year=rep(2020:2023,each=.N))]

older=fread("https://www2.census.gov/programs-surveys/popest/tables/2000-2010/intercensal/national/us-est00int-01.csv")
older=older[6:26,c(1,3:12,14)][,.(age=as.numeric(sub("\\D*(\\d+).*","\\1",sub("Under 5",0,V1))),pop=as.numeric(gsub(",","",unlist(.SD[,-1]))),year=rep(2000:2010,each=.N))]

oldest=fread("https://www2.census.gov/programs-surveys/popest/tables/2000-2009/national/asrh/nc-est2009-01.csv")
oldest=oldest[6:26,c(1,2:11)][,.(age=as.numeric(sub("\\D*(\\d+).*","\\1",sub("Under 5",0,V1))),pop=as.numeric(gsub(",","",unlist(.SD[,-1]))),year=rep(2009:2000,each=.N))]

mult=merge(old[year==2020],new[year==2020,.(age,new=pop)])[,.(ratio=new/pop,age)]
mult=merge(old,mult)[,mult:=(year-2010)/10]
mult=mult[,.(age,year,pop=((1-mult)*1+mult*ratio)*pop)]

won=fread("http://sars2.net/f/wondercanceryearlysingle.csv")[age<85,.(year,pop,age)]
won=rbind(won,fread("http://sars2.net/f/wondercanceryearlyten.csv")[cause==cause[1]&age==85,.(year,pop,age)])

p=cbind(oldest,z="2000-2009")
p=rbind(p,cbind(older,z="2000-2010 intercensal"))
p=rbind(p,cbind(old,z="2010-2020"))
p=rbind(p,cbind(mult,z="2010-2020 multiplied by linear slope"))
p=rbind(p,cbind(new,z="2020-2023"))
p=rbind(p,cbind(won,z="CDC WONDER"))

p[,z:=factor(z,unique(z))]

ages=0:4*20
ages=c(0,25,45,65,75,85)

p=p[,.(pop=sum(pop)),.(z,year,age=agecut(age,ages))]

xstart=2000;xend=2023
p=p[year%in%xstart:xend]

ggplot(p,aes(x=year,y=pop))+
coord_cartesian(clip="off")+
facet_wrap(~age,ncol=2,dir="v",scales="free_y")+
geom_vline(xintercept=c(2009.5,2019.5),color="gray80",linewidth=.23)+
geom_line(aes(color=z,alpha=z),linewidth=.3)+
geom_point(aes(color=z,shape=z,size=z),stroke=.3)+
geom_text(fontface=2,data=p[,max(pop),age],aes(label=paste0("\n   ",age,"   \n"),y=V1),x=xstart-.5,lineheight=.4,hjust=0,vjust=1,size=grid::convertUnit(unit(7,"pt"),"mm"))+
labs(title="Mid-year US resident population estimates by age group",subtitle="Source: www2.census.gov/programs-surveys/popest/datasets and CDC WONDER.",x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart-.5,xend+.5),breaks=seq(xstart-.5,xend+.5,.5),labels=c(rbind("",ifelse(xstart:xend%%2==0,xstart:xend,"")),""),expand=expansion(0))+
scale_y_continuous(labels=kim,breaks=\(x)Filter(\(y)y>min(x)+(max(x)-min(x))*.05,pretty(x,4)),expand=expansion(.04))+
scale_color_manual(values=c("#0055ff","#8844ff","#ff55ff","gray50","#ff5555","black"))+
scale_shape_manual(values=c(16,16,16,16,16,4))+
scale_alpha_manual(values=c(1,1,1,1,1,0))+
scale_size_manual(values=c(.7,.7,.7,.7,.7,1))+
guides(color=guide_legend(ncol=2,byrow=F))+
theme(axis.text=element_text(size=7,color="black"),
  axis.text.x=element_text(angle=90,vjust=.5,hjust=1),
  axis.ticks=element_line(linewidth=.2,color="black"),
  axis.ticks.length=unit(3,"pt"),
  axis.ticks.length.x=unit(0,"pt"),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.justification="left",
  legend.key=element_blank(),
  legend.key.height=unit(8,"pt"),
  legend.key.width=unit(17,"pt"),
  legend.margin=margin(-2,,4),
  legend.position="top",
  legend.spacing.x=unit(2,"pt"),
  legend.spacing.y=unit(1,"pt"),
  legend.text=element_text(size=7,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.23),
  panel.spacing.x=unit(3,"pt"),
  panel.spacing.y=unit(3,"pt"),
  plot.margin=margin(5,5,5,5),
  plot.subtitle=element_text(size=6.9,margin=margin(,,4)),
  plot.title=element_text(size=7.4,face=2,margin=margin(1,,3)),
  strip.background=element_blank(),
  strip.text=element_blank())
ggsave("1.png",width=4,height=3,dpi=350*4)
system("magick 1.png -resize 25% 1.png")

library(data.table);library(ggplot2)

url="https://www2.census.gov/programs-surveys/popest/datasets/"
new=do.call(rbind,lapply(paste0(url,"2020-2023/national/asrh/nc-est2023-alldata-r-file0",1:8,".csv"),\(x)fread(x)))
old=do.call(rbind,lapply(paste0(url,"2010-2020/national/asrh/NC-EST2020-ALLDATA-R-File",sprintf("%02d",1:24),".csv"),\(x)fread(x)))

mult=new[YEAR==2020&MONTH==4.2,.(age=AGE,new=TOT_POP)]
mult=merge(mult,old[YEAR==2020&MONTH==4,.(age=AGE,old=TOT_POP)])[,.(mult=new/old,age)]

oldpop=old[YEAR!=2021&MONTH!=4.2&AGE!=999,.(year=YEAR,month=floor(MONTH),pop=TOT_POP,age=AGE)]
newpop=new[AGE!=999,.(year=YEAR,month=floor(MONTH),pop=TOT_POP,age=AGE)]

proj=fread("https://www2.census.gov/programs-surveys/popproj/datasets/2023/2023-popproj/np2023_d1_mid.csv")
proj=proj[SEX==0&ORIGIN==0&RACE==0&YEAR%in%2024:2025,.(x=c(7,19),age=rep(0:100,each=.N),pop=unlist(.SD[,-(1:5)]))]
proj=rbind(newpop[year==2023&month==12,.(x=0,pop,age)],proj)
newpop=rbind(newpop,proj[,spline(x,pop,xout=1:12),age][,.(age,year=2024,month=x,pop=y)])

oldmod=oldpop[!(year==2020&month>4)]
oldmod=merge(mult,oldmod[,.(year,month,pop,slope=seq(0,1,,.N)),age])
oldmod=oldmod[,.(age,year,month,pop=pop*(slope*mult+(1-slope)*1))]

a=cbind(rbind(oldpop[!(year==2020&month>=4)],newpop),group="Original unmodified population estimates")
a=rbind(a,cbind(rbind(oldmod[!(year==2020&month==4)],newpop),group="Population estimates based on 2010 census multiplied by linear slope"))

a=merge(a,fread("http://sars2.net/f/usdeadmonthly.csv")[!(year==2024&month>9)])

std=fread("https://www2.census.gov/programs-surveys/popest/datasets/2020-2023/national/asrh/nc-est2023-agesex-res.csv")
a=merge(std[SEX==0&AGE!=999,.(age=AGE,std=POPESTIMATE2020)][,std:=std/sum(std)],a)
a[,date:=as.Date(paste(year,month,1,sep="-"))]
a[,group:=factor(group,unique(group))]
a[,pop:=pop*lubridate::days_in_month(date)]

p=a[,.(asmr=sum(dead/pop*std*365e5)),.(date,group)]

p=merge(p[year(date)%in%2011:2019&month(date)%in%5:9,.(date=unique(p$date),base=predict(lm(asmr~date),.(date=unique(p$date)))),group],p)
p[,month:=month(date)]
p=merge(p[year(date)%in%2011:2019,.(monthly=mean(asmr-base)),month],p)

yearly=a[,.(pop=sum(pop),dead=sum(dead)),.(age,std,group,year=year(date))]
yearly=yearly[,.(asmr=sum(dead/pop*std*365e5)),.(year,group)]
yearly=merge(yearly[year%in%2011:2019,.(year=2011:2024,base=predict(lm(asmr~year),.(year=2011:2024))),group],yearly)

pct=yearly[,.(pct=(sum(asmr)/sum(base)-1)*100),.(year,group)][,x:=as.Date(paste0(year,"-7-1"))]

lab=c("Actual ASMR","2011-2019 linear baseline","Baseline adjusted for seasonality","Difference from seasonality-adjusted baseline")
p=p[,.(x=date,y=c(asmr,base,base+monthly,asmr-base-monthly),z=factor(rep(lab,each=.N),lab),group)]

pct=merge(p[,.(y=min(y)),group],pct)

xstart=as.Date("2017-1-1");xend=as.Date("2025-1-1")
xbreak=seq(xstart,xend,"6 month");xlab=ifelse(month(xbreak)==7,year(xbreak),"")

ggplot(p,aes(x=x+14,y))+
facet_wrap(~group,ncol=1,dir="v",scales="free")+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray85",linewidth=.25)+
geom_vline(xintercept=as.Date("2020-1-1"),linetype="44",linewidth=.3)+
geom_line(data=p[z!=lab[4]],aes(linetype=z,color=z),linewidth=.3)+
geom_line(data=p[z==z[1]],aes(linetype=z,color=z),linewidth=.3)+
geom_point(aes(alpha=z),stroke=0,size=.6)+
geom_rect(data=p[z==lab[4]],aes(xmin=x,xmax=x%m+%months(1),ymin=pmin(0,y),ymax=pmax(0,y),fill=z),linewidth=.1,color="gray40")+
geom_text(data=pct,aes(x,y=y,label=paste0(ifelse(pct>0,"+",""),sprintf("%.1f",pct),"%")),vjust=.8,size=2.4)+
labs(title="Age-standardized mortality rate per 100,000 person-years in United States",x=NULL,y=NULL)+
scale_x_continuous(expand=c(0,0),limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(expand=c(.06,0,.03,0),breaks=\(x)pretty(x,6))+
scale_color_manual(values=c("black","gray60","gray60","gray40"))+
scale_fill_manual(values=rep("gray60",4))+
scale_linetype_manual(values=c("solid","42","solid","solid"))+
scale_alpha_manual(values=c(1,0,0,0),guide="none")+
guides(fill=guide_legend(order=3,keyheight=unit(1,"pt")),color=guide_legend(order=1),linetype=guide_legend(order=1))+
coord_cartesian(clip="off")+
theme(axis.text=element_text(size=7,color="black"),
  axis.ticks.x=element_line(color=alpha("black",c(1,0))),
  axis.ticks=element_line(linewidth=.3,color="black"),
  axis.ticks.length=unit(3,"pt"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.title=element_text(size=8),
  legend.background=element_blank(),
  legend.key=element_blank(),
  legend.key.height=unit(9,"pt"),
  legend.key.width=unit(19,"pt"),
  legend.position="top",
  legend.justification="left",
  legend.spacing.x=unit(2,"pt"),
  legend.box.margin=margin(),
  legend.spacing.y=unit(0,"pt"),
  legend.margin=margin(,,2),
  legend.text=element_text(size=7,vjust=.5),
  legend.title=element_blank(),
  legend.direction="vertical",
  legend.box="horizontal",
  legend.box.spacing=unit(0,"pt"),
  legend.spacing=unit(0,"pt"),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.3),
  panel.spacing.y=unit(2,"pt"),
  plot.margin=margin(4,4,2,4),
  plot.subtitle=element_text(size=6.5,margin=margin(,,4)),
  plot.title=element_text(size=7.4,face=2,margin=margin(1,,3)),
  strip.background=element_blank(),
  strip.text=element_text(size=7,face=2,margin=margin(,,2)))
ggsave("1.png",width=4.7,height=3.9,dpi=400*4)

sub="\u00a0     The data for deaths was retrieved from CDC WONDER on November 21st 2024 UTC, so some deaths are still missing in 2024 because of a registration delay.
      The yearly total excess mortality percentage is shown above the year.
      ASMR was calculated by single year of age so that the vintage 2023 mid-2020 resident population estimates were used as the standard population.
      The baseline was calculated by doing a linear regression for monthly data in 2011 to 2019. Only May to September were included in the linear regression. In order to calculate the seasonality-adjusted baseline, the average difference between the actual ASMR and baseline during each January of 2011 to 2019 was added to the baseline for each January, and similarly for other months.
      Monthly resident population estimates by single year of age were downloaded from www2.census.
gov/programs-surveys/popest/datasets/2020-2023/national/asrh/nc-est2023-alldata-r-file0{1..8}.csv, and from the corresponding files in the 2010-2020 directory.
      Vintage 2023 population estimates based on the 2020 census are used for April 2020 to December 2023, and vintage 2020 population estimates based on the 2010 census are used for January 2011 to March 2020. In the bottom panel the 2010-based estimates were multiplied by a linear slope so that they were equal to the 2020-based estimates in April 2020 when the two sets of estimates were merged. The 2010-based estimates cover a total of 121 months from April 2010 until April 2020, and for example for the age 80 the ratio between the 2020-based and 2010-based estimates was about 0.946 in April 2020, so the population estimate for age 80 was multiplied by about 1/121*0.946+120/121 in May 2010, by about 2/121*0.946+119/121 in June 2010, and so on. A similar method is used in the Census Bureau's intercensal population estimates, but intercensal estimates for 2010-2020 have not yet been published.
      Monthly population estimates for 2024 were calculated as a spline interpolation of the population estimate on December 2023 and mid-year population projections for 2024 and 2025: www2.census.gov/programs-surveys/popproj/datasets/2023/2023-popproj/np2023_d1_mid.csv."

system(paste0("mar=120;w=`identify -format %w 1.png`;magick 1.png \\( -size $[w-mar*2]x -font Arial -interline-spacing -4 -pointsize 152 caption:'",gsub("'","'\\\\''",sub),"' -gravity southwest -splice $[mar]x70 \\) -append -resize 25% -colorspace gray -colors 16 1.png"))

library(data.table);library(ggplot2)

t=fread("http://sars2.net/f/uspopdeadmonthly.csv")[year%in%2011:2024]
t[,date:=as.integer(factor(date))-84][,cmr:=dead/persondays]
t=merge(t,t[year<2020,.(date=unique(t$date),base=predict(lm(cmr~date),.(date=unique(t$date)))),age])
t=merge(t[year<2020,.(monthly=mean(cmr/base)),.(month,age)],t)
rate=t[date>0,.(base=base*365/12,mult=ifelse(date%in%25:48,cmr/(base*monthly),1)),,.(date,age)]

set.seed(0)
pop=fread("http://sars2.net/f/uspopdead.csv")[year==2020]
people=pop[age>0,rep(age,pop/10)];people=people2=people+runif(length(people))
sim=data.table(month=1:84,dead=0,dead2=0)

set.seed(0)
for(i in 1:nrow(sim)){
  died=runif(length(people))<rate[date==i,base*mult][people]
  people=pmin(100,people[!died]+1/12)
  sim$dead[i]=sum(died)
  print(i)
}

set.seed(0)
for(i in 1:nrow(sim)){
  died=runif(length(people2))<rate[date==i,base][people2]
  people2=pmin(100,people2[!died]+1/12)
  sim$dead2[i]=sum(died)
  print(i)
}

sim$base=sim[month<25,predict(lm(dead~month),sim)]
lab=c("With elevated mortality in 2020-2021","No elevated mortality in 2020-2021","2018-2019 linear trend")
p=melt(sim,id=1)[,.(x=month,y=value,z=`levels<-`(variable,lab))]

sum=sim[month%in%25:48,colSums(.SD)];sum2=sim[month>48,colSums(.SD)]
note=stringr::str_wrap(sprintf("The scenario with elevated mortality in 2020-2021 has %.1f%% more deaths in 2020-2021 and %.1f%% less deaths in 2022-2024.",(sum[2]/sum[3]-1)*100,(sum2[3]/sum2[2]-1)*100),32)

xstart=1;xend=nrow(sim)+1;xbreak=seq(xstart,xend,6);xlab=c(rbind("",2018:2024),"")
ybreak=pretty(c(0,p$y),8);ystart=0;yend=max(p$y)+.02

ggplot(p,aes(x+.5,y))+
geom_vline(xintercept=seq(xstart,xend,12),linewidth=.3,color="gray83",lineend="square")+
geom_hline(yintercept=c(ystart,yend),linewidth=.3,lineend="square")+
geom_vline(xintercept=c(xstart,xend),linewidth=.3,lineend="square")+
geom_line(aes(color=z,linetype=z),linewidth=.35)+
geom_point(aes(color=z,alpha=z),stroke=0,size=.7)+
labs(title="Simulation for deaths per month in one tenth of US population",x=NULL,y=NULL)+
annotate(geom="label",x=42,y=yend*.48,label=note,hjust=0,label.r=unit(2,"pt"),label.padding=unit(3,"pt"),label.size=.2,size=2.4,lineheight=.9)+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak)+
scale_color_manual(values=c("#ff6666","black","gray60"))+
scale_alpha_manual(values=c(1,1,0))+
scale_linetype_manual(values=c("solid","solid","42"))+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=7,color="black"),
  axis.ticks=element_line(linewidth=.3),
  axis.ticks.length=unit(3,"pt"),
  axis.ticks.length.x=unit(0,"pt"),
  legend.direction="vertical",
  legend.justification=c(0,0),
  legend.key=element_blank(),
  legend.key.height=unit(9,"pt"),
  legend.key.width=unit(20,"pt"),
  legend.background=element_rect(color="black",linewidth=.3),
  legend.margin=margin(3,3,3,3),
  legend.position=c(0,0),
  legend.spacing.x=unit(1.5,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=7),
  legend.title=element_blank(),
  panel.background=element_blank(),
  plot.margin=margin(4,4,2,4),
  plot.title=element_text(size=7.5,face="bold",margin=margin(1,,3)))
ggsave("1.png",width=3.7,height=2,dpi=380*4)
sub=paste0("The initial population for each age was one tenth of the mid-2018 US resident population estimates except infants of age zero were excluded, so the total initial population size was about 32 million people. The projection of a linear trend for CMR in 2011-2019 was used for each age throughout the simulation, except in the red scenario the mortality rate for each combination of month and age in 2020-2021 was multiplied by the ratio between the actual CMR for the age during the month and a seasonality-adjusted trend for CMR during the same month.")
sub=paste0(sub," The length of each month is 1/12 years. Births are not simulated. Seasonal variation in mortality or differences in the health of people were not simulated. The same seed was used for both scenarios so the scenarios are identical in 2018 and 2019.")
system(paste0("mar=120;w=`identify -format %w 1.png`;magick 1.png \\( -size $[w-mar*2]x -font Arial -interline-spacing -4 -pointsize 146 caption:'",gsub("'","'\\'",sub),"' -gravity southwest -splice $[mar]x70 \\) -append -resize 25% -colors 256 1.png"))

library(data.table);library(ggplot2)

months=12*60;minage=0;maxage=Inf

covid=do.call(rbind,lapply(2020:2022,\(i)fread(paste0(i,".csv.gz"))[ucod=="U071"&restatus!=4,.(age,month=(year-2020)*12+monthdth)]))
covid=covid[age!=9999][,age:=ifelse(age<2000,age%%1000,0)][age>=minage&age<=maxage,.(pop=.N),.(month,age=age*12+6)]

pop=fread("http://sars2.net/f/uspopdead.csv")[year==2019]
risk=pop[age<100,.(age,risk=1-(1-dead/pop)^(1/12))]
risk=c(risk$risk,risk[age>=90,exp(predict(lm(log(risk)~age),.(age=100:120)))])
risk=spline(seq(6,,12,121),risk,xout=0:(121*12-1))$y

all=fread("http://sars2.net/f/uspopdeadmonthly.csv")[year%in%2011:2022]
all[,date:=(year-2020)*12+month][,cmr:=dead/persondays]
all=merge(all,all[year<2020,.(date=unique(all$date),base=predict(lm(cmr~date),.(date=unique(all$date)))),age])
all=merge(all[year<2020,.(monthly=mean(cmr/base)),.(month,age)],all)
allcause=all[age>=minage&age<=maxage&date%in%1:36,.(pop=sum(dead)-sum(base*persondays*monthly)),,.(age=age*12+6,month=date)]

dead0=double();cur=pop[age>=minage&age<=maxage,.(age=age*12+6,pop,month=1)]
for(i in 1:months){
  died=cur[,risk[pmin(age+1,length(risk))]*pop]
  dead0[i]=sum(died)
  cur=cur[,.(month=i+1,age=age+1,pop=pop-died)][,.(pop=sum(pop)),.(month,age)]
}

dead1=double();prev=NULL
for(i in 1:months){
  cur=rbind(prev,covid[month==i])[,.(pop=sum(pop)),.(month,age)]
  died=cur[,risk[pmin(age+1,length(risk))]*pop]
  dead1[i]=sum(died)
  prev=cur[,.(month=i+1,age=age+1,pop=pop-died)]
}

dead2=double();prev=NULL
for(i in 1:months){
  cur=rbind(prev,allcause[month==i])[,.(pop=sum(pop)),.(month,age)]
  died=cur[,risk[pmin(age+1,length(risk))]*pop]
  dead2[i]=sum(died)
  prev=cur[,.(month=i+1,age=age+1,pop=pop-died)]
}

pfe=c(-3.5,-4.9,-5.9,-5.9,-4.7,-2.9,-1.7,-1)
pfe=spline(seq(1,351,50),pfe,xout=seq(1,351,1/7))
pfe=setDT(pfe)[,.(y=mean(y)),.(x=(x*7)%/%(365.25/12)+15)]

lab="My model based on UCD COVID deaths in 2020-2022"
lab=c(lab,"My model based on all-cause excess deaths in 2020-2022")
lab=c(lab,"Approximation of Ethical Skeptic's PFE arrival function")

p=data.table(x=1:months,y=-c(dead1,dead2)/dead0*100)
p[,z:=factor(rep(lab[1:2],each=.N/2),lab[1:2])]
p=rbind(p,pfe[,z:=lab[3]])

xstart=1;xend=months+1;xbreak=seq(xstart,xend,12*2);xlab=seq(2020,,2,length(xbreak))
ystart=floor(min(p$y));yend=0;ybreak=ystart:yend

ggplot(p,aes(x,y))+
geom_hline(yintercept=c(ystart,yend),linewidth=.3,lineend="square")+
geom_vline(xintercept=c(xstart,xend),linewidth=.3,lineend="square")+
geom_line(aes(color=z),linewidth=.35)+
labs(title="PFE-adjusted baseline models",x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak,labels=\(x)paste0(x,"%"))+
coord_cartesian(ylim=c(ystart,yend),clip="off",expand=F)+
scale_color_manual(values=c("#f6f","#f66","black"))+
theme(axis.text=element_text(size=7,color="black"),
  axis.text.x=element_text(angle=90,vjust=.5,hjust=1),
  axis.text.y=element_text(margin=margin(,1)),
  axis.ticks=element_line(linewidth=.3),
  axis.ticks.length.x=unit(3,"pt"),
  legend.direction="vertical",
  legend.justification=c(1,.5),
  legend.key=element_blank(),
  legend.key.height=unit(9,"pt"),
  legend.key.width=unit(17,"pt"),
  legend.background=element_rect(color="black",linewidth=.3),
  legend.margin=margin(2,2,2,2),
  legend.position=c(1,.5),
  legend.spacing.x=unit(1.5,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=6.8),
  legend.title=element_blank(),
  panel.background=element_blank(),
  plot.margin=margin(4,7,2,4),
  plot.title=element_text(size=7.3,face="bold",margin=margin(1,,3)))
ggsave("1.png",width=3.8,height=1.8,dpi=380*4)
system("magick 1.png -resize 25% -colors 256 1.png")

library(data.table)

kimi=\(x){na=is.na(x);x[na]=0;e=floor(log10(ifelse(x==0,1,abs(x))));e2=pmax(e,0)%/%3+1;x[]=ifelse(abs(x)<1e3,round(x),paste0(sprintf(paste0("%.",ifelse(e%%3==0,1,0),"f"),x/1e3^(e2-1)),c("","k","M","B","T")[e2]));x[na]="NA";x}

t=fread("http://sars2.net/f/wonderbyicd.csv.gz")[type=="MCD"&year!=2024&code%like%"D6"]
t=t[code%in%t[,sum(dead),code][V1>=1000,code]]
t=t[order(code)]
t=t[,label:=paste0(code,": ",sub(" \\[.*","",cause))][,label:=factor(label,unique(label))]

m=t[,xtabs(dead~label+year)]
maxcolor=max(m);exp=.6
pal=hsv(c(21,21,12,6,3,0,0)/36,c(0,.5,.5,.5,.5,.5,0),rep(1:0,c(6,1)))

pheatmap::pheatmap(m^exp,filename="1.png",
  display_numbers=ifelse(m<10,m,kimi(m)),
  cluster_rows=F,cluster_cols=F,legend=F,cellwidth=17,cellheight=15,fontsize=9,fontsize_number=8,
  border_color=NA,na_col="gray90",
  number_color=ifelse(abs(m^exp)>maxcolor^exp*.8,"white","black"),
  breaks=seq(0,maxcolor^exp,,256),
  colorRampPalette(pal)(256))

system("w=`identify -format %w 1.png`;pad=20;magick -pointsize 44 -font Arial-Bold \\( -size $[w-pad] caption:'CDC WONDER: Yearly deaths by underlying cause of death (ICD codes starting with D6 with at least 1,000 total deaths)' -splice $[pad]x24 \\) 1.png -append 1.png")

library(data.table);library(ggplot2)

t=fread("http://sars2.net/f/wonderthrombocytopenia.csv")
t=t[,.(x=as.Date(paste(year,month,1,sep="-")),y=dead,z=factor(cause,unique(cause)),type)]
t[,y:=y/lubridate::days_in_month(x)]

p=rbind(t[type=="mcd"][,type:=1],dcast(t,z+x~type,value.var="y")[,.(z,x,type=2,y=mcd-nafill(and,,0))])

p=rbind(p,p[type==1&year(x)%in%2014:2019,{i=p[year(x)>=2014,x];.(x=i,y=predict(lm(y~x),.(x=i)),type=3)},z])
p[,type:=factor(type,,c("MCD COVID not subtracted","MCD COVID subtracted","2014-2019 linear trend"))]

xstart=as.Date("1999-1-1");xend=as.Date("2025-1-1")
xbreak=seq(xstart,xend,"6 month");xlab=ifelse(month(xbreak)==7,year(xbreak),"")
ybreak=pretty(p$y,6);ystart=0;yend=max(p$y)

lab=p[rowid(z)==1][order(z),.(z,x=xend-150,y=c(32.3,9.7,8,11.4,34))]

ggplot(p,aes(x+15,y))+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray90",linewidth=.25)+
geom_hline(yintercept=c(ystart,yend),linewidth=.25,lineend="square")+
geom_vline(xintercept=c(xstart,xend),linewidth=.25,lineend="square")+
geom_line(aes(y=ifelse(y<0,NA,y),color=z,alpha=type,linetype=type),linewidth=.3)+
geom_label(data=lab,aes(label=z,color=z,x=x),size=2.4,hjust=1,fill=alpha("white",.85),label.r=unit(0,"pt"),label.padding=unit(1,"pt"),label.size=0,show.legend=F)+
labs(title="CDC WONDER: Monthly MCD deaths divided by number of days in month",,x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak)+
scale_color_manual(guide="none",values=c(hsv(0,1,.6),hsv(30/36,.6,.6),hsv(12/36,1,.7),"black",hsv(22/36,1,.7)))+
scale_alpha_manual(values=c(1,.4,1))+
scale_linetype_manual(values=c("solid","solid","42"))+
coord_cartesian(clip="off",expand=F)+
guides(linetype=guide_legend(order=2),alpha=guide_legend(order=2))+
theme(axis.text=element_text(size=7,color="black"),
  axis.text.x=element_text(angle=90,vjust=.5,hjust=1),
  axis.text.y=element_text(margin=margin(,1)),
  axis.ticks=element_line(linewidth=.25,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(3,"pt"),
  legend.background=element_rect(fill=alpha("white",1),color="black",linewidth=.25),
  legend.box="vertical",
  legend.box.just="left",
  legend.box.spacing=unit(0,"pt"),
  legend.justification=c(0,1),
  legend.key=element_blank(),
  legend.margin=margin(2,3,2,2),
  legend.key.height=unit(8,"pt"),
  legend.key.width=unit(17,"pt"),
  legend.position=c(0,1),
  legend.spacing.x=unit(1,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=7,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  plot.margin=margin(4,4,4,4),
  plot.title=element_text(size=7,face=2,margin=margin(,,3)))
ggsave("1.png",width=4.2,height=2.8,dpi=380*4)
system("magick 1.png -resize 25% -colors 256 1.png")

t=fread("http://sars2.net/f/wonderthrombostate.csv")

disp=t[,tapply(dead,.(state,year),c)]
disp=disp[rowSums(is.na(disp))<2,]
m=disp/apply(disp,1,max,na.rm=T)

pal=sapply(255:0,\(i)rgb(i,i,i,maxColorValue=255))

step=ceiling(nrow(m)/2)
for(i in 1:2){
  rows=((i-1)*step):min(nrow(m),i*step)
  pheatmap::pheatmap(m[rows,],filename=paste0("i",i,".png"),display_numbers=disp[rows,],
    cluster_rows=F,cluster_cols=F,legend=F,cellwidth=15,cellheight=15,fontsize=9,fontsize_number=8,
    border_color=NA,na_col="white",
    number_color=ifelse(!is.na(m[rows,])&m[rows,]>.4,"white","black"),
    breaks=seq(0,1,,256),
    colorRampPalette(pal)(256))
}

system("magick i[12].png +append 1.png;w=`identify -format %w 1.png`;pad=26;magick -pointsize 42 -font Arial \\( -size $[w-pad*2] caption:'CDC WONDER: Yearly deaths by state with MCD D69.5 (secondary thrombocytopenia). Each row has its own color scale where the maximum value is shown in black. States with data missing for more than one year are omitted. The data was retrieved on December 2nd 2024 UTC so some deaths in 2024 are still missing.' -splice $[pad]x20 \\) 1.png -append 1.png")

library(data.table);library(ggplot2)

t=fread("http://sars2.net/f/wondernaturalsingle.csv")
t=t[age%in%2:17&date!="2024-11",.(dead=sum(dead)),date]

d=t[,.(x=as.Date(paste0(date,"-1")),dead)]
d[,dead:=dead/lubridate::days_in_month(x)]
d$base=d[year(x)%in%2013:2019,predict(lm(dead~x),d)]
d$base=d$base+d[year(x)%in%2013:2019,mean(dead-base),.(month=month(x))]$V1[month(d$x)]
p=melt(d,id=1,,"z","y")
p[,z:=factor(z,unique(z),c("Actual deaths","2013-2019 linear baseline"))]

xstart=as.Date("2015-1-1");xend=as.Date("2025-1-1")
p=p[x>=xstart&x<=xend]
xbreak=seq(xstart+182,xend,"year");xlab=year(xbreak)
ybreak=pretty(p$y,7);ystart=ybreak[1];yend=max(ybreak)

ggplot(p,aes(x=x+14,y))+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray83",linewidth=.25)+
geom_line(aes(color=z),linewidth=.3)+
geom_point(aes(alpha=z,color=z),stroke=0,size=.6)+
labs(x=NULL,y=NULL,title="CDC WONDER: Monthly deaths in ages 2-17 with UCD A-Q (natural causes\nwithout R codes, U codes, or COVID), divided by number of days in month",subtitle="The data was retrieved on December 5th 2024 UTC, so many deaths in 2024 are
still missing because of a registration delay.")+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak)+
scale_color_manual(values=c("black","blue"))+
scale_alpha_manual(values=c(1,0))+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=7,color="black"),
  axis.text.y=element_text(margin=margin(,1.5)),
  axis.ticks=element_line(linewidth=.25,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(2,"pt"),
  legend.background=element_blank(),
  legend.box.background=element_rect(color="black",linewidth=.25),
  legend.justification=c(0,1),
  legend.key=element_blank(),
  legend.key.height=unit(10,"pt"),
  legend.key.width=unit(17,"pt"),
  legend.position=c(0,1),
  legend.spacing.x=unit(1,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.margin=margin(3,4,3,3),
  legend.text=element_text(size=7,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.3),
  plot.margin=margin(4,4,3,4),
  plot.subtitle=element_text(size=7,margin=margin(,,3)),
  plot.title=element_text(size=7.2,face=2,margin=margin(1,,3)))
ggsave("1.png",width=3.9,height=2.4,dpi=400*4)
system("magick 1.png -resize 25% PNG8:1.png")

library(data.table);library(ggplot2)

ma=\(x,b=1,f=b){x[]=rowMeans(embed(c(rep(NA,b),x,rep(NA,f)),f+b+1),na.rm=T);x}

t=fread("http://sars2.net/f/wonder217naturalweekly.csv")[date<="2024-11-1"]
t[,dead:=ma(dead,3,2)]

t=merge(t,t[year<2020,.(base=mean(dead),base2=mean(dead)),week])
slope=t[year<2020,mean(dead),year][,predict(lm(V1~year),.(year=2018:2024))]
slope2=t[year<2020&week%in%15:35,mean(dead),year][,predict(lm(V1~year),.(year=2018:2024))]
t[,base:=base*(slope/mean(slope[1:2]))[factor(t$year)]]
t[,base2:=base2*(slope2/mean(slope2[1:2]))[factor(t$year)]]

t$base3=t[year<2020,predict(lm(dead~date),t)]
t$base3=t$base3+t[year<2020,mean(dead-base3),week]$V1[t$week]

lab=c("Actual deaths","Baseline with slope determined by total on all weeks","Baseline with slope determined by total on weeks 15 to 35","Linear regression of weekly data with week number residuals added")
p=t[,.(x=date,y=c(dead,base,base2,base3),z=factor(rep(lab,each=.N),lab))]

p$facet="Deaths"
p=rbind(p,merge(p[z!=z[1]],p[z==z[1],.(x,actual=y)])[,.(x,y=actual-y,z,facet="Excess deaths")])
p[,facet:=factor(facet,unique(facet))]

xstart=as.Date("2018-1-1");xend=as.Date("2025-1-1")
xbreak=seq(xstart,xend,"6 month");xlab=ifelse(month(xbreak)==7,year(xbreak),"")

ggplot(p,aes(x=x,y))+
facet_wrap(~facet,dir="v",scales="free")+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray83",linewidth=.4)+
geom_segment(data=tail(p,1),x=xstart,xend=xend,y=0,yend=0,linewidth=.4)+
geom_line(aes(color=z),linewidth=.5)+
geom_text(data=p[,.(y=max(y)),facet],x=pmean(xstart,xend),aes(label=facet),size=3.8,fontface=2,vjust=1.4)+
labs(x=NULL,y=NULL,title="CDC WONDER: Weekly deaths in ages 2-17 with UCD A-Q (natural causes\nwithout R codes, U codes, or COVID)",subtitle="The data was retrieved on January 3rd 2025 UTC, so many deaths in 2024 are
still missing because of a registration delay.")+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab,expand=expansion(0,0))+
scale_y_continuous(breaks=\(x)pretty(x,7),expand=expansion(.03,0))+
scale_color_manual(values=c("black","blue","#8888ff","#ff6666","#ff66ff"))+
coord_cartesian(clip="off")+
theme(axis.text=element_text(size=11,color="black"),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(5,"pt"),
  legend.background=element_blank(),
  legend.box.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.direction="vertical",
  legend.justification="left",
  legend.key=element_blank(),
  legend.key.height=unit(11,"pt"),
  legend.key.width=unit(26,"pt"),
  legend.margin=margin(,,6,-2),
  legend.position="top",
  legend.spacing.x=unit(3,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.4),
  panel.spacing=unit(2,"pt"),
  plot.margin=margin(6,6,1,5),
  plot.title=element_text(size=11.5,face=2,margin=margin(2,,4)),
  strip.background=element_blank(),
  strip.text=element_blank())
ggsave("1.png",width=6.3,height=4.8,dpi=300*4)

sub="\u00a0    The deaths are plotted as a moving average where the window extends 3 weeks backwards and 2 weeks forwards. The baselines were also calculated based on the moving average of deaths.
     In the bright blue baseline the average of weekly deaths was calculated in 2018 and 2019, a linear regression of the averages was projected to 2018-2024, and multipliers for each year were calculated by dividing the value of the projection during the year with the average value of the projection in the years 2018 and 2019. And then the baseline for each week 1 was calculated by multiplying the average deaths on weeks 1 of 2018 and 2019 with the yearly multiplier, and similarly for other week numbers. The light blue baseline was calculated in the same way except the yearly averages were calculated based on weeks 15 to 35 only.
     The red line was calculated by first doing a linear regression of weekly data in 2018-2019, and then the average difference between the actual deaths and the linear regression on weeks 1 of 2018 and 2019 was added to the baseline for each week 1, and similarly for other week numbers.."

system(paste0("f=1.png;mar=100;w=`identify -format %w $f`;magick \\( $f \\) \\( -size $[w-mar*2]x -font Arial -interline-spacing -3 -pointsize $[42*4] caption:'",gsub("'","'\\\\'",sub),"' -gravity southwest -splice $[mar]x80 \\) -append -resize 25% -colors 256 1.png"))

icd=fread("http://sars2.net/f/wonderbyicd.csv.gz")[rowid(code)==1,.(cause,icd=sub("\\.","",code))]

# sars2.net/stat.html#Download_fixed_width_and_CSV_files_for_the_NVSS_data_used_at_CDC_WONDER
t=do.call(rbind,lapply(2013:2022,\(i)fread(paste0(i,".csv.gz"))[age%in%1002:1017,.(ucd=ucod,year)]))

a=t[,.(dead=.N),.(icd=ucd,year)]
a=merge(do.call(CJ,lapply(a[,-3],unique)),a,all.x=T)[is.na(dead),dead:=0]

a=a[a[year%in%2013:2019,.(year=2013:2022,base=predict(lm(dead~year),.(year=2013:2022))),icd],on=names(a)[1:2]]
a=merge(icd,a)

m=a[icd%like%"[A-Q]",tapply(dead-base,.(paste0(icd,": ",cause),year),c)]
m=m[order(-rowSums(m[,ncol(m)-2:0])),][1:40,]

maxcolor=max(abs(m))
exp=.8
pal=hsv(rep(c(21,0)/36,5:4),c(1,.8,.6,.3,0,.3,.6,.8,1),c(.3,.65,1,1,1,1,1,.65,.3))

pheatmap::pheatmap(sign(m)*abs(m)^exp,filename="1.png",display_numbers=kimi(m),
  cluster_rows=F,cluster_cols=F,legend=F,cellwidth=14,cellheight=14,fontsize=9,fontsize_number=8,
  border_color=NA,na_col="gray90",
  number_color=ifelse(!is.na(m)&abs(m)^exp>.7*(maxcolor^exp),"white","black"),
  breaks=seq(-maxcolor^exp,maxcolor^exp,,256),
  colorRampPalette(pal)(256))
system("convert 1.png -extent 2000x 1.png;w=`identify -format %w 1.png`;pad=26;magick -pointsize 44 -font Arial \\( -size $[w-pad] caption:'NVSS, ages 2 to 17: Excess number of deaths by underlying cause of death relative to 2015-2019 linear trend, top 40 causes shown sorted by 2020-2022 total' -splice $[pad]x24 \\) 1.png -append 1.png")

library(data.table);library(ggplot2)

order="RI MA VT ME CT NJ NM NY MD HI VA PA NC NH DE WA CA SD CO FL OR MN IL AZ NV TX KS UT WI OK NE AK SC IA AR MI ND MO KY GA MT WV OH AL TN IN ID LA MS WY"
order=state.name[match(strsplit(order," ")[[1]],state.abb)]

t=fread("http://sars2.net/f/wonderstatemalignant.csv")[!cause%like%"COVID"&age=="0-64"&date!="2024-11"]
g=c("15 most vaccinated states","Middle 20 states","15 least vaccinated states")
t$group=g[2]
t[state%in%order[1:15],group:=g[1]]
t[state%in%tail(order,15),group:=g[3]]
d=t[,.(dead=sum(dead)),,.(group=factor(group,g),date=as.Date(paste0(date,"-1")))]

d[,dead:=dead/lubridate::days_in_month(date)]
d[,month:=month(date)]
d$base=d[year(date)%in%2014:2019&month%in%4:9,predict(lm(dead~date),.(date=unique(d$date))),group]$V1
d$base2=d[year(date)%in%2018:2019&month%in%4:9,predict(lm(dead~date),.(date=unique(d$date))),group]$V1
d=merge(d[year(date)%in%2014:2019,.(monthly=mean(dead-base)),.(group,month)],d)
d=merge(d[year(date)%in%2018:2019,.(monthly2=mean(dead-base2)),.(group,month)],d)
lab=c("Actual deaths","2014-2019 linear trend","2018-2019 linear trend")
p=d[,.(x=date,y=c(dead,base+monthly,base2+monthly2),z=factor(rep(lab,each=.N),lab),group)]

p=merge(p[year(x)==2019,.(base=mean(y)),group],p)[,y:=y/base*100]

xstart=as.Date("2010-1-1");xend=as.Date("2025-1-1")
p=p[x%in%xstart:xend]
xbreak=seq(xstart,xend,"6 month");xlab=ifelse(month(xbreak)==7,year(xbreak),"")
ybreak=pretty(p$y,6);ystart=0;yend=max(p$y)
yend=124.85733
ggplot(p,aes(x+15,y))+
facet_wrap(~group,ncol=1,strip.position="top")+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray90",linewidth=.25)+
geom_hline(yintercept=100,linewidth=.25,linetype="42",color="gray50")+
geom_line(aes(y=ifelse(y<0,NA,y),color=z),linewidth=.3)+
geom_text(data=p[rowid(group)==1],aes(label=group),y=yend,x=pmean(xstart,xend),vjust=1.7,size=2.4,fontface=2)+
labs(title="CDC WONDER, ages 0-64: Monthly deaths with MCD malignant neoplasms divided\nby number of days in month, shown as percentage of average value in year 2019",subtitle=stringr::str_wrap("The baseline was calcualted by first doing a linear regression for deaths in April to September of each year during the fitting period, then the regression was projected over the entire period, and then the average difference between actual deaths and the baseline during each January of the fitting period was added to the baseline for each January, and similarly for other months. The baseline was not adjusted for population size or age. The Statista dataset 1202065 was used to rank states by the percentage of vaccinated people as of April 2023.",88),x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(85.54061,124.85733),breaks=pretty,labels=\(x)paste0(x,"%"))+
scale_color_manual(values=c("black","blue","#aaf"))+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=7,color="black"),
  axis.text.y=element_text(margin=margin(,1)),
  axis.ticks=element_line(linewidth=.25,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(3,"pt"),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.justification="right",
  legend.key=element_blank(),
  legend.key.height=unit(7,"pt"),
  legend.key.width=unit(19,"pt"),
  legend.margin=margin(,,4),
  legend.position="top",
  legend.spacing.x=unit(1,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=7,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.25),
  panel.spacing=unit(0,"pt"),
  plot.margin=margin(3,3,3,3),
  plot.subtitle=element_text(size=6.8,margin=margin(,,2)),
  plot.title=element_text(size=7,face=2,margin=margin(1,,3)),
  strip.background=element_blank(),
  strip.text=element_blank())
ggsave("1.png",width=4.3,height=3.8,dpi=380*4)
system("magick 1.png -resize 25% -colors 256 1.png")

library(data.table);library(ggplot2)

# http://sars2.net/stat.html#Download_fixed_width_and_CSV_files_for_the_NVSS_data_used_at_CDC_WONDER
t=do.call(rbind,lapply(2010:2022,\(i)fread(paste0(i,".csv.gz"))[restatus!=4]))

a=t[rowSums(t=="C419",na.rm=T)>0]
a=a[year>=2011&age<65]
a=a[,.(dead=.N),.(age=ifelse(age%%1000%in%c(1,9),age%%1000,0),year,month=monthdth)]

pop=fread("http://sars2.net/f/uspopdeadmonthly.csv")
a=merge(pop[,.(age,year,month,pop=persondays)],a)
a=merge(pop[year==2020&month==4,.(age,std=pop/sum(pop))],a)
a[,x:=as.Date(paste(year,month,1,sep="-"))]

d=a[,.(y=sum(dead),group="Raw deaths"),x]
d=rbind(d,a[,.(y=sum(dead/pop*std*365e5),group="ASMR per 100,000 person-years"),x])
d[,group:=factor(group,unique(group))]

d$base=d[year(x)<2020,predict(lm(y~x),d[rowid(x)==1]),group]$V1
d[,month:=month(x)]
d=merge(d[year(x)<2020,.(monthly=mean(y-base)),.(group,month)],d)

lab=c("Actual mortality","2011-2019 baseline","Excess deaths")
p=d[,.(x=rep(x,3),y=c(y,base+monthly,y-base-monthly),z=factor(rep(lab,each=.N),lab),group)]

xstart=as.Date("2011-1-1");xend=as.Date("2023-1-1")
xbreak=seq(xstart+182,xend,"year");xlab=year(xbreak)
ybreak=pretty(p$y,7);ystart=ybreak[1];yend=max(ybreak)

pct=d[,.(pct=(sum(y)/sum(base+monthly)-1)*100),.(group,year=year(x))]
pct[,x:=as.Date(paste0(year,"-7-1"))]
pct=merge(p[,.(y=min(y)),group],pct)

ggplot(p,aes(x=x+14,y))+
facet_wrap(~group,ncol=1,dir="v",scales="free")+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray85",linewidth=.25)+
geom_vline(xintercept=as.Date("2020-1-1"),linetype="44",linewidth=.3)+
geom_line(data=p[z!=lab[3]],aes(linetype=z,color=z),linewidth=.3)+
geom_line(data=p[z==z[1]],aes(linetype=z,color=z),linewidth=.3)+
geom_point(aes(alpha=z),stroke=0,size=.6)+
geom_rect(data=p[z==lab[3]],aes(xmin=x,xmax=x%m+%months(1),ymin=pmin(0,y),ymax=pmax(0,y),fill=z),linewidth=.1,color="gray40")+
geom_text(data=pct,aes(x,y=y,label=paste0(ifelse(pct>0,"+",""),sprintf("%.1f",pct),"%")),vjust=.7,size=2.4,color="gray40")+
labs(title="NVSS: Deaths in ages 0-64 with MCD C41.9 (Bone and articular cartilage,
unspecified - Malignant neoplasms)",subtitle="The percentage above the year shows the yearly total excess mortality. ASMR was calculated by single year of age so that the vintage 2023 resident population estimates for April 2020 were used as the standard population. The 2010-based population estimates were multiplied by a linear slope to get rid of a sudden jump when they were merged with the 2020-based population estimates."|>stringr::str_wrap(92),x=NULL,y=NULL)+
scale_x_continuous(expand=c(0,0),limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(expand=c(.06,0,.03,0),breaks=\(x)pretty(x,6))+
scale_color_manual(values=c("black","gray60",NA))+
scale_fill_manual(values=c(NA,NA,"gray60"))+
scale_linetype_manual(values=c("solid","solid",NA))+
scale_alpha_manual(values=c(1,0,0),guide="none")+
guides(color=guide_legend(order=1),linetype=guide_legend(order=1),fill=guide_legend(order=2,keyheight=unit(5,"pt")))+
coord_cartesian(clip="off")+
theme(axis.text=element_text(size=7,color="black"),
  axis.ticks=element_line(linewidth=.3,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(3,"pt"),
  axis.title=element_text(size=8),
  legend.background=element_blank(),
  legend.box="horizontal",
  legend.box.margin=margin(),
  legend.box.spacing=unit(0,"pt"),
  legend.direction="horizontal",
  legend.justification="right",
  legend.key=element_blank(),
  legend.key.height=unit(5,"pt"),
  legend.key.width=unit(19,"pt"),
  legend.margin=margin(,,3),
  legend.position="top",
  legend.spacing.x=unit(2,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=7,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.3),
  panel.spacing.y=unit(2,"pt"),
  plot.margin=margin(4,4,4,4),
  plot.subtitle=element_text(size=7,margin=margin(1,,3)),
  plot.title=element_text(size=7.4,face=2,margin=margin(1,,3)),
  strip.background=element_blank(),
  strip.text=element_text(size=7,face=2,margin=margin(,,2)))
ggsave("1.png",width=4.5,height=3.9,dpi=400*4)
system("magick 1.png -resize 25% -colorscale gray -colors 16 1.png")

library(data.table);library(ggplot2)

v=do.call(rbind,lapply(2011:2019,\(i)fread(paste0(i,".csv.gz"))[restatus!=4,.(.I,year,month=monthdth,age,ucd=ucod,mcd=unlist(.SD,,F)),.SDcols=patterns("enicon_")][mcd%like%"C34"]))
v=v[,age:=pmin(100,ifelse(age%/%1000%in%c(1,9),age%%1000,0))]

a=v[ucd%like%"C34"&rowid(I,year)==1,.(type=3,y=.N),.(age,year,month)]
a=rbind(a,v[mcd%like%"C34"][rowid(I,year)==1,.(type=1,y=.N),.(age,year,month)])
dim=lapply(a[,-5],unique);dim$year=1999:2025
a=merge(do.call(CJ,dim),a,all.x=T)[is.na(y),y:=0]
a=merge(fread("http://sars2.net/f/uspopdeadmonthly.csv"),a)
a[,x:=year*12+month]
a=merge(a[year%in%2011:2019,.(x=unique(a$x),base=predict(lm(y/persondays~x),.(x=unique(a$x)))),.(age,type)],a)

t=fread("http://sars2.net/f/wonderlungcancerallage.csv")[!(year==2024&month>=11)]
p=dcast(t,year+month~type,value.var="dead")
p=p[,.(x=year*12+month,year,month,y=c(mcd,mcd-nafill(and,,0),ucd),type=rep(1:3,each=.N),group=1)]
p=p[!(x<t[type=="and",min(year*12+month)]&type==2)]
p[,y:=y/ifelse(month==2,28+(year%%4==0),(31:30)[rep(1:2,,13)[-8]][month])]

base=merge(p,a[year>=2011,.(base=sum(base*pop)),.(x,type)])
base2=p[type!=2&year>=2011]
base2$base=p[year%in%2018:2019&month%in%4:9,predict(lm(y~x),base2[type==3]),type]$V1
base2[,base:=base+.SD[year%in%2018:2019,mean(y)-mean(base)]]
p=rbind(p,rbind(base[,group:=2],base2[,group:=3])[,.(year,month,x,y=base,group,type)])

p[,type:=factor(type,,c("MCD C34","MCD C34 but not COVID","UCD C34"))]
p[,group:=factor(group,,c("Actual deaths","Baseline from 2011-2019 trend in CMR by age","Baseline from 2018-2019 trend in raw deaths"))]

xstart=min(p$x);xend=2025*12;p=p[x%in%xstart:xend]
xbreak=seq(xstart+6,xend,12)
ybreak=pretty(p$y,7);ystart=ybreak[1];yend=max(ybreak)

ggplot(p,aes(x,y))+
geom_vline(xintercept=seq(xstart-.5,xend+.5,12),color="gray90",linewidth=.4)+
geom_hline(yintercept=c(ystart,yend),linewidth=.4,lineend="square")+
geom_vline(xintercept=c(xstart-.5,xend+.5),linewidth=.4,lineend="square")+
geom_line(aes(color=type,linetype=group),linewidth=.5)+
geom_point(aes(color=type,alpha=group),stroke=0,size=.8)+
labs(title="NVSS: Monthly deaths divided by number of days in month for cause C34
(malignant neoplasm of bronchus and lung)",x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart-.5,xend+.5),breaks=xbreak,labels=substr(xbreak%/%12,3,4))+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak)+
scale_color_manual(values=c("#44f","#bbf","black"))+
scale_linetype_manual(values=c("solid","42","11"))+
scale_alpha_manual(values=c(1,0,0))+
coord_cartesian(clip="off",expand=F)+
guides(color=guide_legend(order=1),linetype=guide_legend(order=2),alpha=guide_legend(order=2))+
theme(axis.text=element_text(color="black",size=11),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(4,"pt"),
  legend.background=element_rect(color="black",linewidth=.4),
  legend.justification=c(0,0),
  legend.key=element_blank(),
  legend.key.height=unit(10,"pt"),
  legend.key.width=unit(28,"pt"),
  legend.position=c(0,0),
  legend.spacing.x=unit(2,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.margin=margin(4,5,4,4),
  legend.text=element_text(size=11,vjust=.7),
  legend.title=element_blank(),
  panel.background=element_blank(),
  plot.margin=margin(5,5,4,5),
  plot.title=element_text(size=11.6,face=2,margin=margin(4,,4),lineheight=.9))
ggsave("1.png",width=6.4,height=3.5,dpi=300*4)
sub="\u00a0     The dashed baseline was calculated by multiplying the linear trend in CMR for each age in 2011-2019 with the population size of the age. The slope of the dotted baseline was determined based on April to September only but the intercept was determined based on all months.
      Monthly resident population estimates by single year of age are from www2.census.gov/
programs-surveys/popest/datasets/2020-2023/national/asrh/nc-est2023-alldata-r-file0{1..8}.csv, and from the corresponding files in the 2010-2020 directory. In order to get rid of a sudden jump in population size after the switch from the 2010-2020 to 2020-2023 estimates, the 2010-based estimates were multiplied by a slope so that they matched the 2020-based estimates in April 2020 when the estimates were merged. Vintage 2023 population projections were used for 2024."
system(paste0("f=1.png;mar=70;w=`identify -format %w $f`;magick \\( $f -chop 0x60 \\) \\( -size $[w-mar*2]x -font Arial -interline-spacing -3 -pointsize $[43*4] caption:'",gsub("'","'\\\\'",sub),"' -gravity southwest -splice $[mar]x80 \\) -append -resize 25% -colors 256 1.png"))

library(data.table);library(ggplot2)

v=do.call(rbind,lapply(2011:2019,\(i)fread(paste0(i,".csv.gz"))[restatus!=4,.(.I,year,month=monthdth,age,ucd=ucod,mcd=unlist(.SD,,F)),.SDcols=patterns("enicon_")][mcd%like%"C"]))
v=v[,age:=pmin(100,ifelse(age%/%1000%in%c(1,9),age%%1000,0))]

a=v[ucd%like%"C"&rowid(I,year)==1,.(type=3,y=.N),.(age,year,month)]
a=rbind(a,v[mcd%like%"C"][rowid(I,year)==1,.(type=1,y=.N),.(age,year,month)])
dim=lapply(a[,-5],unique);dim$year=1999:2025
a=merge(do.call(CJ,dim),a,all.x=T)[is.na(y),y:=0]
a=merge(fread("http://sars2.net/f/uspopdeadmonthly.csv"),a)
a[,x:=year*12+month]
a=merge(a[year%in%2011:2019,.(x=unique(a$x),base=predict(lm(y/persondays~x),.(x=unique(a$x)))),.(age,type)],a)

t=fread("http://sars2.net/f/wondermalignantallage.csv")[!(year==2024&month>=11)]
p=dcast(t,year+month~type,value.var="dead")
p=p[,.(x=year*12+month,year,month,y=c(mcd,mcd-nafill(and,,0),ucd),type=rep(1:3,each=.N),group=1)]
p=p[!(x<t[type=="and",min(year*12+month)]&type==2)]
p[,y:=y/ifelse(month==2,28+(year%%4==0),(31:30)[rep(1:2,,13)[-8]][month])]

base=merge(p,a[year>=2011,.(base=sum(base*pop)),.(x,type)])
base2=p[type!=2&year>=2011]
base2$base=p[year%in%2018:2019&month%in%4:9,predict(lm(y~x),base2[type==3]),type]$V1
base2[,base:=base+.SD[year%in%2018:2019,mean(y)-mean(base)]]
p=rbind(p,rbind(base[,group:=2],base2[,group:=3])[,.(year,month,x,y=base,group,type)])

p[,type:=factor(type,,c("MCD C00-C97","MCD C00-C97 but not COVID","UCD C00-C97"))]
p[,group:=factor(group,,c("Actual deaths","Baseline from 2011-2019 trend in CMR by age","Baseline from 2018-2019 trend in raw deaths"))]

xstart=min(p$x);xend=2024*12;p=p[x%in%xstart:xend]
xbreak=seq(xstart+6,xend,12)
ybreak=pretty(p$y,7);ystart=min(p$y)-10;yend=max(p$y)+10

note="Ethical Skeptic calculates excess MCD cancer deaths using the equivalent of the light blue line but Akston uses the equivalent of the dark blue line"|>stringr::str_wrap(40)

note2="Ethical Skeptic's baseline without PFE adjustment is similar to my dotted baseline. Compared to the PFE adjustment used by ES, my dashed baseline employs a more accurate method to account for the reduction in population size due to excess deaths since 2020."|>stringr::str_wrap(30)

ggplot(p,aes(x,y))+
geom_vline(xintercept=seq(xstart-.5,xend+.5,12),color="gray90",linewidth=.4)+
geom_hline(yintercept=c(ystart,yend),linewidth=.4,lineend="square")+
geom_vline(xintercept=c(xstart-.5,xend+.5),linewidth=.4,lineend="square")+
geom_line(aes(color=type,linetype=group),linewidth=.5)+
geom_point(aes(color=type,alpha=group),stroke=0,size=.8)+
annotate(geom="label",x=2009*12+6,y=yend-10,vjust=1,label=note,label.r=unit(2,"pt"),label.padding=unit(3,"pt"),label.size=.2,size=3.6,lineheight=.82,hjust=0)+
annotate(geom="label",x=1999*12+9,y=yend-10,vjust=1,label=note2,label.r=unit(2,"pt"),label.padding=unit(3,"pt"),label.size=.2,size=3.6,lineheight=.82,hjust=0)+
labs(title="NVSS, malignant neoplasms (C00-C97): Monthly deaths divided by days in month",x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart-.5,xend+.5),breaks=xbreak,labels=substr(xbreak%/%12,3,4))+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak)+
scale_color_manual(values=c("#44f","#bbf","black"))+
scale_linetype_manual(values=c("solid","42","11"))+
scale_alpha_manual(values=c(1,0,0))+
coord_cartesian(clip="off",expand=F)+
guides(color=guide_legend(order=1),linetype=guide_legend(order=2),alpha=guide_legend(order=2))+
theme(axis.text=element_text(color="black",size=11),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(4,"pt"),
  legend.background=element_rect(color="black",linewidth=.4),
  legend.box="horizontal",
  legend.box.just="bottom",
  legend.box.spacing=unit(0,"pt"),
  legend.direction="vertical",
  legend.justification=c(1,0),
  legend.key=element_blank(),
  legend.key.height=unit(12,"pt"),
  legend.key.width=unit(28,"pt"),
  legend.margin=margin(4,5,4,4),
  legend.position="top",
  legend.spacing.x=unit(0,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=11,vjust=.7),
  legend.title=element_blank(),
  panel.background=element_blank(),
  plot.margin=margin(6,6,5,6),
  plot.title=element_text(size=11.3,face=2,margin=margin(4,,4),lineheight=.9))
ggsave("1.png",width=6.9,height=4.1,dpi=300*4)
sub="\u00a0     The dashed baseline was calculated by multiplying the linear trend in CMR for each age in 2011-2019 with the population size of the age. The slope of the dotted baseline was determined based on April to September only but the intercept was determined based on all months.
      Monthly resident population estimates by single year of age are from www2.census.gov/
programs-surveys/popest/datasets/2020-2023/national/asrh/nc-est2023-alldata-r-file0{1..8}.csv, and from the corresponding files in the 2010-2020 directory. In order to get rid of a sudden jump in population size after the switch from the 2010-2020 to 2020-2023 estimates, the 2010-based estimates were multiplied by a slope so that they matched the 2020-based estimates in April 2020 when the estimates were merged."
system(paste0("f=1.png;mar=70;w=`identify -format %w $f`;magick \\( $f -chop 0x60 \\) \\( -size $[w-mar*2]x -font Arial -interline-spacing -3 -pointsize $[43*4] caption:'",gsub("'","'\\\\'",sub),"' -gravity southwest -splice $[mar]x80 \\) -append -resize 25% -colors 256 1.png"))

library(data.table);library(ggplot2)

t=fread("http://sars2.net/f/wondernatural.csv")[cause=="A-Q"&date<="2024-06"]
p=na.omit(t[,.(dead=sum(dead)),.(x=as.Date(paste0(date,"-1")))])

p[,dead:=dead/days_in_month(x)*30]

p$linear=p[year(x)<2020,predict(lm(dead~x),p)]
p[,excess:=cumsum(dead-linear)]

p2=copy(p)[,x:=as.integer(x)]
fit=nls(excess~A*sin(B*x+C)+D,p2[year(p$x)<2020],
  list(A=(max(p2$excess)-min(p2$excess))/2,B=2*pi/(max(p2$x)-min(p2$x)),C=0,D=mean(p2$excess)),
  nls.control(maxiter=100,tol=1e-6))
p$sine=predict(fit,p2)

p=melt(p,id=1,,"z","y")
p[,z:=factor(z,unique(z),c("Deaths","1999-2019 linear trend in deaths","Cumulative excess relative to linear trend","Sine wave fitted to pre-2020 cumulative excess"))]

xstart=as.Date("1999-1-1");xend=as.Date("2025-1-1")
p=p[x>=xstart&x<=xend]
xbreak=seq(xstart+182,xend,"year");xlab=year(xbreak)
ybreak=pretty(p$y,7);ystart=ybreak[1];yend=max(ybreak)

ggplot(p,aes(x=x+14,y))+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray83",linewidth=.4)+
geom_hline(yintercept=0,linewidth=.4)+
geom_line(aes(color=z),linewidth=.5)+
labs(x=NULL,y=NULL,title="CDC WONDER: Monthly deaths with underlying cause A-Q",subtitle="Divided by number of days in month and multiplied by 30. Excludes UCD external causes, COVID, and R and X codes."|>stringr::str_wrap(74))+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak,labels=\(x)ifelse(x==0,0,paste0(x/1e3,"k")))+
scale_color_manual(values=c("black","#6666ff","gray60","#ff66ff","#ff8888","#ff66ff"))+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=11,color="black"),
  axis.text.x=element_text(angle=90,vjust=.5,hjust=1),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(5,"pt"),
  legend.background=element_blank(),
  legend.box.background=element_rect(color="black",linewidth=.4),
  legend.direction="vertical",
  legend.justification=c(0,1),
  legend.key=element_blank(),
  legend.key.height=unit(12,"pt"),
  legend.key.width=unit(26,"pt"),
  legend.margin=margin(4,4,4,4),
  legend.position=c(0,1),
  legend.spacing.x=unit(3,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.3),
  plot.margin=margin(6,6,5,5),
  plot.subtitle=element_text(size=11,margin=margin(,,5)),
  plot.title=element_text(size=11.5,face=2,margin=margin(2,,4)))
ggsave("1.png",width=6.3,height=4,dpi=300*4)
system("magick 1.png -resize 25% PNG8:1.png")

library(data.table);library(ggplot2)

xstart=1999;xend=2080

t=fread("http://sars2.net/f/uspopdead.csv")
p=t[,.(y=sum(dead),z="Actual deaths"),.(x=year)]
cmr=t[year==2010,dead/pop]
pop=t[year==2010,pop]
births=pop[1]
sim=data.table(x=NA,y=NA)
for(year in 2010:xend){
  died=pop*cmr;pop=pop-died
  pop=c(births,pop[1:99],sum(pop[100:101]))
  sim=rbind(sim,list(year,sum(died)))
}

p=rbind(p,p[,.(x=xstart:xend,y=p[x<2020,predict(lm(y~x),.(x=xstart:xend))],z="1999-2019 linear trend")])
p=rbind(p,sim[,z:="Modeled deaths"])
p[,z:=factor(z,unique(z))]

excess=merge(p[z==levels(z)[1]],p[z==levels(z)[2],.(x,y2=y)])[x<2020,.(x,y=cumsum(y-y2))]
fit=nls(y~A*sin(B*x+C)+D,excess,
  list(A=excess[,max(y)-min(y)]/2,B=2*pi/excess[,max(x)-min(x)],C=0,D=mean(excess$y)),
  nls.control(maxiter=100,tol=1e-6))

pred=diff(predict(fit,.(x=1998:2080)))[-1]
p=rbind(p,p[z==levels(z)[2]][,y:=y+pred][,z:="1999-2019 plus sine"])

xbreak=seq(2000,xend,5)
ybreak=pretty(p$y,7);ystart=ybreak[1];yend=max(ybreak)

note=stringr::str_wrap("The model started in 2010 so that the population size of each age was set to the vintage 2020 mid-2010 resident population estimates. Each year the same fraction of people for each age killed as the 2010 CMR value for the age, the age of the remaining people was incremented by one, and the same number of people aged zero were added as in the mid-2010 resident population estimates.",40)

ggplot(p,aes(x=x,y))+
geom_hline(yintercept=0,linewidth=.4)+
geom_line(aes(color=z),linewidth=.5)+
labs(x=NULL,y=NULL,title="Model for yearly deaths in United States")+
scale_x_continuous(limits=c(xstart-.5,xend+.5),breaks=xbreak)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak,labels=\(x)paste0(x/1e6,"M"))+
scale_color_manual(values=c("black","#6666ff","#ff3333","#ff55ff"))+
annotate(geom="label",x=2039,y=32e5,vjust=1,label=note,label.r=unit(2,"pt"),label.padding=unit(3,"pt"),label.size=.2,size=3.6,lineheight=.82,hjust=0)+
annotate(geom="text",x=2017,y=2.4e6,label="The magenta line is\nsimilar to Akston's\nbaseline with\nlong-term correction",hjust=0,size=3.6,lineheight=.93,color="#ff55ff")+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=11,color="black"),
  axis.text.x=element_text(angle=90,vjust=.5,hjust=1),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length=unit(5,"pt"),
  legend.background=element_blank(),
  legend.box.background=element_rect(color="black",linewidth=.4),
  legend.box.spacing=unit(0,"pt"),
  legend.direction="vertical",
  legend.justification=c(0,1),
  legend.key=element_blank(),
  legend.key.height=unit(14,"pt"),
  legend.key.width=unit(26,"pt"),
  legend.position=c(0,1),
  legend.spacing.x=unit(3,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.margin=margin(5,5,5,5),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.3),
  plot.margin=margin(6,8,5,5),
  plot.title=element_text(size=11.5,face=2,margin=margin(1,,5)))
ggsave("1.png",width=6.3,height=4,dpi=300*4)
system("magick 1.png -resize 25% PNG8:1.png")

library(data.table);library(ggplot2)

cul=\(x,y)y[cut(x,c(y,Inf),,T,F)]
agecut=\(x,y)cut(x,c(y,Inf),paste0(y,c(paste0("-",y[-1]-1),"+")),T,F)

t=fread("http://sars2.net/f/cancerincidence.csv")

pop=fread("http://sars2.net/f/uspopdead.csv")[year>=2010,.(pop=sum(pop)),.(year,age=cul(age,unique(t$age)))]
pop2=fread("https://www2.census.gov/programs-surveys/popest/datasets/2000-2010/intercensal/national/us-est00int-alldata.csv")
pop=rbind(pop,pop2[MONTH==7&YEAR<2010&AGE!=999,.(pop=sum(TOT_POP)),.(year=YEAR,age=cul(AGE,unique(t$age)))])

t=merge(pop[year==2020,.(age,std=pop/sum(pop))],merge(t,pop))

ages=c(0,30,50,70,85)

p=t[,.(y=sum(cases/pop*std*1e5)),.(x=year,group=agecut(age,ages))]
p=rbind(p,t[,.(y=sum(cases/pop*std*1e5),group="Total"),.(x=year)])

xstart=2000;xend=2021;xbreak=xstart:xend
base=p[x%in%2010:2019,.(x=xbreak,y=predict(lm(y~x),.(x=xbreak))),group]
p$z="Age-standardized rate"
p=rbind(p,base[x%in%2010:2019][,z:="2010-2019 linear trend"],base[,z:="Trend projection"])
p[,z:=factor(z,unique(z))]

expand=p[,.(min=min(y),max=max(y)),group]
rat=expand[,max(max/min)]*1.02
expand=expand[,{x=(rat*min-max)/(1+rat);.(y=c(min-x,max+x),group)}]

ggplot(p,aes(x,y))+
facet_wrap(~group,ncol=2,dir="v",scales="free_y")+
geom_line(aes(linetype=z,color=z),linewidth=.6)+
geom_point(aes(alpha=z),stroke=0,size=1.2)+
geom_line(data=p[z==z[1]],aes(linetype=z,color=z),linewidth=.6)+
geom_text(data=rbind(expand,p[,.(y,group)])[,max(y),group],aes(y=V1,label=group,x=pmean(xstart,xend)),fontface=2,size=3.5,vjust=1.7)+
geom_point(data=expand,x=xstart,alpha=0)+
labs(title="CDC WONDER: Age-standardized yearly new cancer cases per 100,000 people",subtitle=paste0("Source: wonder.cdc.gov/cancer.html. The age-standardized rate was calculated by 5-year age groups up to ages 85+ so that the 2020 resident population estimates were used as the standard population. The population estimates are from census.gov/programs-surveys/popest/datasets. In order to get rid of a sudden jump at the point when the 2010-based and 2020-based population estimates were merged, the 2010-based estimates were multiplied by a linear slope so that they were equal to the 2020-based estimates in 2020. Intercensal population estimates are used for 2000-2009. All panels have the same ratio between the maximum and minimum point of the y-axis.")|>stringr::str_wrap(88),x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart-.5,xend+.5),breaks=xbreak)+
scale_y_continuous(breaks=\(x)pretty(x,4))+
coord_cartesian(clip="off",expand=F)+
scale_color_manual(values=c("black","#5555ff","#5555ff"))+
scale_alpha_manual(values=c(1,0,0))+
scale_linetype_manual(values=c("solid","solid","11"))+
theme(axis.text=element_text(size=11,color="black"),
  axis.text.x=element_text(angle=90,vjust=.5,hjust=1),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length=unit(4,"pt"),
  axis.title=element_text(size=8),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.key=element_blank(),
  legend.key.height=unit(13,"pt"),
  legend.key.width=unit(26,"pt"),
  legend.position="top",
  legend.justification="left",
  legend.spacing.x=unit(2,"pt"),
  legend.margin=margin(,,6),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.4),
  panel.spacing.x=unit(4,"pt"),
  panel.spacing.y=unit(3,"pt"),
  plot.margin=margin(7,7,7,7),
  plot.subtitle=element_text(size=11,margin=margin(,,4)),
  plot.title=element_text(size=11.5,face=2,margin=margin(4,,6)),
  strip.background=element_blank(),
  strip.text=element_blank())
ggsave("1.png",width=7,height=6,dpi=300*4)
system("magick 1.png -resize 25% -colors 256 1.png")

library(data.table);library(ggplot2)

download.file("https://data.cdc.gov/api/views/3yf8-kanr/rows.csv?accessType=DOWNLOAD","Weekly_Counts_of_Deaths_by_State_and_Select_Causes__2014-2019.csv")

t=fread("Weekly_Counts_of_Deaths_by_State_and_Select_Causes__2014-2019.csv")

t=t[t[[1]]=="United States",.(year=`MMWR Year`,week=`MMWR Week`,dead=`All  Cause`,date=as.IDate(`Week Ending Date`,"%m/%d/%Y"))]

t=rbind(t,fread("http://sars2.net/f/wonderdeadweekly.csv")[year>2019&!(year==2024&week>=25)&week!=99])

t$trough=t[year<2020&week%in%21:31,predict(lm(dead~date),t)]
base=merge(t,t[year<2020&week%in%21:31,.(week=26,base=mean(dead)),year])
t$trough2=predict(lm(base~date,base),t)

t$trough3=t[year<2020&week%in%16:36,predict(lm(dead~date),t)]
base2=merge(t,t[year<2020&week%in%16:36,.(week=26,base=mean(dead)),year])
t$trough4=predict(lm(base~date,base2),t)

lab=c("Actual deaths","Linear regression of weeks 21-31 of 2014-2019 (each week included in regression)","Linear regression of weeks 21-31 of 2014-2019 (regression of yearly averages plotted on week 26)","Linear regression of weeks 16-36 of 2014-2019 (each week included in regression)","Linear regression of weeks 16-36 of 2014-2019 (regression of yearly averages plotted on week 26)")
p=t[,.(x=date,y=c(dead,trough,trough2,trough3,trough4),z=factor(rep(lab,each=.N),lab))]

xstart=as.Date("2014-1-1");xend=as.Date("2025-1-1")
xbreak=seq(xstart,xend,"6 month");xlab=ifelse(month(xbreak)==7,year(xbreak),"")
ybreak=pretty(p$y,7);ystart=ybreak[1];yend=max(ybreak)

ggplot(p,aes(x,y))+
geom_vline(xintercept=c(xstart),color="magenta",linewidth=.2,lineend="square")+
geom_hline(yintercept=c(ystart),color="magenta",linewidth=.2,lineend="square")+
geom_line(aes(color=z,linetype=z),linewidth=.3)+
labs(x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart-.2,xend+.5),breaks=xbreak,labels=xlab)+
scale_y_continuous(breaks=ybreak,labels=\(x)paste0(x/1e3,"k"))+
scale_color_manual(values=hsv(c(30,12,12,22,22)/36,.6,1))+
scale_linetype_manual(values=c("solid","42","11","42","11"))+
coord_cartesian(ylim=c(ystart,yend),clip="off",expand=F)+
theme(axis.text=element_text(size=6,color="magenta"),
  axis.ticks=element_line(linewidth=.2,color="magenta"),
  axis.ticks.x=element_line(color=alpha("magenta",c(1,0))),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.key=element_blank(),
  legend.key.height=unit(7,"pt"),
  legend.key.width=unit(18,"pt"),
  legend.position=c(.01,.46),
  legend.justification=c(0,0),
  legend.spacing.x=unit(2,"pt"),
  legend.margin=margin(,,6),
  legend.text=element_text(color="white",size=5.4,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.grid=element_blank(),
  plot.background=element_rect(fill="transparent",color=NA))
ggsave("1.png",width=7,height=2,dpi=400)

library(data.table);library(ggplot2)

t=fread("http://sars2.net/f/uspopdeadmonthly.csv")
mult=7/(31+30+31+31)
t=t[month%in%5:8,.(dead=sum(dead)*mult,pop=sum(persondays)*mult),.(year,age)]
t=merge(t[year%in%2014:2019,.(year=1999:2024,base=predict(lm(dead/pop~year),.(year=1999:2024))),age],t)
lab=c("Deaths in May to August divided by 31+30+31+31 and multiplied by 7","2014-2019 trend in CMR by age multiplied by population size (May to August only)")
p=t[,.(y=c(sum(dead),sum(base*pop)),z=lab),.(x=as.Date(paste0(year,"-7-1")))]
p[,z:=factor(z,lab)]

xstart=as.Date("2014-1-1");xend=as.Date("2024-1-1")
p=p[x%in%xstart:xend]
xbreak=seq(xstart,xend,"6 month");xlab=ifelse(month(xbreak)==7,year(xbreak),"")
ybreak=4:8*1e4;ystart=4e4;yend=8e4

ggplot(p,aes(x,y))+
geom_line(aes(linewidth=z,color=z,linetype=z))+
geom_point(aes(color=z,alpha=z),size=1.4,stroke=.3,shape=1)+
labs(x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak,labels=\(x)paste0(x/1e3,"k"))+
scale_color_manual(values=c("#ff6666",hsv(12/36,.6,1)))+
scale_linetype_manual(values=c("solid","solid"))+
scale_linewidth_manual(values=c(0,.3))+
scale_alpha_manual(values=c(1,0))+
coord_cartesian(clip="off",expand=F)+
theme(axis.line=element_line(linewidth=.3,color="white",lineend="square"),
  axis.text=element_text(size=6,color="white"),
  axis.ticks=element_line(linewidth=.2,color="white"),
  axis.ticks.x=element_line(color=alpha("white",c(1,0))),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.key=element_blank(),
  legend.key.height=unit(7,"pt"),
  legend.key.width=unit(18,"pt"),
  legend.position=c(.01,.6),
  legend.justification=c(0,0),
  legend.spacing.x=unit(2,"pt"),
  legend.margin=margin(,,6),
  legend.text=element_text(color="white",size=5.4,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.grid=element_blank(),
  plot.background=element_rect(fill="transparent",color=NA))
ggsave("1.png",width=7,height=2.05,dpi=400)

library(data.table);library(ggplot2)

t=fread("https://raw.githubusercontent.com/owid/covid-19-data/refs/heads/master/public/data/excess_mortality/excess_mortality.csv")
o=fread("owid-covid-data.csv")

t2=t[,.(location,date,cum=cum_excess_per_million_proj_all_ages)]
t2=na.omit(t2[location%in%o[continent=="Europe",unique(location)]])
cum=t2[date%like%"2023-12"][rev(order(date))][rowid(location)==1,.(location,cum2=cum)]
cum=merge(cum,t2[date%like%"2022-12"][rev(order(date))][rowid(location)==1,.(location,cum)])

vax=na.omit(o[,.(date,location,vax=people_vaccinated_per_hundred)])[order(-date)][rowid(location)==1]
p=merge(cum,vax)[,.(x=cum2-cum,y=vax,location)]

east=strsplit("Bulgaria,Croatia,Czechia,Estonia,Hungary,Latvia,Poland,Romania,Serbia,Slovenia,Slovakia,Russia,Kosovo,Montenegro,Albania,North Macedonia,Bosnia and Herzegovina",",")[[1]]
lab=c("Not eastern bloc","Eastern bloc")
p[,z:=factor(ifelse(location%in%east,lab[2],lab[1]),lab)]

xbreak=p[,pretty(x,7)];ybreak=p[,pretty(y,7)]

ggplot(p,aes(x,y))+
coord_cartesian(clip="off",expand=F)+
geom_smooth(method="lm",formula=y~x,linewidth=.5,se=F,color="black",linetype="42")+
geom_vline(xintercept=0,linewidth=.4,color="gray60")+
geom_hline(yintercept=0,linewidth=.4,color="gray60")+
geom_point(aes(color=z),size=.6)+
ggrepel::geom_text_repel(aes(label=location,color=z),size=3,max.overlaps=Inf,segment.size=.3,min.segment.length=.2,box.padding=.07,show.legend=F)+
labs(title=sprintf("OWID: Excess deaths in December 2022 to December 2023 vs\npercentage of vaccinated people (r ≈ %.2f)",p[,cor(x,y)]
),x="Excess deaths per million in December 2022 to December 2023",y="Percentage of vaccinated people at OWID")+
scale_x_continuous(limits=range(xbreak),breaks=xbreak)+
scale_y_continuous(limits=range(ybreak),breaks=ybreak,labels=\(x)paste0(x,"%"))+
scale_color_manual(values=c("#5555ff","#ff4444"))+
theme(axis.text=element_text(size=11,color="black"),
  axis.text.x=element_text(margin=margin(3)),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length=unit(4,"pt"),
  axis.title=element_text(size=11),
  axis.title.x=element_text(margin=margin(3)),
  axis.title.y=element_text(margin=margin(,2)),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.justification="right",
  legend.key=element_blank(),
  legend.key.height=unit(13,"pt"),
  legend.key.width=unit(26,"pt"),
  legend.position="top",
  legend.spacing.x=unit(2,"pt"),
  legend.margin=margin(,,6),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.4),
  panel.grid.major=element_blank(),
  plot.margin=margin(4,18,4,4),
  plot.subtitle=element_text(size=11),
  plot.title=element_text(size=11,face=2,margin=margin(1,,4),hjust=0))
ggsave("1.png",width=6.2,height=4,dpi=300*4)
system("mogrify 1.png -trim -resize 25% -bordercolor white -border 30 -colors 256 1.png")

library(data.table);library(ggplot2)

o=fread("owid-covid-data.csv")

isoweek=\(year,week,weekday=7){d=as.Date(paste0(year,"-1-4"));d-(as.integer(format(d,"%w"))+6)%%7-1+7*(week-1)+weekday}
q=\(...)as.character(substitute(...()))
ma=\(x,b=1,f=b){x[]=rowMeans(embed(c(rep(NA,b),x,rep(NA,f)),f+b+1),na.rm=T);x}

euro=fread("demo_r_mwk_05.tsv")
meta=fread(text=euro[[1]],header=F)
pick=meta[,V2%like%"^Y"&V3=="T"];meta=meta[pick];euro=euro[pick]
meta[,age:=as.integer(sub("\\D+(\\d+).*","\\1",sub("Y_LT5",0,V2)))]
eu=meta[,.(loc=V5,age,date=rep(names(euro)[-1],each=.N),dead=as.integer(sub("\\D+","",unlist(euro[,-1],,F))))]
eu[,year:=as.integer(substr(date,1,4))][,week:=as.integer(substr(date,7,8))]
eu[,date:=isoweek(year,week,4)]
eu=eu[loc%in%eu[year%in%2013:2023,.N,loc][N==max(N),loc]]
eu=eu[year>=2013]

pop=fread("http://sars2.net/f/eurostatpopdead.csv.gz")[year%in%2013:2023]
a=pop[,.(pop=sum(pop)),.(location,name,date=as.Date(paste0(year,"-1-1")),age=pmin(age,90)%/%5*5)]
a=a[location%in%na.omit(a)[,.N,location][N==max(N),location]]
a=a[,spline(date,pop,xout=unique(eu$date),method="natural"),.(location,age,name)]

a=merge(eu,a[,.(date=`class<-`(x,"Date"),pop=y,loc=location,name,age)])
a=merge(pop[location=="EU27_2020"&year==2020,.(pop=sum(pop)),.(age=pmin(age,90)%/%5*5)][,.(age,std=pop/sum(pop))],a)
a=a[loc%in%na.omit(a)[year%in%2013:2023,.N,loc][N==max(N),loc]]

a=a[,.(dead=sum(dead),asmr=sum(dead/pop*std/7*365e5)),,.(name,date,year,week)]
slope=a[year%in%2013:2019&week%in%15:35,mean(asmr),.(year,name)]
slope=slope[,{x=.(year=2013:2024,slope=predict(lm(V1~year),.(year=2013:2024)));x$slope=x$slope/x$slope[x$year==2016];x},name]
weekly=a[year%in%2013:2019,.(weekly=mean(asmr)),.(week,name)]
a=merge(slope,merge(weekly,a))
a=merge(a[year%in%2016:2019,.(deadbase=mean(dead)),.(name,week)],a)
p=a[year==2023,.(pct2=(sum(asmr)/sum(slope*weekly)-1)*100),name]

vax=o[rev(order(date))][,vax:=people_vaccinated_per_hundred][!is.na(vax)][rowid(location)==1,.(vax,name=location)]
p=na.omit(merge(vax,p))[,.(x=pct2,y=vax,name)]

east=q(Estonia,Latvia,Poland,Czechia,Slovakia,Hungary,Montenegro,Slovenia,Croatia,Serbia,Romania,Bulgaria)
p[,z:=ifelse(name%in%east,"Eastern bloc","Not eastern bloc")]
p[,z:=factor(z,unique(z))]

xbreak=pretty(p$x,8);ybreak=pretty(p$y,8)

ggplot(p,aes(x,y))+
geom_smooth(method="lm",formula=y~x,linewidth=.5,se=F,color="black",linetype="42")+
coord_cartesian(clip="off",expand=F)+
geom_vline(xintercept=0,linewidth=.4,color="gray60")+
geom_hline(yintercept=0,linewidth=.4,color="gray60")+
geom_point(aes(color=z),size=.6)+
ggrepel::geom_text_repel(aes(label=name,color=z),size=3,max.overlaps=Inf,segment.size=.3,min.segment.length=.2,box.padding=.07,show.legend=F)+
labs(title=sprintf("Excess ASMR in 2023 vs percentage of vaccinated people (r ≈ %.2f)",p[,cor(x,y)]
),x="Excess ASMR in 2023 at Eurostat",y="Percentage of vaccinated people at OWID")+
scale_x_continuous(limits=range(xbreak),breaks=xbreak,labels=\(x)paste0(x,"%"))+
scale_y_continuous(limits=range(ybreak),breaks=ybreak,labels=\(x)paste0(x,"%"))+
scale_color_manual(values=c("#5555ff","#ff4444"))+
theme(axis.text=element_text(size=11,color="black"),
  axis.text.x=element_text(margin=margin(3)),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length=unit(4,"pt"),
  axis.title=element_text(size=11),
  axis.title.x=element_text(margin=margin(3)),
  axis.title.y=element_text(margin=margin(,2)),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.key=element_blank(),
  legend.key.height=unit(13,"pt"),
  legend.key.width=unit(26,"pt"),
  legend.position="top",
  legend.justification="right",
  legend.spacing.x=unit(2,"pt"),
  legend.margin=margin(,,6),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.4),
  panel.grid.major=element_blank(),
  plot.margin=margin(4,18,4,4),
  plot.title=element_text(size=11,face=2,margin=margin(1,,4),hjust=1))
ggsave("1.png",width=6.2,height=4,dpi=300*4)
system("mogrify 1.png -trim -resize 25% -bordercolor white -border 30 -colors 256 1.png")

download.file("https://covid.ourworldindata.org/data/owid-covid-data.csv","owid-covid-data.csv")

library(data.table);library(ggplot2)

o=fread("owid-covid-data.csv")

t=fread("http://sars2.net/f/eurostatpopdead.csv.gz")[year%in%2013:2022]
t=merge(t,t[location=="EU27_2020"&year==2020,.(age,std=pop/sum(pop))])
t=t[location%in%na.omit(t)[,.N,location][N==max(N),location]]
t=t[!location%in%c("DE_TOT")]

a=t[,.(dead=sum(dead),asmr=sum(dead/pop*std)*1e5),.(location,name,year)]
a$base=a[year<2020,predict(lm(asmr~year),.(year=2013:2022)),location]$V1
a$base2=a[year%in%2015:2019,predict(lm(dead~year),.(year=2013:2022)),location]$V1
a=merge(a,a[year%in%2015:2019,.(base3=mean(dead)),location])
p=a[year==2022,.(y1=(asmr/base-1)*100,y2=(dead/base2-1)*100,y3=(dead/base3-1)*100),name]

p=merge(p,o[year(date)==2022,.(vaxpct=mean(people_vaccinated_per_hundred,na.rm=T)),.(name=location)])
p=na.omit(p)

lab=c("2013-2019 linear regression for ASMR","2015-2019 linear regression for raw deaths","2015-2019 average for raw deaths")
p=p[,.(name,x=vaxpct,y=c(y1,y2,y3),facet=factor(rep(lab,each=.N),lab))]

q=\(...)as.character(substitute(...()))
east=q(Bulgaria,Croatia,Czechia,Estonia,Hungary,Latvia,Poland,Romania,Serbia,Slovenia,Slovakia)
lab=c("Eastern bloc","Not eastern bloc")
p$z=factor(ifelse(p$name%in%east,lab[1],lab[2]),lab)

levels(p$facet)=p[,cor(x,y),facet][,sprintf("Baseline is %s (r ≈ %.2f)",facet,V1)]

xbreak=pretty(p$x,7);ybreak=pretty(p$y,7);ylim=extendrange(p$y,,.02)

ggplot(p,aes(x,y))+
facet_wrap(~facet,dir="v",scales="free_x")+
coord_cartesian(clip="off",expand=F)+
geom_vline(xintercept=0,linewidth=.4,color="gray60")+
geom_hline(yintercept=0,linewidth=.4,color="gray60")+
geom_smooth(method="lm",formula=y~x,linewidth=.5,se=F,color="black",linetype="42")+
geom_label(data=p[rowid(facet)==1],aes(label=stringr::str_wrap(facet,40)),x=xbreak[1],y=ylim[2],lineheight=.9,hjust=0,vjust=1,size=3.7,label.r=unit(0,"pt"),label.padding=unit(6,"pt"),label.size=.4)+
geom_point(aes(color=z),size=.7)+
ggrepel::geom_text_repel(aes(label=name,color=z),size=3,max.overlaps=Inf,segment.size=.3,min.segment.length=.2,box.padding=.07,show.legend=F)+
labs(title="Correlation between percentage of vaccinated people in 2022 and\nexcess mortality percent in 2022",x="Average percentage of vaccinated people throughout 2022 at OWID",y="Excess percentage of deaths in 2022 at Eurostat")+
scale_x_continuous(breaks=xbreak,limits=range(xbreak),labels=\(x)paste0(x,"%"))+
scale_y_continuous(breaks=ybreak,limits=ylim,labels=\(x)paste0(x,"%"))+
scale_color_manual(values=c("#ff4444","#5555ff"))+
theme(axis.text=element_text(size=11,color="black"),
  axis.text.x=element_text(margin=margin(3)),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length=unit(4,"pt"),
  axis.title=element_text(size=11),
  axis.title.x=element_text(margin=margin(3)),
  axis.title.y=element_text(margin=margin(,2)),
  legend.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.key=element_blank(),
  legend.key.height=unit(13,"pt"),
  legend.key.width=unit(16,"pt"),
  legend.position="top",
  legend.justification="right",
  legend.spacing.x=unit(2,"pt"),
  legend.margin=margin(,,6),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.4),
  panel.grid.major=element_blank(),
  panel.spacing=unit(4,"pt"),
  plot.margin=margin(6,22,5,5),
  plot.subtitle=element_text(size=11),
  plot.title=element_text(size=11.5,face=2,margin=margin(1,,5)),
  strip.background=element_blank(),
  strip.text=element_blank())
ggsave("1.png",width=6,height=8,dpi=300*4)

sub="Source: covid.ourworldindata.org/data/owid-covid-data.csv and Eurostat datasets demo_magec and demo_pjan. ASMR was calculated by single year of age up to ages 100+ so that the 2020 EU population estimates were used as the standard population. Countries that had deaths or population estimates missing for any combination of age or year in 2013-2022 were omitted. Romania was omitted because it was missing vaccination data at OWID."
system(paste0("mogrify -trim 1.png;w=`identify -format %w 1.png`;magick 1.png \\( -size $[w]x -font Arial -interline-spacing -3 -pointsize $[44*4] caption:'",gsub("'","'\\\\'",sub),"' -splice x100 \\) -append -resize 25% -bordercolor white -border 30 -colors 256 1.png"))

download.file("https://ec.europa.eu/eurostat/api/dissemination/sdmx/2.1/data/demo_r_mwk_ts?format=TSV","demo_r_mwk_ts.tsv")

library(data.table);library(ggplot2)

isoweek=\(year,week,weekday=7){d=as.Date(paste0(year,"-1-4"));d-(as.integer(format(d,"%w"))+6)%%7-1+7*(week-1)+weekday}

euro=fread("demo_r_mwk_ts.tsv")
euro=euro[euro[[1]]%like%"^W,T"]
euro[[1]]=sub(".*,","",euro[[1]])
eu=melt(euro,id=1)
names(eu)=c("iso","date","dead")
eu[,dead:=as.integer(sub("\\D+","",dead))]
eu[,year:=as.integer(substr(date,1,4))]
eu[,week:=as.integer(substr(date,7,8))]

d=na.omit(eu[iso=="RO"][,x:=isoweek(year,week,4)])

yearly=d[,.(year=2015:2024)]
yearly$mult=d[year%in%2015:2019&week%in%10:40,sum(dead),year][,predict(lm(V1~year),yearly)]
yearly[,mult:=mult/mult[year==2017]]
d=merge(d,d[year%in%2015:2019,.(weekly=mean(dead)),week])
d=merge(yearly,d)

p=rbind(d[,.(x,y=dead,z="Actual deaths")],d[,.(x,y=weekly*mult,z="2015-2019 linear baseline")])
p[,z:=factor(z,unique(z))]

xstart=as.Date("2015-1-1");xend=as.Date("2025-1-1")
p=p[x%in%xstart:xend]
ylim=c(0,max(p$y)*1.02);ybreak=pretty(ylim,7)
xbreak=seq(xstart+182,xend,"year");xlab=year(xbreak)

pct=d[,.(pct=(sum(dead)/sum(weekly*mult)-1)*100),.(x=as.Date(paste0(year,"-7-1")))]

ggplot(p,aes(x=x,y))+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray84",linewidth=.4)+
geom_hline(yintercept=0,linewidth=.4)+
geom_line(aes(color=z),linewidth=.6)+
geom_text(data=pct,aes(x,y=0,label=paste0(ifelse(pct>0,"+",""),sprintf("%.1f",pct),"%")),vjust=-.7,size=3.5,color="gray40")+
labs(x=NULL,y=NULL,title="Weekly deaths in Romania")+
scale_x_continuous(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=ylim,breaks=ybreak)+
scale_color_manual(values=c("black","#6666ff"))+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=11,color="black"),
  axis.text.x=element_text(margin=margin(4)),
  axis.text.y=element_text(margin=margin(,3)),
  axis.ticks=element_line(linewidth=.4,color="black"),
  axis.ticks.length.x=unit(0,"pt"),
  axis.ticks.length.y=unit(4,"pt"),
  legend.background=element_rect(fill="white",color="black",linewidth=.4),
  legend.box.background=element_blank(),
  legend.box.spacing=unit(0,"pt"),
  legend.direction="vertical",
  legend.justification=c(0,1),
  legend.key=element_blank(),
  legend.key.height=unit(13,"pt"),
  legend.key.width=unit(26,"pt"),
  legend.margin=margin(5,6,5,5),
  legend.position=c(0,1),
  legend.spacing.x=unit(2,"pt"),
  legend.spacing.y=unit(0,"pt"),
  legend.text=element_text(size=11,vjust=.5),
  legend.title=element_blank(),
  panel.background=element_blank(),
  panel.border=element_rect(fill=NA,linewidth=.4),
  panel.spacing=unit(2,"pt"),
  plot.margin=margin(5,5,2,5),
  plot.subtitle=element_text(margin=margin(,,5)),
  plot.title=element_text(size=12,face=2,margin=margin(2,,5)))
ggsave("1.png",width=5.8,height=3.2,dpi=300*4)

sub="\u00a0    Source: Eurostat dataset demo_r_mwk_ts. The percentage above the year shows the yearly total excess deaths, where the year 2024 only consists of weeks 1-43.
     The baseline was each week was calculated by multiplying average deaths for the same week number in 2015-2019 with a yearly multiplier. The yearly multiplier was calculated by doing a linear regression for total deaths on weeks 10-40 of 2015-2019, and then the ratio between the projected regression was divided by the value of the regression in 2017."
system(paste0("f=1.png;mar=100;w=`identify -format %w $f`;magick \\( $f \\) \\( -size $[w-mar*2]x -font Arial -interline-spacing -3 -pointsize $[45*4] caption:'",gsub("'","'\\\\'",sub),"' -gravity southwest -splice $[mar]x80 \\) -append -resize 25% -colors 256 1.png"))