# -------------------------------------------------------------------------------------------- #
# R SCRIPT FOR SARS-CoV-2 INFECTION FATALITY RISK ESTIMATES AND CONFIDENCE INTERVALS IN SPAIN #
# #
# This code reproduces the analyses performed in Pastor-Barriuso R, Pérez-Gómez B, Hernán MA, #
# et al. Infection fatality risk for SARS-CoV-2 in community dwelling population of Spain: #
# nationwide seroepidemiological study. BMJ 2020;371:m4509. #
# #
# The supplementary methods in the web appendix of this paper provide further details on the #
# calculation of point and interval estimates of infection fatality risks and standardized #
# infection fatality risk ratios. #
# #
# -------------------------------------------------------------------------------------------- #
# --------------------------------------------------------------------------------------------
# Non-institutionalized population by sex and 10-year age group (0-9,..., 70-79, >=80 years)
# as of July 15, 2020 (Spanish National Institute of Statistics)
data <- matrix(c(0,0,46887075, # Overall
1,0,23006948, # All men
1,1, 2205512,
1,2, 2557872,
1,3, 2479132,
1,4, 2978715,
1,5, 3916706,
1,6, 3493827,
1,7, 2598212,
1,8, 1783664,
1,9, 993308,
6,0,23880127, # All women
6,1, 2078268,
6,2, 2396681,
6,3, 2404106,
6,4, 3012421,
6,5, 3877832,
6,6, 3563505,
6,7, 2803368,
6,8, 2138101,
6,9, 1605845),
ncol=3,byrow=T,dimnames=list(NULL,c("sex","age.cat","pop")))
# --------------------------------------------------------------------------------------------
# Seroprevalence of SARS-CoV-2 with chemiluminiscent microparticle immunoassay by sex and age group (ENE-COVID study, April 27–June 22, 2020)
# Note: Seropositivity for SARS-CoV-2 with immunoassay in at least one round (61,098 participants)
# Note: Standard errors and 95% confidence intervals for SARS-CoV-2 seroprevalences accounted for the ENE-COVID complex survey design
prev <- matrix(c(0,0,0.0491785,0.0017181,0.0459153,0.0526607, # Overall
1,0,0.0480349,0.0020331,0.0442 ,0.0521844, # All men
1,1,0.032491 ,0.0085854,0.0192738,0.0542704,
1,2,0.0364694,0.0050465,0.0277637,0.0477708,
1,3,0.0575589,0.0060172,0.0468298,0.0705642,
1,4,0.0469104,0.0048414,0.0382763,0.057376 ,
1,5,0.0533675,0.0040497,0.0459585,0.0618934,
1,6,0.0526599,0.0041078,0.0451596,0.0613259,
1,7,0.0489133,0.0046771,0.0405125,0.058949 ,
1,8,0.0468757,0.0058305,0.0366764,0.0597354,
1,9,0.0459045,0.0083506,0.0320389,0.0653654,
6,0,0.0502719,0.0019933,0.0465025,0.0543296, # All women
6,1,0.0423326,0.0099907,0.0265351,0.0668889,
6,2,0.0438495,0.0055863,0.0341074,0.0562123,
6,3,0.0571521,0.0061229,0.0462607,0.0704184,
6,4,0.0520345,0.0045337,0.0438257,0.0616816,
6,5,0.0533281,0.0040707,0.0458839,0.0619017,
6,6,0.0517353,0.0039205,0.0445629,0.0599897,
6,7,0.0500889,0.0045885,0.0418162,0.059896 ,
6,8,0.0462329,0.0052633,0.0369377,0.0577271,
6,9,0.0499314,0.0078711,0.036566 ,0.0678381),
ncol=6,byrow=T,dimnames=list(NULL,c("sex","age.cat","prev","prev.se","prev.ll","prev.ul")))
# Merge population with seroprevalence data
data <- merge(data,prev,by=c("sex","age.cat"))
rm(prev)
# --------------------------------------------------------------------------------------------
# Sensitivity analysis: Correction of estimated SARS-CoV-2 seroprevalences for the immunoassay sensitivity and specificity
# Note: Pooled estimates of sensitivity and specificity of the immunoassay (SARS-CoV-2 IgG for use with ARCHITECT, Abbott Laboratories)
# against RT-PCR were obtained from a bivariate random-effects meta-analysis of 23 studies of diagnostic test accuracy
# Logit-transformed pooled estimates of sensitivity and specificity, together with their covariance matrix
logit.accur <- c(2.266413,4.897884)
names(logit.accur) <- c("sens","spec")
var.logit.accur <- matrix(c(0.018066444,0.004085381,
0.004085381,0.026223817),nrow=2,dimnames=list(c("sens","spec"),c("sens","spec")))
# Back-transform the estimates and covariance matrix to the original scale
accur <- exp(logit.accur)/(1+exp(logit.accur))
mat <- diag(exp(logit.accur)/(1+exp(logit.accur))^2)
var.accur <- mat%*%var.logit.accur%*%mat
rm(logit.accur,var.logit.accur,mat)
# Corrected seroprevalences for sensitivity and specificity
data$corr.prev <- (data$prev+accur[2]-1)/(accur[1]+accur[2]-1)
# Standard errors of corrected seroprevalences (delta method)
mat <- cbind(-data$corr.prev,1-data$corr.prev)
data$corr.prev.se <- sqrt(data$prev.se^2+diag(mat%*%var.accur%*%t(mat)))/(accur[1]+accur[2]-1)
rm(accur,var.accur,mat)
# 95% confidence intervals for corrected seroprevalences (based on logit transform)
data$corr.prev.ll <- 1/(1+exp(-(log(data$corr.prev/(1-data$corr.prev))-qt(0.975,1500-183)*data$corr.prev.se/(data$corr.prev*(1-data$corr.prev)))))
data$corr.prev.ul <- 1/(1+exp(-(log(data$corr.prev/(1-data$corr.prev))+qt(0.975,1500-183)*data$corr.prev.se/(data$corr.prev*(1-data$corr.prev)))))
# --------------------------------------------------------------------------------------------
# Estimated number of SARS-CoV-2 infections in non-institutionalized population by sex and age group
data$infect <- data$pop*data$prev
data$infect.se <- data$pop*data$prev.se
data$infect.ll <- data$pop*data$prev.ll
data$infect.ul <- data$pop*data$prev.ul
cbind(format(data$pop/1000,digits=1,nsmall=1,trim=T,big.mark=","), # Population (thousands)
paste(format(100*data$prev,digits=0,nsmall=1,trim=T)," (", # Seroprevalence (%) (95% CI)
format(100*data$prev.ll,digits=0,nsmall=1,trim=T)," to ",
format(100*data$prev.ul,digits=0,nsmall=1,trim=T),")",sep=""),
paste(format(data$infect/1000,digits=0,nsmall=1,trim=T,big.mark=",")," (", # Infections (thousands) (95% CI)
format(data$infect.ll/1000,digits=0,nsmall=1,trim=T,big.mark=",")," to ",
format(data$infect.ul/1000,digits=0,nsmall=1,trim=T,big.mark=","),")",sep=""))
data$corr.infect <- data$pop*data$corr.prev
data$corr.infect.se <- data$pop*data$corr.prev.se
data$corr.infect.ll <- data$pop*data$corr.prev.ll
data$corr.infect.ul <- data$pop*data$corr.prev.ul
cbind(format(data$pop/1000,digits=1,nsmall=1,trim=T,big.mark=","),
paste(format(100*data$corr.prev,digits=0,nsmall=1,trim=T)," (",
format(100*data$corr.prev.ll,digits=0,nsmall=1,trim=T)," to ",
format(100*data$corr.prev.ul,digits=0,nsmall=1,trim=T),")",sep=""),
paste(format(data$corr.infect/1000,digits=0,nsmall=1,trim=T,big.mark=",")," (",
format(data$corr.infect.ll/1000,digits=0,nsmall=1,trim=T,big.mark=",")," to ",
format(data$corr.infect.ul/1000,digits=0,nsmall=1,trim=T,big.mark=","),")",sep=""))
# --------------------------------------------------------------------------------------------
# Number of deaths with confirmed COVID-19 in non-institutionalized population up to July 15, 2020 by sex and age group (SiViES)
death <- matrix(c(0,0,19228, # Overall
1,0,12317, # All men
1,1,3,
1,2,3,
1,3,18,
1,4,48,
1,5,192,
1,6,705,
1,7,1904,
1,8,4145,
1,9,5299,
6,0,6911, # All women
6,1,2,
6,2,3,
6,3,17,
6,4,29,
6,5,103,
6,6,318,
6,7,749,
6,8,1986,
6,9,3704),
ncol=3,byrow=T,dimnames=list(NULL,c("sex","age.cat","death.covid")))
# Merge death data
data <- merge(data,death,by=c("sex","age.cat"))
rm(death)
format(data$death.covid,digits=1,nsmall=0,trim=T,big.mark=",")
# --------------------------------------------------------------------------------------------
# Number of excess all-cause deaths in non-institutionalized population by sex and age group (March 13 - May 22, 2020, MoMo)
death <- matrix(c(0,0,24778, # Overall
1,0,15480, # All men
1,1,32,
1,2,0,
1,3,0,
1,4,3,
1,5,168,
1,6,601,
1,7,2065,
1,8,5114,
1,9,7497,
6,0,9298, # All women
6,1,11,
6,2,22,
6,3,10,
6,4,71,
6,5,91,
6,6,369,
6,7,875,
6,8,2646,
6,9,5203),
ncol=3,byrow=T,dimnames=list(NULL,c("sex","age.cat","death.excess")))
# Merge death data
data <- merge(data,death,by=c("sex","age.cat"))
rm(death)
format(data$death.excess,digits=1,nsmall=0,trim=T,big.mark=",")
# --------------------------------------------------------------------------------------------
# Infection fatality risks in non-institutionalized population by sex and age group
# Note: CIs based on delta methods accounting for both binomial variance in the number of deaths and estimated variance in the number of infections
# IFR based on estimated SARS-CoV-2 seroprevalence
ifr.covid <- data$death.covid/data$infect
logit.ifr.covid.var <- 1/data$death.covid + 1/(data$infect-data$death.covid) + data$infect.se^2/(data$infect-data$death.covid)^2
paste(format(100*ifr.covid,digits=0,nsmall=2,trim=T)," (",
format(100/(1+exp(-(log(ifr.covid/(1-ifr.covid))-qt(0.975,1500-183)*sqrt(logit.ifr.covid.var)))),digits=0,nsmall=2,trim=T)," to ",
format(100/(1+exp(-(log(ifr.covid/(1-ifr.covid))+qt(0.975,1500-183)*sqrt(logit.ifr.covid.var)))),digits=0,nsmall=2,trim=T),")",sep="")
death.excess <- data$death.excess + 0.5*(data$death.excess==0) # Continuity correction for IFR with zero deaths
infect <- data$infect + 1*(data$death.excess==0)
ifr.excess <- death.excess/infect
logit.ifr.excess.var <- 1/death.excess + 1/(infect-death.excess) + data$infect.se^2/(infect-death.excess)^2
rm(death.excess,infect)
paste(format(100*ifr.excess,digits=0,nsmall=2,trim=T)," (",
format(100/(1+exp(-(log(ifr.excess/(1-ifr.excess))-qt(0.975,1500-183)*sqrt(logit.ifr.excess.var)))),digits=0,nsmall=2,trim=T)," to ",
format(100/(1+exp(-(log(ifr.excess/(1-ifr.excess))+qt(0.975,1500-183)*sqrt(logit.ifr.excess.var)))),digits=0,nsmall=2,trim=T),")",sep="")
# IFR based on corrected prevalence of SARS-CoV-2 infection
corr.ifr.covid <- data$death.covid/data$corr.infect
logit.corr.ifr.covid.var <- 1/data$death.covid + 1/(data$corr.infect-data$death.covid) + data$corr.infect.se^2/(data$corr.infect-data$death.covid)^2
paste(format(100*corr.ifr.covid,digits=0,nsmall=2,trim=T)," (",
format(100/(1+exp(-(log(corr.ifr.covid/(1-corr.ifr.covid))-qt(0.975,1500-183)*sqrt(logit.corr.ifr.covid.var)))),digits=0,nsmall=2,trim=T)," to ",
format(100/(1+exp(-(log(corr.ifr.covid/(1-corr.ifr.covid))+qt(0.975,1500-183)*sqrt(logit.corr.ifr.covid.var)))),digits=0,nsmall=2,trim=T),")",sep="")
death.excess <- data$death.excess + 0.5*(data$death.excess==0) # Continuity correction for IFR with zero deaths
corr.infect <- data$corr.infect + 1*(data$death.excess==0)
corr.ifr.excess <- death.excess/corr.infect
logit.corr.ifr.excess.var <- 1/death.excess + 1/(corr.infect-death.excess) + data$corr.infect.se^2/(corr.infect-death.excess)^2
rm(death.excess,corr.infect)
paste(format(100*corr.ifr.excess,digits=0,nsmall=2,trim=T)," (",
format(100/(1+exp(-(log(corr.ifr.excess/(1-corr.ifr.excess))-qt(0.975,1500-183)*sqrt(logit.corr.ifr.excess.var)))),digits=0,nsmall=2,trim=T)," to ",
format(100/(1+exp(-(log(corr.ifr.excess/(1-corr.ifr.excess))+qt(0.975,1500-183)*sqrt(logit.corr.ifr.excess.var)))),digits=0,nsmall=2,trim=T),")",sep="")
rm(corr.ifr.covid,logit.corr.ifr.covid.var,corr.ifr.excess,logit.corr.ifr.excess.var)
# --------------------------------------------------------------------------------------------
# Figure 2: Infection fatality risks based on estimated SARS-CoV-2 seroprevalence
age <- c(seq(5,75,10),90)
ifr.covid.m <- ifr.covid[3:11]
logit.ifr.covid.m.var <- logit.ifr.covid.var[3:11]
ifr.covid.w <- ifr.covid[13:21]
logit.ifr.covid.w.var <- logit.ifr.covid.var[13:21]
ifr.excess.m <- ifr.excess[3:11]
logit.ifr.excess.m.var <- logit.ifr.excess.var[3:11]
ifr.excess.w <- ifr.excess[13:21]
logit.ifr.excess.w.var <- logit.ifr.excess.var[13:21]
pdf(file="Figure 2.pdf",width=8.5/2.54,height=12/2.54,pointsize=6)
par(mfrow=c(2,1),mar=c(3.25,3.75,1,0),oma=c(0,0,0.35,0),tcl=-0.35)
plot(age,ifr.covid.m,type="n",xlim=c(0,90),ylim=c(0,0.22),xlab="",ylab="",axes=F)
polygon(c(age,rev(age)),
c( 1/(1+exp(-(log(ifr.covid.m/(1-ifr.covid.m))-qt(0.975,1500-183)*sqrt(logit.ifr.covid.m.var)))),
rev(1/(1+exp(-(log(ifr.covid.m/(1-ifr.covid.m))+qt(0.975,1500-183)*sqrt(logit.ifr.covid.m.var)))))),
density=-1,col=adjustcolor("deepskyblue2",alpha.f=0.25),border=NA)
polygon(c(age,rev(age)),
c( 1/(1+exp(-(log(ifr.covid.w/(1-ifr.covid.w))-qt(0.975,1500-183)*sqrt(logit.ifr.covid.w.var)))),
rev(1/(1+exp(-(log(ifr.covid.w/(1-ifr.covid.w))+qt(0.975,1500-183)*sqrt(logit.ifr.covid.w.var)))))),
density=-1,col=adjustcolor("red3",alpha.f=0.25),border=NA)
lines(age,ifr.covid.m,lty=1,lwd=2,col="deepskyblue2")
lines(age,ifr.covid.w,lty=1,lwd=2,col="red3")
axis(1,at=seq(0,90,10),labels=seq(0,90,10),cex.axis=1,font.axis=1,mgp=c(3,0.50,0),lty=1,lwd=1)
axis(2,at=seq(0,0.22,0.02),labels=seq(0,22,2),cex.axis=1,font.axis=1,mgp=c(3,0.65,0),las=1,adj=1,lty=1,lwd=1)
mtext(side=1,"Age (years)",cex=1,font=1,line=1.75)
mtext(side=2,"Infection fatality risk for SARS-CoV-2 (%)",cex=1,font=1,las=3,line=2.25)
mtext(side=2,at=0.235,"A",cex=1.25,font=1,las=1,line=2.35)
mtext(side=2,at=0.235,"Deaths with confirmed COVID-19",cex=1,font=1,las=1,adj=0,line=0)
legend(0,0.22,inset=0,bty="n",legend=c("Men","Women"),
lty=1,lwd=2,col=c("deepskyblue2","red3"),cex=1,y.intersp=1.25)
plot(age,ifr.excess.m,type="n",xlim=c(0,90),ylim=c(0,0.22),xlab="",ylab="",axes=F)
polygon(c(age,rev(age)),
c( 1/(1+exp(-(log(ifr.excess.m/(1-ifr.excess.m))-qt(0.975,1500-183)*sqrt(logit.ifr.excess.m.var)))),
rev(1/(1+exp(-(log(ifr.excess.m/(1-ifr.excess.m))+qt(0.975,1500-183)*sqrt(logit.ifr.excess.m.var)))))),
density=-1,col=adjustcolor("deepskyblue2",alpha.f=0.25),border=NA)
polygon(c(age,rev(age)),
c( 1/(1+exp(-(log(ifr.excess.w/(1-ifr.excess.w))-qt(0.975,1500-183)*sqrt(logit.ifr.excess.w.var)))),
rev(1/(1+exp(-(log(ifr.excess.w/(1-ifr.excess.w))+qt(0.975,1500-183)*sqrt(logit.ifr.excess.w.var)))))),
density=-1,col=adjustcolor("red3",alpha.f=0.25),border=NA)
lines(age,ifr.excess.m,lty=1,lwd=2,col="deepskyblue2")
lines(age,ifr.excess.w,lty=1,lwd=2,col="red3")
axis(1,at=seq(0,90,10),labels=seq(0,90,10),cex.axis=1,font.axis=1,mgp=c(3,0.50,0),lty=1,lwd=1)
axis(2,at=seq(0,0.22,0.02),labels=seq(0,22,2),cex.axis=1,font.axis=1,mgp=c(3,0.65,0),las=1,adj=1,lty=1,lwd=1)
mtext(side=1,"Age (years)",cex=1,font=1,line=1.75)
mtext(side=2,"Infection fatality risk for SARS-CoV-2 (%)",cex=1,font=1,las=3,line=2.25)
mtext(side=2,at=0.235,"B",cex=1.25,font=1,las=1,line=2.35)
mtext(side=2,at=0.235,"Excess deaths from all causes",cex=1,font=1,las=1,adj=0,line=0)
dev.off()
rm(age,ifr.covid.m,logit.ifr.covid.m.var,ifr.covid.w,logit.ifr.covid.w.var,ifr.excess.m,logit.ifr.excess.m.var,ifr.excess.w,logit.ifr.excess.w.var)
# --------------------------------------------------------------------------------------------
# Sex- and age-standardized infection fatality risk to the entire non-institutionalized population
# Infection fatality risks by sex and age group
ifr.covid <- ifr.covid[c(3:11,13:21)]
logit.ifr.covid.var <- logit.ifr.covid.var[c(3:11,13:21)]
ifr.covid.var <- ifr.covid^2*(1-ifr.covid)^2*logit.ifr.covid.var
names(ifr.covid) <- names(ifr.covid.var) <- c(paste("men.age",1:9,sep=""),paste("women.age",1:9,sep=""))
rm(logit.ifr.covid.var)
ifr.excess <- ifr.excess[c(3:11,13:21)]
logit.ifr.excess.var <- logit.ifr.excess.var[c(3:11,13:21)]
ifr.excess.var <- ifr.excess^2*(1-ifr.excess)^2*logit.ifr.excess.var
names(ifr.excess) <- names(ifr.excess.var) <- c(paste("men.age",1:9,sep=""),paste("women.age",1:9,sep=""))
rm(logit.ifr.excess.var)
# Sex and age distribution of the entire non-institutionalized population
weight <- data$pop[c(3:11,13:21)]/sum(data$pop[c(3:11,13:21)])
# Standardized infection fatality risk
std.ifr.covid <- sum(weight*ifr.covid)
std.ifr.covid.var <- sum(weight^2*ifr.covid.var)
paste(format(100*std.ifr.covid,digits=0,nsmall=2,trim=T)," (",
format(100/(1+exp(-(log(std.ifr.covid/(1-std.ifr.covid))-qt(0.975,1500-183)*sqrt(std.ifr.covid.var)/(std.ifr.covid*(1-std.ifr.covid))))),
digits=0,nsmall=2,trim=T)," to ",
format(100/(1+exp(-(log(std.ifr.covid/(1-std.ifr.covid))+qt(0.975,1500-183)*sqrt(std.ifr.covid.var)/(std.ifr.covid*(1-std.ifr.covid))))),
digits=0,nsmall=2,trim=T),")",sep="")
std.ifr.excess <- sum(weight*ifr.excess)
std.ifr.excess.var <- sum(weight^2*ifr.excess.var)
paste(format(100*std.ifr.excess,digits=0,nsmall=2,trim=T)," (",
format(100/(1+exp(-(log(std.ifr.excess/(1-std.ifr.excess))-qt(0.975,1500-183)*sqrt(std.ifr.excess.var)/(std.ifr.excess*(1-std.ifr.excess))))),
digits=0,nsmall=2,trim=T)," to ",
format(100/(1+exp(-(log(std.ifr.excess/(1-std.ifr.excess))+qt(0.975,1500-183)*sqrt(std.ifr.excess.var)/(std.ifr.excess*(1-std.ifr.excess))))),
digits=0,nsmall=2,trim=T),")",sep="")
rm(weight,std.ifr.covid,std.ifr.covid.var,std.ifr.excess,std.ifr.excess.var)
# --------------------------------------------------------------------------------------------
# Supplementary Figure 2: Sex- and age-standardized infection fatality risk ratios for each geographical region compared with the entire country
# Note: The distribution of SARS-CoV-2 infections in each geographical region was used as the standard
# Non-institutionalized population by geographical region, sex, and 10-year age group (July 15, 2020, Spanish National Institute of Statistics)
pop <- matrix(c(160037,189548,182954,245670, 350343,322458,275689,202845,124855,149704,177763,175717,245989, 353456, 336640,301253,246144,214043,
204927,234087,214489,257145, 363386,338015,266666,187440,107257,193808,219228,209587,259448, 356136, 340464,283803,222321,176866,
334698,365813,360468,443609, 566727,470271,325263,227048,121302,316744,346240,368903,466977, 589093, 509615,385267,291298,206450,
234322,277158,277359,332787, 431128,431445,337080,225195,146897,220585,260530,262176,317939, 414166, 422791,330793,256273,225301,
668280,775215,745336,878899,1172561,999300,739262,520055,278445,627606,719454,717355,895926,1151293,1011581,809626,623625,440159,
514572,599242,571772,664357, 832856,749957,533956,344717,177576,485569,563560,542107,666193, 818523, 761412,566625,411557,287402,
88676,116809,126754,156248, 199705,182381,120296, 76364, 36976, 84252,109906,128261,159949, 195165, 181002,126001, 86883, 55624),
ncol=18,byrow=T,dimnames=list(c("northwest","northeast","Madrid","central","east","south","Canary"),
c(paste("men.age",1:9,sep=""),paste("women.age",1:9,sep=""))))
# Seroprevalence of SARS-CoV-2 with immunoassay by geographical region, sex, and age group (ENE-COVID study, April 27–June 22, 2020)
prev.region <- c(0.015859,0.0386743,0.1150137,0.0870589,0.0422418,0.0198615,0.0141864) # Overall seroprevalence by geographical region
names(prev.region) <- c("northwest","northeast","Madrid","central","east","south","Canary")
prev <- matrix(c(0.0116323,0.0070508,0.0202856,0.0178784,0.0106334,0.0131712,0.0217481,0.0104676,0.0235730,0.0000000,0.0261002,0.0091958,0.0199994,0.0188421,0.0231923,0.0073796,0.0180255,0.0158169,
0.0107299,0.0240349,0.0416992,0.0401015,0.0530551,0.0466910,0.0379835,0.0227586,0.0480078,0.0218220,0.0434575,0.0332139,0.0410294,0.0412433,0.0447347,0.0447346,0.0438495,0.0102263,
0.0000000,0.0528818,0.1437742,0.0951984,0.1228826,0.1469668,0.1163952,0.1430424,0.1038248,0.1317784,0.1144435,0.1448796,0.1214396,0.1154044,0.1122389,0.1055483,0.1145744,0.1288093,
0.0575273,0.1025736,0.0856711,0.0722530,0.0871421,0.0916653,0.0704588,0.0899935,0.0582883,0.0963974,0.0795505,0.0940754,0.0973675,0.1031803,0.1007662,0.0866384,0.0806763,0.0716813,
0.0653336,0.0386236,0.0442795,0.0466920,0.0458188,0.0335747,0.0492345,0.0350766,0.0377887,0.0253335,0.0213040,0.0537833,0.0413752,0.0463914,0.0406075,0.0490818,0.0362487,0.0642789,
0.0091484,0.0142962,0.0255706,0.0221848,0.0255576,0.0235305,0.0188744,0.0104829,0.0215379,0.0271066,0.0227767,0.0135637,0.0214193,0.0190164,0.0223774,0.0206512,0.0126024,0.0136301,
0.0642236,0.0039325,0.0350148,0.0205724,0.0097205,0.0093357,0.0193720,0.0000000,0.0490143,0.0000000,0.0063997,0.0059447,0.0171654,0.0177938,0.0154786,0.0200427,0.0000000,0.0000000),
ncol=18,byrow=T,dimnames=list(c("northwest","northeast","Madrid","central","east","south","Canary"),
c(paste("men.age",1:9,sep=""),paste("women.age",1:9,sep=""))))
prev.se <- matrix(c(0.0116642,0.0049889,0.0094242,0.0071710,0.0035438,0.0047243,0.0068321,0.0053291,0.0142583,0.0000000,0.0094828,0.0047892,0.0086455,0.0054200,0.0058487,0.0031312,0.0071460,0.0096394,
0.0088660,0.0085721,0.0129972,0.0105181,0.0088907,0.0081669,0.0093992,0.0067129,0.0154903,0.0165767,0.0124666,0.0109274,0.0101349,0.0086040,0.0069975,0.0100971,0.0105668,0.0040388,
0.0000000,0.0193462,0.0250401,0.0197442,0.0159828,0.0218566,0.0250514,0.0297866,0.0442606,0.0703434,0.0258742,0.0250464,0.0207181,0.0174886,0.0187184,0.0200884,0.0261549,0.0480963,
0.0195881,0.0234049,0.0130497,0.0107801,0.0091349,0.0102049,0.0087564,0.0135757,0.0177023,0.0292757,0.0147122,0.0151489,0.0123934,0.0104522,0.0102818,0.0099819,0.0129207,0.0159895,
0.0266228,0.0095426,0.0106762,0.0105564,0.0081929,0.0064074,0.0091210,0.0105903,0.0155433,0.0157792,0.0066435,0.0121481,0.0072974,0.0075712,0.0070937,0.0098174,0.0089616,0.0154839,
0.0065717,0.0057070,0.0072875,0.0058878,0.0063038,0.0057467,0.0048379,0.0055460,0.0103896,0.0124651,0.0079414,0.0044557,0.0059636,0.0042647,0.0044520,0.0056149,0.0052939,0.0078844,
0.0583800,0.0039479,0.0221202,0.0197326,0.0041727,0.0052434,0.0098695,0.0000000,0.0463344,0.0000000,0.0063752,0.0058667,0.0090457,0.0091558,0.0080209,0.0122432,0.0000000,0.0000000),
ncol=18,byrow=T,dimnames=list(c("northwest","northeast","Madrid","central","east","south","Canary"),
c(paste("men.age",1:9,sep=""),paste("women.age",1:9,sep=""))))
df <- c(155,193,96,279,276,258,60) # Survey design degrees of freedom
names(df) <- c("northwest","northeast","Madrid","central","east","south","Canary")
# Estimated number of SARS-CoV-2 infections in non-institutionalized population by geographical region, sex, and age group
infect <- pop*prev
infect.se <- pop*prev.se
rm(pop,prev,prev.se)
# Number of deaths with confirmed COVID-19 in non-institutionalized population up to July 15, 2020 by geographical region (SiViES)
death.covid <- c(519,1518,7469,3180,5427,972,143)
names(death.covid) <- c("northwest","northeast","Madrid","central","east","south","Canary")
# Number of excess all-cause deaths in non-institutionalized population by geographical region (March 13 - May 22, 2020, MoMo)
death.excess <- c(526,1631,8164,4248,9048,1054,107)
names(death.excess) <- c("northwest","northeast","Madrid","central","east","south","Canary")
# Standardized infection fatality risk ratios for each geographical region compared with the entire country
# Note: CIs based on delta methods accounting for all relevant sources of uncertainty
expdth.covid <- (infect%*%ifr.covid)[,1]
expdth.covid.var <- (infect.se^2%*%ifr.covid^2 + infect^2%*%ifr.covid.var)[,1]
srr.covid <- death.covid/expdth.covid
log.srr.covid.var <- 1/death.covid - 1/apply(infect,1,sum) + expdth.covid.var/expdth.covid^2
rm(ifr.covid,ifr.covid.var,expdth.covid.var)
expdth.excess <- (infect%*%ifr.excess)[,1]
expdth.excess.var <- (infect.se^2%*%ifr.excess^2 + infect^2%*%ifr.excess.var)[,1]
srr.excess <- death.excess/expdth.excess
log.srr.excess.var <- 1/death.excess - 1/apply(infect,1,sum) + expdth.excess.var/expdth.excess^2
rm(infect,infect.se,ifr.excess,ifr.excess.var,expdth.excess.var)
# Supplementary Figure 2
region <- c("Northwest","Northeast","Madrid","Central","East","South","Canary Islands")
pdf(file="Supplementary Figure 2.pdf",width=18.5/2.54,height=5.5/2.54,pointsize=6)
par(mfrow=c(1,2),mar=c(1.85,14.25,5.10,0),oma=c(0,14.5,0,0.25),tcl=-0.35)
plot(log(srr.covid),7:1,type="n",xlim=log(c(0.25,2)),ylim=c(0,7),xlab="",ylab="",axes=F)
symbols(log(srr.covid),7:1,squares=1/sqrt(log.srr.covid.var),
inches=0.085,add=T,lty=1,lwd=1,fg="black",bg="black")
segments((log(srr.covid)-qt(0.975,df)*sqrt(log.srr.covid.var))[-7],(7:1)[-7],
(log(srr.covid)+qt(0.975,df)*sqrt(log.srr.covid.var))[-7],(7:1)[-7],lty=1,lwd=1)
arrows(log(0.25),(7:1)[7],
log(2),(7:1)[7],code=3,length=0.035,lty=1,lwd=1)
abline(v=0,lty=2,lwd=1)
axis(1,at=log(c(0.25,0.5,1,2)),labels=c(0.25,0.5,1,2),cex.axis=1,font.axis=1,mgp=c(3,0.50,0),lty=1,lwd=1)
mtext(side=2,at=7:1,region,cex=1,font=1,las=1,adj=0,line=28.25)
mtext(side=2,at=8,expression(bold("Region")),cex=1,font=1,las=1,adj=0,line=28.25)
mtext(side=2,at=7:1,format(100*prev.region,digits=0,nsmall=1,trim=T),cex=1,font=1,las=1,adj=1,line=17.75)
mtext(side=2,at=8.5,expression(bold("SARS-CoV-2")),cex=1,font=1,las=1,adj=0.5,line=18.25)
mtext(side=2,at=8,expression(bold("seroprevalence (%)")),cex=1,font=1,las=1,adj=0.5,line=18.25)
mtext(side=2,at=7:1,paste(format(death.covid,trim=T,big.mark=","),"/",sep=""),cex=1,font=1,las=1,adj=1,line=10.75)
mtext(side=2,at=7:1,format(expdth.covid,digits=0,nsmall=0,trim=T,big.mark=","),cex=1,font=1,las=1,adj=0,line=10.75)
mtext(side=2,at=8.5,expression(bold("Observed/")),cex=1,font=1,las=1,adj=0.5,line=10.75)
mtext(side=2,at=8,expression(bold("expected")),cex=1,font=1,las=1,adj=0.5,line=10.75)
mtext(side=2,at=7:1,paste(format(srr.covid,digits=0,nsmall=2)," (",
format(exp(log(srr.covid)-qt(0.975,df)*sqrt(log.srr.covid.var)),digits=0,nsmall=2)," to ",
format(exp(log(srr.covid)+qt(0.975,df)*sqrt(log.srr.covid.var)),digits=0,nsmall=2),")",sep=""),
cex=1,font=1,las=1,adj=1,line=0.5)
mtext(side=2,at=8.5,expression(bold("Standardized infection fatality risk ratio")),cex=1,font=1,las=1,adj=0.5,line=-3.5)
mtext(side=2,at=8,expression(bold("for SARS-CoV-2 (95% confidence interval)")),cex=1,font=1,las=1,adj=0.5,line=-3.5)
mtext(side=2,at=9.5,expression(bold("Deaths with confirmed COVID-19")),cex=1,font=1,las=1,adj=0.5,line=0.75)
plot(log(srr.excess),7:1,type="n",xlim=log(c(0.25,2)),ylim=c(0,7),xlab="",ylab="",axes=F)
symbols(log(srr.excess),7:1,squares=1/sqrt(log.srr.excess.var),
inches=0.085*max(1/sqrt(log.srr.excess.var))/max(1/sqrt(log.srr.covid.var)),add=T,lty=1,lwd=1,fg="black",bg="black")
segments((log(srr.excess)-qt(0.975,df)*sqrt(log.srr.excess.var))[-7],(7:1)[-7],
(log(srr.excess)+qt(0.975,df)*sqrt(log.srr.excess.var))[-7],(7:1)[-7],lty=1,lwd=1)
arrows(log(0.25),(7:1)[7],
(log(srr.excess)+qt(0.975,df)*sqrt(log.srr.excess.var))[7],(7:1)[7],code=1,length=0.035,lty=1,lwd=1)
abline(v=0,lty=2,lwd=1)
axis(1,at=log(c(0.25,0.5,1,2)),labels=c(0.25,0.5,1,2),cex.axis=1,font.axis=1,mgp=c(3,0.50,0),lty=1,lwd=1)
mtext(side=2,at=7:1,paste(format(death.excess,trim=T,big.mark=","),"/",sep=""),cex=1,font=1,las=1,adj=1,line=10.75)
mtext(side=2,at=7:1,format(expdth.excess,digits=0,nsmall=0,trim=T,big.mark=","),cex=1,font=1,las=1,adj=0,line=10.75)
mtext(side=2,at=8.5,expression(bold("Observed/")),cex=1,font=1,las=1,adj=0.5,line=10.75)
mtext(side=2,at=8,expression(bold("expected")),cex=1,font=1,las=1,adj=0.5,line=10.75)
mtext(side=2,at=7:1,paste(format(srr.excess,digits=0,nsmall=2)," (",
format(exp(log(srr.excess)-qt(0.975,df)*sqrt(log.srr.excess.var)),digits=0,nsmall=2)," to ",
format(exp(log(srr.excess)+qt(0.975,df)*sqrt(log.srr.excess.var)),digits=0,nsmall=2),")",sep=""),
cex=1,font=1,las=1,adj=1,line=0.5)
mtext(side=2,at=8.5,expression(bold("Standardized infection fatality risk ratio")),cex=1,font=1,las=1,adj=0.5,line=-3.5)
mtext(side=2,at=8,expression(bold("for SARS-CoV-2 (95% confidence interval)")),cex=1,font=1,las=1,adj=0.5,line=-3.5)
mtext(side=2,at=9.5,expression(bold("Excess deaths from all causes")),cex=1,font=1,las=1,adj=0.5,line=0.75)
dev.off()