R/covid.deaths.R

Defines functions shade covid.deaths IN

Documented in covid.deaths

## https://data.cdc.gov/NCHS/Provisional-COVID-19-Death-Counts-by-Sex-Age-and-W/vsak-wrfu

IN<-function(x,y){
	x<-toupper(x)
	y<-toupper(y)
	x%in%y
}

covid.deaths<-function(
	age.group=c(
	"Under 1 year",
	"1-4 years",
	"5-14 years",
	"15-24 years",
	"25-34 years",
	"35-44 years",
	"45-54 years",
	"55-64 years",
	"65-74 years",
	"75-84 years",
	"85 years and over"),
	sex=c("Male","Female"),
	all.causes=TRUE,
	cumulative=FALSE,
	data=list(),
	xlim=c(60,366+365+365),
	bg="transparent",
	plot=c("standard","smooth","bar"),
	show=c("raw","per.capita","percent","percent.of.covid.deaths"),
	split.groups=TRUE,
	show.total.deaths=TRUE,
	palette=c("new","original"),
	...){
	if(length(age.group)>0){
		if(hasArg(ylim)){ 
			ylim<-list(...)$ylim
			y.lim<-ylim
		} else ylim<-NULL
		plot<-plot[1]
		show<-show[1]
		palette<-palette[1]
		if(!is.null(data$CovidDeaths)) CovidDeaths<-data$CovidDeaths
		else CovidDeaths<-read.csv("https://liamrevell.github.io/data/Provisional_COVID-19_Death_Counts_by_Sex__Age__and_Week.csv")
		if(!is.null(data$Age.Pop)) Age.Pop<-data$Age.Pop
		else Age.Pop<-read.csv("https://liamrevell.github.io/data/US_Population_by_Age.csv")
		if(length(age.group)>0 && length(sex)>0){
			DD<-CovidDeaths
			dates<-as.Date(DD$End.Week,format="%m/%d/%Y")
			ii<-intersect(which(dates>="2021-01-09"),which(dates<"2022-01-08")) ## &&dates<"2022-01-08"
			jj<-intersect(which(dates>="2022-01-08"),which(dates<"2023-01-07"))
			DD$MMWR.Week[ii]<-DD$MMWR.Week[ii]+53
			DD$MMWR.Week[jj]<-DD$MMWR.Week[jj]+53+52
			ii<-which(IN(DD$Age.Group,age.group))
			DD<-DD[ii,]
			ii<-which(DD$Sex%in%sex)
			DD<-DD[ii,]
			mmwr<-sort(unique(DD$MMWR.Week))
			mmwr<-mmwr[6:length(mmwr)]
			cd<-sapply(mmwr,function(week,DD) 
				sum(DD$COVID.19.Deaths[which(DD$MMWR.Week==week)]),
				DD=DD)
			td<-sapply(mmwr,function(week,DD) 
				sum(DD$Total.Deaths[which(DD$MMWR.Week==week)]),
				DD=DD)-cd
			CD<-setNames(vector(mode="list",length=length(age.group)),age.group)
			TD<-setNames(vector(mode="list",length=length(age.group)),age.group)
			for(i in 1:length(age.group)){
				CD[[i]]<-matrix(NA,length(cd),length(sex),dimnames=list(mmwr,sex))
				TD[[i]]<-matrix(NA,length(td),length(sex),dimnames=list(mmwr,sex))
				for(j in 1:length(sex)){
					ii<-which(IN(DD$Age.Group,age.group[i]))
					jj<-which(DD$Sex%in%sex[j])
					CD[[i]][,j]<-if(cumulative) cumsum(DD[intersect(ii,jj),"COVID.19.Deaths"][mmwr]) else
						DD[intersect(ii,jj),"COVID.19.Deaths"][mmwr]
					TD[[i]][,j]<-if(cumulative) cumsum(DD[intersect(ii,jj),"Total.Deaths"][mmwr]) else
						DD[intersect(ii,jj),"Total.Deaths"][mmwr]
				}
			}
			if(cumulative){ 
				td<-cumsum(td)
				cd<-cumsum(cd)
			}
			if(show=="per.capita"){
				ii<-which(Age.Pop$Age.Group%in%age.group)
				Age.Pop<-Age.Pop[ii,]
				ii<-which(Age.Pop$Sex%in%sex)
				Age.Pop<-Age.Pop[ii,]
				pop<-sum(Age.Pop$Total.Population)
				cd<-cd/(pop/1e6)
				td<-td/(pop/1e6)
				ii<-which(Age.Pop$Age.Group%in%age.group)
				jj<-which(Age.Pop$Sex%in%sex)
				Pop<-Age.Pop[intersect(ii,jj),]
				for(i in 1:length(age.group)){
					for(j in 1:length(sex)){
						ii<-intersect(which(Age.Pop$Age.Group==age.group[i]),
							which(Age.Pop$Sex==sex[j]))
						CD[[i]][,j]<-CD[[i]][,j]/(Pop[ii,"Total.Population"]/1e6)
					}
				}
			} else if(show=="percent"){
				Tot<-td+cd
				cd<-cd/Tot*100
				td<-td/Tot*100
				for(i in 1:length(age.group)){
					for(j in 1:length(sex)){
						CD[[i]][,j]<-CD[[i]][,j]/Tot*100
					}
				}
			} else if(show=="percent.of.covid.deaths"){
				for(i in 1:length(age.group)){
					for(j in 1:length(sex)){
						CD[[i]][,j]<-CD[[i]][,j]/cd*100
						CD[[i]][is.nan(CD[[i]][,j]),j]<-0
					}
				}
				Tot<-td+cd
				cd<-cd/Tot*100
				td<-td/Tot*100
			}
			if(split.groups) 
				if(show%in%c("per.capita","percent")&&plot%in%c("bar","standard")) 
					plot<-"smooth"
			ms<-cumsum(c(0,31,29,31,30,31,30,31,31,30,31,30,31,
				31,28,31,30,31,30,31,31,30,31,30,31,
				31,28,31,30,31,30,31,31,30,31,30,31))
			mm<-c("Jan '20","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec",
				"Jan '21","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec",
				"Jan '22","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec",
				"Jan '23")
			xx<-seq(from=35.5,by=7,length.out=length(td))
			if(plot!=FALSE){
				if(show.total.deaths)
					par(mfrow=c(2,1),mar=c(5.1,5.1,3.1,2.1),bg=bg)
				else par(mar=c(5.1,5.1,3.1,2.1),bg=bg)
			}
			pp<-if(show=="per.capita") "/ 1M population" else 
				if(show=="percent") "as % of all deaths" else 
				if(show=="percent.of.covid.deaths") "as % of COVID deaths" else
				""	
			qq<-if(show=="per.capita") "/ 1M" else 
				if(show=="percent"||show=="percent.of.covid.deaths") "%" else ""
			aa<-age.group
			if(any(aa=="Under 1 year")) aa[aa=="Under 1 year"]<-"<1 years"
			if(any(aa=="85 years and over")) aa[aa=="85 years and over"]<-">85 years"
			aa<-gsub(" years","",aa)
			rr<-if(split.groups) "" else 
				paste("(",paste(aa,collapse=", "),")",sep="")
			if(show=="percent.of.covid.deaths"){
				tmp<-cd
				if(cumulative==FALSE){
					cd<-cd/cd*100
					cd[is.nan(cd)]<-0
				} else cd<-cd/max(cd)*100
			} else tmp<-cd
			if(is.null(ylim)){
				y.lim<-if(split.groups&&plot=="smooth")
					c(0,max(sapply(CD,max))) else c(0,1.2*max(cd))
				if(split.groups&&plot!="smooth"&&show=="percent.of.covid.deaths")
					y.lim<-c(0,120)
			}
			if(plot!=FALSE)
				plot(NA,xlim=xlim,ylim=y.lim,bty="n",axes=FALSE,
					xlab="",ylab="")
			if(palette=="new"){
				cols<-hcl.colors(22,palette="Temps")
				## cols<-cols[sort(c(seq(1,32,by=3),seq(2,32,by=3)))]
			} else if(palette=="original"){
				cols<-colorRampPalette(colors=brewer.pal("YlOrRd",n=8))(52)
				cols<-cols[sort(c(seq(1,52,by=5),seq(2,52,by=5)))]
			}
			Cols<-matrix(cols,11,2,byrow=TRUE,
				dimnames=list(c("Under 1 year","1-4 years","5-14 years",
				"15-24 years","25-34 years","35-44 years","45-54 years",
				"55-64 years","65-74 years","75-84 years","85 years and over"),
				c("Female","Male")))
			Cols<-Cols[age.group,sex,drop=FALSE]
			if(plot=="standard"){
				if(!split.groups){
					polygon(x=c(xx,max(xx),min(xx)),y=c(cd,0,0),border=FALSE,
						col=palette()[2])
				} else {
					tots<-rep(0,length(cd))
					for(i in length(age.group):1){
						for(j in length(sex):1){
							polygon(x=c(xx,xx[length(xx):1]),
								y=c(CD[[i]][,j]+tots,tots[length(tots):1]),
								border=FALSE,col=Cols[names(CD)[i],j])
							tots<-tots+CD[[i]][,j]
						}
					}
				}
			} else if(plot=="bar"){
				if(!split.groups){
					for(i in 1:length(cd)) polygon(xx[i]+c(-3.25,3.25,3.25,-3.25),
						c(0,0,cd[i],cd[i]),border=FALSE,col=palette()[2])
				} else {
					tots<-rep(0,length(cd))
					for(i in length(age.group):1){
						for(j in length(sex):1){
							for(k in 1:length(cd)){
								polygon(xx[k]+c(-3.25,3.25,3.25,-3.25),
									c(0,0,CD[[i]][k,j],CD[[i]][k,j])+tots[k],
									border=FALSE,col=Cols[names(CD)[i],j])
								tots[k]<-tots[k]+CD[[i]][k,j]
							}
						}
					}
				}
			} else if(plot=="smooth"){
				cds<-predict(loess(cd~xx,span=0.1))
				if(!split.groups){
					if(show!="percent.of.covid.deaths")
						lines(xx,cds,lwd=2,lty="dashed",col=palette()[2])
					else
						lines(xx,cd,lwd=2,lty="dashed",col=palette()[2])
				} else {
					for(i in length(age.group):1){
						for(j in length(sex):1){
							lines(xx,predict(loess(CD[[i]][,j]~xx,span=0.1)),
								lwd=2,lty="dashed",col=Cols[names(CD)[i],j])
						}
					}
				}
			}
			if(plot!=FALSE){
				Args<-list(...)
				Args$side<-2
				Args$labels<-FALSE
				if(show=="percent.of.covid.deaths") 
					Args$at<-seq(0,100,by=20)
				h<-do.call(axis,Args)
				Args$at<-h
				Args$labels<-relabel.axis(h)
				do.call(axis,Args)
				abline(h=h,col=grey(0.75),lwd=1,lty="dotted")
				Args$side<-1
				Args$at<-ms
				Args$labels<-mm
				v<-do.call(axis,Args)
				title(ylab=if(cumulative) paste("cumulative deaths",qq) else paste("weekly deaths",qq))
				if(show.total.deaths){
					if(cumulative)
						mtext(paste("a) cumulative COVID-19 deaths",pp,rr),adj=0,line=1,cex=1.2)
					else
						mtext(paste("a) weekly COVID-19 deaths",pp,rr),adj=0,line=1,cex=1.2)
				}
			}
			cd<-tmp
			if(plot!=FALSE){
				if(!split.groups){
					legend(x="topleft","confirmed COVID-19 deaths",
						pch=15,cex=0.9,
						col=palette()[2],
						pt.cex=1.5,bty="n",xpd=TRUE,
						xjust=0.5,yjust=1)
				} else {
					BG<-if(show=="percent.of.covid.deaths") 
						make.transparent("white",0.75) else "transparent"
					if(ncol(Cols)==2){
						i<-1
						xy<-legend(x="topleft",rep("",nrow(Cols)+1),
							pch=15,cex=0.8,
							col=c("transparent",Cols[,i]),pt.cex=1.5,xpd=TRUE,
							xjust=0.5,yjust=1,bg="transparent",box.col="transparent")
						text(xy$text$x[1]+1,xy$text$y[1],
							if(colnames(Cols)[i]=="Female") "F" else "M",cex=0.8,pos=2)
					} else i<-0
					xy<-legend(x=if(i==1) xy$rect$left+strwidth("W") else "topleft",
						y=if(i==1) xy$rect$top else NULL,
						legend=c("",rownames(Cols)),
						pch=15,cex=0.8,
						col=c("transparent",Cols[,i+1]),
						pt.cex=1.5,bg=BG,box.col="transparent")
					text(xy$text$x[1]+1,xy$text$y[1],if(colnames(Cols)[i+1]=="Female") "F" else "M",
						cex=0.8,pos=2)
					if(ncol(Cols)==2){
						i<-1
						xy<-legend(x="topleft",rep("",nrow(Cols)+1),
							pch=15,cex=0.8,
							col=c("transparent",Cols[,i]),pt.cex=1.5,xpd=TRUE,
							xjust=0.5,yjust=1,bg="transparent",box.col="transparent")
						text(xy$text$x[1]+1,xy$text$y[1],
							if(colnames(Cols)[i]=="Female") "F" else "M",cex=0.8,pos=2)
					}
				}
			}
			if(show.total.deaths){
				if(plot!=FALSE){
					if(is.null(ylim)) y.lim<-c(0,1.2*max(td+cd))
					plot(NA,xlim=xlim,ylim=y.lim,bty="n",axes=FALSE,
						xlab="",ylab="")
					if(plot=="standard"){
						polygon(x=c(xx,max(xx),min(xx)),y=c(td,0,0),border=FALSE,
							col=palette()[4])
						polygon(x=c(xx,xx[length(xx):1]),y=c(td+cd,td[length(td):1]),
							border=FALSE,col=palette()[2])	
					} else if(plot=="bar"){
						for(i in 1:length(td)) polygon(xx[i]+c(-3.25,3.25,3.25,-3.25),
							c(0,0,td[i],td[i]),border=FALSE,col=palette()[4])
						for(i in 1:length(cd)) polygon(xx[i]+c(-3.25,3.25,3.25,-3.25),
							c(0,0,cd[i],cd[i])+td[i],border=FALSE,col=palette()[2])
					} else if(plot=="smooth"){
						if(show!="percent.of.covid.deaths")
							lines(xx,cds,lwd=2,lty="dashed",col=palette()[2])
						else
							lines(xx,cd,lwd=2,lty="dashed",col=palette()[2])
						tds<-predict(loess(td~xx,span=0.1))
						lines(xx,tds,lwd=2,lty="dashed",col=palette()[4])
					}
					Args<-list(...)
					Args$side<-2
					Args$labels<-FALSE
					if(show=="percent"||show=="percent.of.covid.deaths") 
						Args$at<-seq(0,100,by=20)
					h<-do.call(axis,Args)
					Args$at<-h
					Args$labels<-relabel.axis(h)
					do.call(axis,Args)
					abline(h=h,col=grey(0.75),lwd=1,lty="dotted")
					Args$side<-1
					Args$at<-ms
					Args$labels<-mm
					v<-do.call(axis,Args)
					legend(x="topleft",c("non COVID-19 deaths",
						"confirmed COVID-19 deaths"),pch=15,cex=0.9,
						col=palette()[c(4,2)],
						pt.cex=1.5,bty="n",xpd=TRUE,
						xjust=0.5,yjust=1)
					title(ylab=if(cumulative) paste("cumulative deaths",qq) else 
						paste("weekly deaths",qq))
					pp<-if(show=="percent.of.covid.deaths") "as % of all deaths" else pp
					if(cumulative)
						mtext(paste("b) cumulative COVID-19 and non-COVID deaths",pp,rr),adj=0,line=1,cex=1.2)
					else
						mtext(paste("b) weekly COVID-19 and non-COVID deaths",pp,rr),adj=0,line=1,cex=1.2)
				}
			}
		}
		invisible(list(CD=CD,TD=TD))
	}
}

shade<-function(rgb,deviation=10){
	new<-rgb+sample(c(-1,1),3,replace=TRUE)*deviation
	while(any(new<0)||any(new>255))
		new<-rgb+sample(c(-1,1),3,replace=TRUE)*deviation
	new
}
liamrevell/covid19.Explorer documentation built on Feb. 19, 2023, 3:09 p.m.