Links to earlier parts: czech.html, czech2.html, czech3.html.
Kirsch told me on Twitter: "You have all failed to explain why in every age group where we have sufficient data. Moderna was significantly more deadly than Pfizer. it wasn't comorbidities because those didn't matter. So why was the death rate higher? Why can't you answer that simple question?" [https://x.com/stkirsch/status/1818053265357160918]
I told him that in the plot below if you look at the age groups 50-59 and above where there's a sufficient sample size of deaths, the Moderna-Pfizer ratio was close to 1.0 among people who got the first dose in the 4th quarter of 2021. So did Moderna vaccines get less deadly in the 4th quarter?
library(data.table)
agecut=\(x,y)cut(x,c(y,Inf),paste0(y,c(paste0("-",y[-1]-1),"+")),T,F)
quar=\(x){u=unique(x);paste0(substr(u,1,4),"Q",(as.numeric(substr(u,6,7))-1)%/%3+1)[match(x,u)]}
ages=c(0,2:9*10)
t=fread("http://sars2.net/f/czbucketskeep.csv.gz")[dose<=1][,age:=pmin(age,95)]
t=merge(t[dose==1&vaxmonth>="2021-01"&vaxmonth<="2022-03"],t[,.(base=sum(dead)/sum(alive)),age])[,base:=base*alive]
t=t[,.(base=sum(base),dead=sum(dead),alive=sum(as.double(alive))),.(type,vaxmonth=quar(vaxmonth),age=agecut(age,ages))]
t=rbind(t,t[,.(base=sum(base),dead=sum(dead),alive=sum(alive),vaxmonth="Total"),.(age,type)])
t=rbind(t,t[,.(base=sum(base),dead=sum(dead),alive=sum(alive),age="Total"),.(vaxmonth,type)])
a=t[,tapply(dead/base,.(type,age,vaxmonth),c)]
m1=a["Moderna",,];m2=a["Pfizer",,]
m=(m1-m2)/ifelse(m1>m2,m2,m1)
pop=t[type%in%c("Moderna","Pfizer"),xtabs(alive~age+vaxmonth)]/365
disp=matrix(sprintf(ifelse(m1/m2>10,"%.1f","%.2f"),m1/m2),nrow(m))
hide=pmin(m1,m2)==0;m[hide]=NA;disp[is.na(m)]=""
exp=.7;m=abs(m)^exp*sign(m);maxcolor=10^exp
pheatmap::pheatmap(m,filename="i1.png",display_numbers=disp,
gaps_row=nrow(m)-1,gaps_col=ncol(m)-1,
cluster_rows=F,cluster_cols=F,legend=F,cellwidth=19,cellheight=19,fontsize=9,fontsize_number=8,
border_color=NA,na_col="gray90",
number_color=ifelse((abs(m)>.55*maxcolor)&!is.na(m),"white","black"),
breaks=seq(-maxcolor,maxcolor,,256),
colorRampPalette(hsv(rep(c(7/12,0),5:4),c(.9,.75,.6,.3,0,.3,.6,.75,.9),c(.4,.65,1,1,1,1,1,.65,.4)))(256))
rec=fread("Czech/data/CR_records.csv",showProgress=F)
pop=rec[grep("comirnaty|spikevax",ignore.case=T,OckovaciLatka_1),.N,.(age=agecut(2021-Rok_narozeni,ages),vaxq=quar(Datum_1))]
pop=rbind(pop,pop[,.(N=sum(N),age="Total"),vaxq])
pop=rbind(pop,pop[,.(N=sum(N),vaxq="Total"),age])
pop=xtabs(N~age+vaxq,pop)[rownames(m),colnames(m)]
exp2=.6;maxcolor2=max(pop[-nrow(pop),])
kimi=\(x){e=floor(log10(ifelse(x==0,1,abs(x))));e2=pmax(e,0)%/%3+1;p=!is.na(x)
x[p]=paste0(sprintf(paste0("%.",ifelse(e[p]%%3==0,1,0),"f"),x[p]/1e3^(e2[p]-1)),c("","k","M","B","T")[e2[p]]);x}
exp2=.6;maxcolor2=max(pop[-nrow(pop),-ncol(pop)])
pheatmap::pheatmap(pop^exp2,filename="i2.png",display_numbers=kimi(ifelse(pop==0,0,pop)),
gaps_row=nrow(m)-1,gaps_col=ncol(m)-1,
cluster_rows=F,cluster_cols=F,legend=F,cellwidth=19,cellheight=19,fontsize=9,fontsize_number=8,
border_color=NA,na_col="white",
number_color=ifelse(pop^exp2>maxcolor2^exp2*.45,"white","black"),
breaks=seq(0,maxcolor2^exp2,,256),
sapply(seq(1,0,,256),\(i)rgb(i,i,i)))
cap="Czech record-level data: Age-normalized Moderna-Pfizer mortality ratio by vaccination quarter and age group. Only first doses are included. The mortality rate was calculated from the day of the first dose up to the end of 2022. A baseline number of deaths was calculated for each single year of age by multiplying the person-days of the age by the total mortality rate for the age in 2020 to 2022 among all people who are included in the record-level data, and the baseline deaths for all ages within an age group were added together. The ratio was calculated as (Moderna_deaths / Moderna_baseline) / (Pfizer_deaths / Pfizer_baseline). The right side shows the number of people each quarter who received either a Pfizer or Moderna vaccine for their first dose."
system(paste0("magick i[12].png +append 0.png;w=`identify -format %w 0.png`;pad=24;convert -pointsize 42 -interline-spacing -3 -font Arial \\( -size $[w-pad*2]x caption:'",cap,"' -splice $[pad]x20 \\) 0.png -append 1.png"))
Kirsch has been saying that the Moderna-Pfizer ratio is consistently around 1.3 for all months of vaccination that have a sufficiently large sample size of vaccinated people (even though my age-normalized Moderna-Pfizer ratio was about 1.02 for people who got the first dose in November 2021, and there's about about 300,000 people who got a first dose from a Pfizer or Moderna vaccine in November 2021 so I think it's a sufficiently large sample size).
But Kirsch has generally been looking at first doses only. The plot below shows the results for the third dose instead. I used the same R script as in this section: czech2.html#Triangle_plot_of_Moderna_Pfizer_ratio.
The total Moderna-Pfizer ratio is below 1.0 for people who got the third dose during each month from March 2022 to August 2022, but it shouldn't just be random noise because there's a total of about 300,000 people who received a third shot from a Moderna or Pfizer vaccine during those months. So there appears to be a similar effect as with the first dose, where the late vaccinees have a lower Moderna-Pfizer ratio than the early vaccinees:
The monthly number of 4th doses peaked in October 2022. But for people who got a 4th dose in October 2022, the Moderna-Pfizer ratio was about 0.67:
I had forgotten to make these triangle plots for the third and fourth doses, but now I think these plots have become yet another nail in the coffin of Kirsch's theory that his Moderna-Pfizer ratio proved that Moderna vaccines were more deadly. Now of course it's possible that Moderna vaccines were actually more deadly, but I don't think Kirsch's Moderna-Pfizer ratio is the smoking gun piece of evidence that he makes it out to be.
In the file prehled-ockovacich-mist.csv which has information about vaccination sites, the sites are classified under these 5 types: [https://ockodoc.mzcr.cz/wp-content/uploads/2024/01/231230_manual_crs.pdf]
| Abbreviation | Translation | Czech |
|---|---|---|
| OČM | standard vaccination site | standardní očkovací místo |
| DOČM | vaccination site for children | očkovací místo pro děti |
| VOČM | large-capacity vaccination site | velkokapacitní očkovací místo |
| SMPL | vaccination site for self-payers | očkovací místo pro samoplátce |
| WALKIN | mobile vaccination site | mobilní očkovací místo |
The file ockovani-profese.csv contains about 19 million rows with one row for each administered dose, but when I tried joining the file with the file for information about vaccination sites, only about 7 million rows could be joined by the vaccinator ID. The joined rows didn't include any row for a large-capacity vaccination site. And I further excluded vaccination sites for children and I excluded vaccines given to ages below 18, because children have little effect on overall mortality but the pediatric vaccination sites distorted the overall percentages of site types. But this shows the distribution of the remaining 3 types of vaccination sites:
dlf=\(x,y,...){if(missing(y))y=sub(".*/","",x);for(i in 1:length(x))download.file(x[i],y[i],quiet=T,...)}
dlf("https://onemocneni-aktualne.mzcr.cz/api/v2/covid-19/ockovani-profese.csv") # about 4 GiB
dlf("https://onemocneni-aktualne.mzcr.cz/api/v2/covid-19/prehled-ockovacich-mist.csv") # about 20 kiB
pro=fread("ockovani-profese.csv")
mist=fread("prehled-ockovacich-mist.csv")
d=pro[,.N,.(vaxtype=vakcina,vaccinator=zarizeni_kod,age=vekova_skupina)]
d=d[!d$age%in%c("0-4","05-11","12-15","16-17")]
name1=unique(d$vaxtype);name2=rep("Other",length(name1))
name2[name1==""]=""
name2[grep("comirnaty",ignore.case=T,name1)]="Pfizer"
name2[grep("spikevax",ignore.case=T,name1)]="Moderna"
name2[grep("nuvaxovid",ignore.case=T,name1)]="Novavax"
name2[name1=="COVID-19 Vaccine Janssen"]="Janssen"
name2[name1=="VAXZEVRIA"]="AstraZeneca"
d$vaxtype=name2[match(d$vaxtype,name1)]
d=merge(unique(mist[,.(sitetype=ockovaci_misto_typ,vaccinator=nrpzs_kod)],by="vaccinator"),d)
m=d[sitetype!="DOČM",xtabs(N~vaxtype+sitetype)]
m=m[order(-rowSums(m)),]
round(m/rowSums(m)*100)
sitetype vaxtype OČM SMPL WALKIN # percentage of site types by manufacturer Pfizer 54 37 9 Moderna 66 20 14 Janssen 67 29 4 AstraZeneca 46 42 12 Novavax 65 32 3 Other 42 46 11
So Pfizer has a higher percentage of vaccination sites for self-payers than Moderna (SMPL), but Moderna has a higher percentage of standard vaccination sites (OČM) and walk-in sites. Which is yet another line of evidence that there are systematic biases in the distribution of vaccine doses by different manufacturers.
In this blog post where Hans-Joachim Kremer looked at US VAERS data, he got a higher ratio of deaths per adverse event reports for Pfizer than Moderna vaccines, so he pointed out that it conflicted with Kirsch's analysis of the Czech data: https://tkp.at/2024/08/06/spikevax-gefaehrlicher-als-comirnaty-steve-kirsch-greift-daneben/.
However it could be that Moderna had both a higher number of adverse event reports per dose than Pfizer and a higher number of deaths per dose than Pfizer, but the ratio of deaths per adverse event reports could've still been higher for Pfizer than Moderna. So he could've repeated the calculation so that he would've used doses administered instead of adverse event reports as the denominator.
I took the number of adverse event reports in Czechia per vaccine manufacturer from this spreadsheet: https://github.com/PalackyUniversity/batch-dependent-safety/blob/main/data/SUKL-batches-AE.xlsx. The spreadsheet only included the number of doses released but not the number of doses administered, so I calculated the number of doses administered from the record-level data. I only included doses administered up to July 4th 2023, which was the date when the authors who published the batch data wrote they received the data. But anyway, Moderna got a higher ratio than Pfizer for both deaths per doses administered and adverse event reports per doses administered:
rename=\(x){u=unique(x);o=rep("Other",length(u))
o[grep("comirnaty|pfizer",u,T)]="Pfizer"
o[grep("moderna|spikevax",u,T)]="Moderna"
o[grep("nuvaxovid|novavax",u,T)]="Novavax"
o[grep("janssen",u,T)]="Janssen"
o[grep("vaxzevria|astrazeneca",u,T)]="AstraZeneca"
o[match(x,u)]}
rec=fread("curl -Ls github.com/skirsch/Czech/raw/main/data/CR_records.csv.xz|xz -dc")
k=grep("Datum_",names(rec));k2=k+3
d=data.table(date=`class<-`(unlist(rec[,..k],,F),"Date"),type=rename(unlist(rec[,..k2],,F)))[!is.na(date)]
d=d[date<="2023-7-4",.(doses=.N),type]
f="SUKL-batches-AE.xlsx"
download.file("https://github.com/PalackyUniversity/batch-dependent-safety/raw/main/data/SUKL-batches-AE.xlsx",f)
s2=data.table(readxl::read_excel(f,sheet=2))
s2=s2[,.(reports=sum(no_reports),serious=sum(serious),dead=sum(deaths)),.(type=rename(name))]
o=merge(s2,d)
o[,deadpermildose:=dead/doses*1e6]
o[,reportspermildose:=reports/doses*1e6]
o[,deadpermilreport:=dead/reports*1e6]
dplyr::mutate_if(o[type!="Other"],is.double,round)[order(-doses)]|>print(r=F)
type reports serious dead doses deadpermildose reportspermildose deadpermilreport
Pfizer 10712 5183 134 15680944 9 683 12509
Moderna 1575 583 26 1643159 16 959 16508
AstraZeneca 1512 695 44 887601 50 1703 29101
Janssen 555 241 10 415077 24 1337 18018
Novavax 14 4 0 11190 0 1251 0
When I used Pfizer as the baseline, Moderna got about 1.9 times higher deaths per dose and about 1.4 times higher adverse event reports per dose, so the ratio of deaths per reports also ended up being higher for Moderna than Pfizer:
> round(o[type=="Moderna",-1]/o[type=="Pfizer",-1],3) reports serious dead doses deadpermildose reportspermildose deadpermilreport 1: 0.147 0.112 0.194 0.105 1.852 1.403 1.32
The spreadsheet didn't include information about age so I wasn't able to look at only specific age groups or normalize my calculation by age. In the record-level data the average age is about 7 years higher for people who got a Moderna vaccine for the first dose than people who got a Pfizer vaccine for the first dose. But earlier I got a correlation of only about 0.18 between the average age of people under a batch and the rate of deaths per batch that had been reported as adverse events: czech3.html#Batch_study_by_authors_from_Palack_University.
But anyway, next I tried looking at US VAERS data: https://vaers.hhs.gov/data/datasets.html. I got the number of doses administered per manufacturer in the United States from this JSON file: https://covid.cdc.gov/covid-data-tracker/COVIDData/getAjaxData?id=vaccination_data. It only included the cumulative number of doses administered up to May 10th 2023 which was when the data stopped being updated, when there were about 367 million doses administered for Pfizer and 232 million for Moderna. But I didn't find any source which showed the number of doses grouped by both vaccine manufacturer and age, so I wasn't able to adjust my calculation for age.
Moderna again got a higher rate of both deaths per dose and reports per dose than Pfizer. But unlike in the case of the Czech data, the Moderna-Pfizer ratio for deaths per dose was now lower than the ratio for reports per dose, so Pfizer got a higher ratio of deaths per reports:
names=read.csv(text="name,vaers,doses
Pfizer,PFIZER\\BIONTECH,366979906
Moderna,MODERNA,232147784
Janssen,JANSSEN,19007537
Novavax,NOVAVAX,89195
Unknown,UNKNOWN MANUFACTURER,890835")
# # this code is not needed because I added the number of doses from the API to the table above
# con=url("https://covid.cdc.gov/covid-data-tracker/COVIDData/getAjaxData?id=vaccination_data")
# j=jsonlite::fromJSON(readLines(con))$vaccination_data;close(con)
# names$doses=j[1,paste0("Administered_",ifelse(names$name=="Unknown","Unk_Manuf",names$name))]
va=do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSDATA.csv"))))
va=merge(va,do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSVAX.csv")))))
o=va[VAX_TYPE=="COVID19",.(reports=.N,dead=sum(DIED=="Y")),.(vaers=VAX_MANU)]
o=merge(o,names)[,.(name,reports,dead,doses)]
o[,deadpermildose:=dead/doses*1e6]
o[,reportspermildose:=reports/doses*1e6]
o[,deadpermilreport:=dead/reports*1e6]
o[order(-doses)]|>dplyr::mutate_if(is.double,round)|>print(r=F)
name reports dead doses deadpermildose reportspermildose deadpermilreport Pfizer 481402 10439 366979906 28 1312 21685 Moderna 469655 9821 232147784 42 2023 20911 Janssen 73454 1943 19007537 102 3864 26452 Unknown 7363 82 890835 92 8265 11137 Novavax 492 5 89195 56 5516 10163
This shows that the Moderna-Pfizer ratio is about 1.53 for reports per dose, 1.39 for deaths per dose, and 0.91 for deaths per report:
> round(o[name=="Moderna",-1]/o[name=="Pfizer",-1],3) reports dead doses deadpermildose reportspermildose deadpermilreport 1: 0.976 0.941 0.633 1.487 1.542 0.964
However ideally I should've adjusted these calculations for age. In VAERS the ratio of deaths per reports seems to approximately double between each 10-year age group from 30-39 up to 90+ (even though the ratio remains similar in ages 10-19, 20-29, and 30-39):
> va=do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSDATA.csv"))))
> va=merge(va,do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSVAX.csv")))))
> agecut=\(x,y)cut(x,c(y,Inf),paste0(y,c(paste0("-",y[-1]-1),"+")),T,F)
> o=va[!is.na(AGE_YRS)&VAX_TYPE=="COVID19",.(dead=sum(DIED=="Y"),reports=.N),.(age=agecut(AGE_YRS,0:9*10))]
> o[order(age)][,deadper10kreport:=round(dead/reports*1e4)][]|>print(r=F)
age dead reports deadper10kreport
0-9 20 16389 12
10-19 107 54171 20
20-29 201 86426 23
30-39 362 135659 27
40-49 672 140188 48
50-59 1572 155711 101
60-69 3786 162128 234
70-79 5910 117297 504
80-89 5445 48294 1127
90+ 2446 12197 2005
Kremer divided the VAERS data into only two groups for ages 18-64 and 65+, so I wanted to check if I would get different results for more fine-grained age groups. However when I calculated a Moderna-Pfizer ratio for deaths per reports, the ratio was below 1 in all age 5-year groups above 50-54, and some of the younger age groups where the ratio was above 1 had wide confidence intervals:
> o=va[!is.na(AGE_YRS)&VAX_TYPE=="COVID19",.(dead=sum(DIED=="Y"),reports=.N),.(age=agecut(AGE_YRS,0:19*5),type=VAX_MANU)][order(age)] > o1=o[type=="MODERNA"];o2=o[type=="PFIZER\\BIONTECH"] > rate=(o1$dead/o1$reports)/(o2$dead/o2$reports) > ci=exp(log(rate)+outer(qnorm(.975)*sqrt(1/o1$dead-1/o1$reports+1/o2$dead-1/o2$reports),c(-1,1))) > o1[,.(age,rate,low=ci[,1],high=ci[,2])]|>mutate_if(is.double,round,2)|>print(r=F) age rate low high 0-4 1.47 0.40 5.47 5-9 0.00 0.00 NaN 10-14 2.31 0.79 6.78 15-19 0.53 0.32 0.89 20-24 0.83 0.52 1.35 25-29 1.22 0.82 1.81 30-34 0.88 0.61 1.28 35-39 1.13 0.86 1.48 40-44 0.87 0.68 1.13 45-49 0.95 0.77 1.18 50-54 1.02 0.86 1.22 55-59 0.88 0.77 1.00 60-64 0.86 0.77 0.95 65-69 0.74 0.68 0.80 70-74 0.72 0.66 0.77 75-79 0.79 0.73 0.85 80-84 0.71 0.66 0.76 85-89 0.78 0.72 0.84 90-94 0.94 0.86 1.03 95+ 0.91 0.80 1.02
Kremer's "proportional rate ratio" values were expressed as the ratio for Pfizer divided by the ratio for Moderna, so it's the opposite of the Moderna-Pfizer ratios that have been calculated by me and Kirsch. His PRR value for deaths divided by reports was about 1.014 for ages 18-64 and 1.278 for ages 65+:
At first I got about 27% more deaths in ages 65+ than Kremer but only about 4% more reports. But then I realized it was probably because I counted reports where the DIED column had the value Y as deaths, but Kremer might have only counted reports where one of the symptom fields included death, which also gave me about 18% fewer deaths in ages 65+ for Moderna and Pfizer combined:
> va=do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSDATA.csv"))))
> va=merge(va,do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSVAX.csv")))))
> symp=do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSSYMPTOMS.csv"))))
> va$dead=va$VAERS_ID%in%symp$VAERS_ID[rowSums(symp[,c(1:5*2)]=="Death")==1]
> vae=va[AGE_YRS>=65&VAX_TYPE=="COVID19"&grepl("MOD|PFI",VAX_MANU)]
> vae[DIED=="Y",.N,VAX_MANU] # DIED column is Y
VAX_MANU N
1: MODERNA 7598
2: PFIZER\\BIONTECH 7520
> vae[dead==T,.N,VAX_MANU] # symptom columns include death
VAX_MANU N
1: MODERNA 6285
2: PFIZER\\BIONTECH 6138
I don't know which approach is better, but here when I counted reports where the symptom columns included death, I got similar PRR values as Kremer:
> young=va[AGE_YRS%in%18:64&VAX_TYPE=="COVID19",.(dead=sum(dead),reports=.N),.(type=VAX_MANU)]
> young=young[match(c("PFIZER\\BIONTECH","MODERNA"),type)]
> young
type dead reports
1: PFIZER\\BIONTECH 1503 283324
2: MODERNA 1458 273126
> young[1,dead/reports]/young[2,dead/reports]
[1] 0.99376 # PRR for ages 18 to 64 (Pfizer divided by Moderna)
> old=va[AGE_YRS>=65&VAX_TYPE=="COVID19",.(dead=sum(dead),reports=.N),.(type=VAX_MANU)]
> old=old[match(c("PFIZER\\BIONTECH","MODERNA"),type)]
> old
type dead reports
1: PFIZER\\BIONTECH 6138 105769
2: MODERNA 6285 142616
> old[1,dead/reports]/old[2,dead/reports]
[1] 1.3168 # PRR for ages 65 and above (Pfizer divided by Moderna)
VAERS also contains reports from European countries, but several fields including the country codes were scrubbed from European reports in November 2022. WelcomeTheEagle posted a version of the VAERS data here that was downloaded before the country codes were removed: https://www.vaersaware.com/post/vaers-nov-11th-downloadable-files. I used the data to calculate an age-normalized rate of adverse event reports in the Czech Republic. I took the number of doses administered in each age group from the record-level data, so that I only included doses administered up to November 11th 2022 which was the date of the old snapshot of VAERS.
In the old VAERS data you can select Czech reports by selecting rows where the first two characters of the SPLTTYPE field are CZ. At first I tried selecting rows where the SYMPTOM_TEXT field matched identifiers like CZ-CZSUKL-21010694, but it missed about 16% of all reports where the SPLTTYPE field started with CZ. There were only two reports where SYMPTOM_TEXT matched \bCZ- but SPLTTYPE didn't start with CZ, but one of them was a Belgian report that was linked to an earlier Czech report of the same patient, and the other report was missing the SPLTTYPE field.
There were only 5,247 rows where the SPLTTYPE column started with CZ and the VAX_TYPE column was COVID19, and about two thirds of the rows were missing the age field. But among the remaining 1,871 rows, the age-standardized rate of adverse event reports per doses administered was about 3.0 times higher for Moderna than Pfizer, and the crude rate that was not adjusted for age was about 2.9 times higher for Moderna than Pfizer:
age=\(x,y){class(x)=class(y)=NULL;(y-x-(y-789)%/%1461+(x-789)%/%1461)%/%365}
ua=\(x,y,...){u=unique(x);y(u,...)[match(x,u)]}
agecut=\(x,y)cut(x,c(y,Inf),paste0(y,c(paste0("-",y[-1]-1),"+")),T,F)
ages=c(0,30,40,50,seq(60,95,5))
rename=\(x){u=unique(x);o=rep("Other",length(u))
o[grep("comirnaty|pfizer",u,T)]="Pfizer"
o[grep("moderna|spikevax",u,T)]="Moderna"
o[grep("nuvaxovid|novavax",u,T)]="Novavax"
o[grep("janssen",u,T)]="Janssen"
o[grep("vaxzevria|astrazeneca",u,T)]="AstraZeneca"
o[match(x,u)]}
va=fread("AllVAERSDataCSVS/NonDomesticVAERSDATA.csv")
va=merge(va,fread("AllVAERSDataCSVS/NonDomesticVAERSVAX.csv"))
vae=va[grepl("^CZ",SPLTTYPE)&VAX_TYPE=="COVID19"]
rec=fread("Czech/data/CR_records.csv")
k=grep("Datum_",names(rec));k2=k+3
set.seed(0);d=rec[,.(birth=rep(ua(paste0(Rok_narozeni,"-1-1"),as.Date)+sample(0:364,.N,T),7))]
d$date=`class<-`(unlist(rec[,..k],,F),"Date")
d$type=rename(unlist(rec[,..k2],,F))
d=d[!is.na(date)&date<="2022-11-11"]
d[,age:=age(birth,date)]
d=d[,.(doses=.N),.(type,age=agecut(pmax(age,0),ages))]
std=fread("http://sars2.net/f/czcensus2021pop.csv")[,.(std=sum(pop)),.(age=agecut(age,ages))][,std:=std/sum(std)]
a=merge(d,vae[!is.na(AGE_YRS),.(reports=.N),.(type=rename(VAX_MANU),age=agecut(AGE_YRS,ages))])
o=merge(std,a)[,.(reports=sum(reports),doses=sum(doses),agestd=sum(reports/doses*std*1e6)),type]
o[,crude:=reports/doses*1e6][order(-doses)]|>print(r=F)
type reports doses agestd crude Pfizer 1369 15508487 86.98 88.27 # Pfizer has about 3.3 times as many reports as Moderna Moderna 415 1639762 258.43 253.09 Janssen 84 397696 205.40 211.22 Other 2 180 2590.57 11111.11
At first when I tried to select Czech VAERS reports by picking rows where the age field was not missing and the SYMPTOM_TEXT field matched \bCZ-, I got 199 reports for Moderna but 1210 for Pfizer. But now when I picked rows where the SPLTTYPE field started with CZ instead, I got 415 reports for Moderna but 1369 for Pfizer, so it approximately doubled the ratio of Moderna reports to Pfizer reports. It might be because many of the reports for Pfizer vaccines included an identifier like CZ-PFIZER INC-2021285400 in the SYMPTOM_TEXT field, but there wasn't a similar identifier for Moderna vaccines. However it also seems to be a sign that there was somehow a different reporting procedure used for Pfizer vaccines and Moderna vaccines.
In the previous block of code I only included VAERS reports where the age field was not missing. But now when I also included the reports where the age field was missing, Pfizer got about 9% higher crude rate of reports per dose than Moderna:
> allage=merge(d[,.(doses=sum(doses)),type],va[,.(reports=.N),.(type=rename(VAX_MANU))])
> allage[,crude:=reports/doses*1e6][order(-doses)]|>print(r=F)
type doses reports crude
Pfizer 15526217 4548 292.9 # Pfizer has about 10.3 times as many reports as Moderna
Moderna 1639762 441 268.9
Janssen 414410 244 588.8
Other 764 14 18324.6
In the output above Pfizer got about 10.3 times as many reports as Moderna but in the block of output before it the ratio was about 3.3. In the spreadsheet from the GitHub account of Palacký University I got an intermediate ratio of about 6.8. The Palacký spreadsheet includes about 2.7 times as many total reports as VAERS, so there might be some kind of a bias for which reports end up being submitted to VAERS:
> s2=data.table(readxl::read_excel("SUKL-batches-AE.xlsx",sheet=2))
> s2[,.(reports=sum(no_reports)),type=rename(name)][order(-reports)]|>print(r=F)
type reports
Pfizer 10712 # Pfizer has about 6.8 times as many reports as Moderna
Moderna 1575
AstraZeneca 1512
Janssen 555
Other 18
Novavax 14
For some reason in the Czech VAERS data, the percentage of reports that was missing the age field was about 6% for Moderna but 70% for Pfizer (even though in many of the reports where the age field was missing, the age was still mentioned in the SYMPTOM_TEXT field):
> va=fread("AllVAERSDataCSVS/NonDomesticVAERSDATA.csv")
> va=merge(va,fread("AllVAERSDataCSVS/NonDomesticVAERSVAX.csv"))
> vae=va[grepl("^CZ",SPLTTYPE)&VAX_TYPE=="COVID19"]
> vae[,.(reports=.N,missingagepct=mean(is.na(AGE_YRS))*100),VAX_MANU][order(-reports)]|>print(r=F)
VAX_MANU reports missingagepct
PFIZER\\BIONTECH 4548 69.899
MODERNA 441 5.896
JANSSEN 244 65.574
UNKNOWN MANUFACTURER 14 78.571
In the whole VAERS data the percentage of reports that was missing the age field was about 8.9% for Moderna but 9.7% for Pfizer (so it's not a major difference, but it might have still introduced a slight bias to Kremer's analysis where he only included reports within a certain age range):
> va=do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSDATA.csv"))))
> va=merge(va,do.call(rbind,lapply(2020:2024,\(i)fread(paste0(i,"VAERSVAX.csv")))))
> va[VAX_TYPE=="COVID19",.(reports=.N,missingagepct=mean(is.na(AGE_YRS))*100),VAX_MANU][order(-reports)]|>print(r=F)
VAX_MANU reports missingagepct
PFIZER\\BIONTECH 481402 9.738
MODERNA 469655 8.871
JANSSEN 73454 18.101
UNKNOWN MANUFACTURER 7363 27.774
NOVAVAX 492 3.862
I also tried downloading the raw EudraVigilance data from here: https://howbad.info (direct link https://www.dropbox.com/scl/fi/mrzkrrg5391kevsrwex0e/eudra-vigilance3.zip?rlkey=eqg628elz9lqufpehh2ni3o9w). I ran this code to get the number of deaths and reports grouped by age, date, and manufacturer:
awk -F\\t '{sub(" .*","",$3);sub("More than 85 Years",85,$7);sub("Not Specified","",$7);sub("[- ].*","",$7);x=""}$12~/SPIKEVAX|MODERNA/{x="Moderna"}$12~/COMIRNATY|TOZINAMERAN/{x="Pfizer"}$12~/COVID-19 VACCINE ASTRAZENECA|VAXZEVRIA/{x="AstraZeneca"}$12~/COVID-19 VACCINE JANSSEN|JCOVDEN/{x="Janssen"}x{k=$3","x","$7;a[k]++;dead[k]+=$11~/Fatal - Results in Death/}END{for(i in a)print i","a[i]","dead[i]}' eudra-vigilance3/*|(echo date,type,age,reports,dead;sort -nk3|sort -sk2,2|sort -sk1,1)>eudra.csv
However in the EudraVigilance data Moderna got a higher rate of deaths per reports than Pfizer in all age groups:
eu=fread("http://sars2.net/f/eudra.csv")
eu=eu[,.(reports=sum(reports),dead=sum(dead)),.(age=ifelse(age<10,0,age),type)][order(age)]
eu=rbind(eu,eu[,.(reports=sum(reports),dead=sum(dead),age="Total"),type])
o=eu[type=="Moderna",.(age,Moderna_dead=dead,Moderna_reports=reports)]
o=cbind(o,eu[type=="Pfizer",.(Pfizer_dead=dead,Pfizer_reports=reports)])
o$Moderna_rate=o[,2]/o[,3]*1e6
o$Pfizer_rate=o[,4]/o[,5]*1e6
o$Moderna_Pfizer_ratio=o[,6]/o[,7]
o
| Age |
Moderna dead |
Moderna reports |
Pfizer dead |
Pfizer reports |
Moderna rate |
Pfizer rate |
Moderna-Pfizer ratio |
|---|---|---|---|---|---|---|---|
| 0 | 13 | 511 | 74 | 6436 | 25440 | 11498 | 2.21 |
| 12 | 26 | 2553 | 158 | 30958 | 10184 | 5104 | 2.00 |
| 18 | 1832 | 292765 | 3034 | 920149 | 6258 | 3297 | 1.90 |
| 65 | 3139 | 53245 | 5937 | 152645 | 58954 | 38894 | 1.52 |
| 85 | 1367 | 6259 | 4104 | 26633 | 218405 | 154095 | 1.42 |
| NA | 280 | 16845 | 1098 | 74500 | 16622 | 14738 | 1.13 |
| Total | 6657 | 372178 | 14405 | 1211321 | 17887 | 11892 | 1.50 |
The next plot also shows that there's a lot of variance in the "proportional rate ratio" across European countries. Even if you only look at the column which shows the total ratio in all age groups and you look at the countries at the top of the plot which have the highest number of reports, the ratio is far below 1 in Austria and Germany but far above 1 in France, Italy, and Netherlands. The countries are sorted by the total number of reports, which was for some reason the highest in Austria:
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}
agecut=\(x,y)cut(x,c(y,Inf),paste0(y,c(paste0("-",y[-1]-1),"+")),T,F)
iso2=fread("https://raw.githubusercontent.com/lukes/ISO-3166-Countries-with-Regional-Codes/master/all/all.csv")
iso2=setNames(iso2$name,iso2$`alpha-2`)
iso2["NL"]="Netherlands"
va=fread("AllVAERSDataCSVS/NonDomesticVAERSDATA.csv")
va=merge(va,fread("AllVAERSDataCSVS/NonDomesticVAERSVAX.csv"))
symp=fread("AllVAERSDataCSVS/NonDomesticVAERSSYMPTOMS.csv")
va$dead=va$VAERS_ID%in%symp$VAERS_ID[rowSums(symp[,c(1:5*2)]=="Death")==1]
va[,country:=substr(SPLTTYPE,1,2)]
eea=strsplit("AT BE BG HR CY CZ DK EE FI FR DE GR HU IE IT LV LT LU MT NL PL PT RO SK SI ES SE IS LI NO"," ")[[1]]
va=va[VAX_TYPE=="COVID19"&country%in%eea]
ages=c(0,40,60,80)
vae=va[,.(dead=sum(dead),reports=.N),.(type=VAX_MANU,country=iso2[country],age=agecut(AGE_YRS,ages))]
vae[is.na(age),age:="Missing"]
vae=vae[country%in%vae[,sum(reports),country][V1>=5e3,country]]
vae[,country:=factor(country,vae[,sum(reports),country][order(-V1),country])]
vae[,type:=factor(type,vae[,sum(reports),type][order(-V1),type])]
vae=rbind(vae,vae[,.(dead=sum(dead),reports=sum(reports),age="Total"),.(type,country)])
vae=rbind(vae,vae[,.(dead=sum(dead),reports=sum(reports),country="Total"),.(type,age)])
a=vae[,tapply(dead/reports*1e5,.(type,country,age),c)]
type1="MODERNA";type2="PFIZER\\BIONTECH"
m1=a[type1,,];m2=a[type2,,];m=(m1-m2)/ifelse(m1>m2,m2,m1)
disp=m1/m2;disp[]=sprintf("%.*f",ifelse(!is.na(disp)&disp>=10,1,2),disp);disp[is.na(m)]="NA"
maxcolor=10
zero=m1>0&m2==0;m[zero]=-maxcolor;disp[zero]="0.00"
inf=m1==0&m2>0;m[inf]=maxcolor;disp[inf]="Inf"
pal=hsv(rep(c(7/12,0),4:5),c(.9,.75,.6,.3,0,.3,.6,.75,.9),c(.4,.65,1,1,1,1,1,.65,.4))
pheatmap::pheatmap(m,filename="i1.png",display_numbers=disp,
gaps_col=ncol(m)-1,gaps_row=nrow(m)-1,cluster_rows=F,cluster_cols=F,legend=F,
cellwidth=19,cellheight=19,fontsize=9,fontsize_number=8,border_color=NA,na_col="gray90",
number_color=ifelse(!is.na(m)&abs(m)>.5*maxcolor,"white","black"),
breaks=seq(-maxcolor,maxcolor,,256),colorRampPalette(pal)(256))
mpop=vae[type%in%c(type1,type2),xtabs(reports~country+age)]
exp2=.6;maxcolor2=max(mpop[-nrow(mpop),-ncol(mpop)])
pheatmap::pheatmap(mpop^exp2,filename="i2.png",display_numbers=kimi(mpop),
gaps_row=nrow(m)-1,gaps_col=ncol(m)-1,cluster_rows=F,cluster_cols=F,legend=F,
cellwidth=19,cellheight=19,fontsize=9,fontsize_number=8,border_color=NA,na_col="white",
number_color=ifelse(mpop^exp2>maxcolor2^exp2*.45,"white","black"),
breaks=seq(0,maxcolor2^exp2,,256),sapply(seq(1,0,,256),\(i)rgb(i,i,i)))
mdead=vae[type%in%c(type1,type2),xtabs(dead~country+age)]
exp3=.6;maxcolor3=max(mdead[-nrow(mdead),-ncol(mdead)])
pheatmap::pheatmap(mdead^exp3,filename="i3.png",display_numbers=kimi(mdead),
gaps_row=nrow(m)-1,gaps_col=ncol(m)-1,cluster_rows=F,cluster_cols=F,legend=F,
cellwidth=19,cellheight=19,fontsize=9,fontsize_number=8,border_color=NA,na_col="white",
number_color=ifelse(mdead^exp3>maxcolor3^exp3*.45,"white","black"),
breaks=seq(0,maxcolor3^exp3,,256),sapply(seq(1,0,,256),\(i)rgb(i,i,i)))
cap="VAERS reports in EEA up to November 2022: Ratio of deaths per reports for Moderna divided by Pfizer. For example 2.0 means that the ratio of deaths per reports was twice as high for Moderna as Pfizer. Only countries with at least 5,000 total reports were included. An old snapshot of VAERS data where country codes had not been removed was downloaded from here: vaersaware.com/post/vaers-nov-11th-downloadable-files."
system(paste0("pad=20;pad2=30;w=`identify -format %w i1.png`;magick -font Arial-Bold -pointsize 42 -size $[w-pad*2]x \\( caption:'Moderna-Pfizer ratio' -splice $[pad]x20 i1.png -append \\) \\( caption:'Reports for Pfizer and Moderna' -splice $[pad]x20 i2.png -append \\) \\( caption:'Deaths for Pfizer and Moderna' -splice $[pad]x20 i3.png -append \\) -shave 14x +append -trim -bordercolor white -border $pad2 0.png;w=`identify -format %w 0.png`;magick -pointsize 42 -font Arial -interline-spacing -2 \\( -size $[w-pad2*2]x caption:'",cap,"' -splice $[pad2]x24 \\) 0.png -append 1.png"))
The raw data for EudraVigilance that can be downloaded from howbad.info doesn't show the country of the reports. But you can see the number of reports for each country from EMA's portal for EudraVigilance, if you go here, click "C", scroll to COVID vaccines, open the page for a vaccine brand, and click "Number of Individual Cases by EEA countries": https://www.adrreports.eu/en/search_subst.html. I compiled the data for all vaccine brands into this CSV file on August 8th 2024 UTC: f/eudracountry.csv.
I took the number of doses administered by manufacturer and country from OWID: https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/vaccinations/vaccinations-by-manufacturer.csv. OWID was missing data for Greece and UK, and in the following analysis I also excluded Iceland because it was missing most administered doses at OWID.
In most countries Moderna got a higher rate of reports per dose than Pfizer (with the exception of Italy, Netherlands, Portugal, and a few small countries that have a small sample size):
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}
owid=fread("https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/vaccinations/vaccinations-by-manufacturer.csv")
owid=unique(owid[order(-date)],by=c("location","vaccine"))
owid=owid[,.(country=location,doses=total_vaccinations,type=vaccine)]
owid$type[owid$type=="Pfizer/BioNTech"]="Pfizer"
owid$type[owid$type=="Oxford/AstraZeneca"]="AstraZeneca"
owid$type[owid$type=="Johnson&Johnson"]="Janssen"
eu=fread("http://sars2.net/f/eudracountry.csv")
eu$country[eu$country=="Czech Republic"]="Czechia"
me=merge(eu,owid)[country!="Iceland"]
me[,country:=factor(country,me[,.(sum(doses)),country][order(-V1),country])]
me[,type:=factor(type,me[,.(sum(doses)),type][order(-V1),type])]
me=rbind(me,me[,.(reports=sum(reports),doses=sum(doses),country="Total"),type])
me=rbind(me,me[,.(reports=sum(reports),doses=sum(doses),type="Total"),country])
m=me[,tapply(reports/doses*1e4,.(country,type),c)]
m[is.infinite(m)]=NA;disp=kimi(m);exp=.7;maxcolor=200
pal=hsv(c(21,21,21,16,11,6,3,0,0,0)/36,c(0,.25,rep(.5,8)),c(rep(1,8),.5,0))
pheatmap::pheatmap(m^exp,filename="i1.png",display_numbers=disp,
gaps_row=nrow(m)-1,gaps_col=ncol(m)-1,
cluster_rows=F,cluster_cols=F,legend=F,cellwidth=20,cellheight=20,fontsize=9,fontsize_number=8,
border_color=NA,na_col="gray90",
number_color=ifelse(!is.na(m)&m^exp>maxcolor^exp*.8,"white","black"),
breaks=seq(0,maxcolor^exp,,256),
colorRampPalette(pal)(256))
m=me[,xtabs(doses~country+type)]
disp=kimi(m);exp=.7;maxcolor=max(m[-nrow(m),-ncol(m)])
pal=sapply(seq(1,0,,256),\(i)rgb(i,i,i))
pheatmap::pheatmap(m^exp,filename="i2.png",display_numbers=disp,
gaps_row=nrow(m)-1,gaps_col=ncol(m)-1,
cluster_rows=F,cluster_cols=F,legend=F,cellwidth=20,cellheight=20,fontsize=9,fontsize_number=8,
border_color=NA,na_col="gray90",
number_color=ifelse(!is.na(m)&m^exp>maxcolor^exp*.45,"white","black"),
breaks=seq(0,maxcolor^exp,,256),
colorRampPalette(pal)(256))
cap="Number of adverse event reports at EudraVigilance per 10,000 doses administered at OWID. The right side shows the number of doses administered. Source: adrreports.eu/en/search_subst.html (letter C, links for COVID vaccines), ourworldindata.org/grapher/covid-vaccine-doses-by-manufacturer (download CSV data)."
system(paste0("magick \\( i1.png -gravity east -chop 10x \\) \\( i2.png -gravity west -chop 10x \\) +append 0.png;w=`identify -format %w 0.png`;pad=20;magick -pointsize 42 -interline-spacing -2 -font Arial \\( -size $[w-pad*2]x caption:'",cap,"' -splice $[pad]x20 \\) 0.png -append 1.png"))
However I was not able to adjust my analysis for differences in the age distribution between vaccine brands. In the VAERS reports that are not missing the age field, most EEA countries have a higher average age for Moderna than Pfizer, even though in most countries the difference isn't as big as in the Czech Republic:
va=fread("AllVAERSDataCSVS/NonDomesticVAERSDATA.csv")
va=merge(va,fread("AllVAERSDataCSVS/NonDomesticVAERSVAX.csv"))
va[,country:=substr(SPLTTYPE,1,2)]
vae=va[VAX_TYPE=="COVID19",.(age=mean(AGE_YRS,na.rm=T),.N),.(type=VAX_MANU,country)]
o=vae[type=="MODERNA",.(country,Moderna_age=age,Moderna_N=N)]
o=merge(o,vae[type=="PFIZER\\BIONTECH",.(country,Pfizer_age=age,Pfizer_N=N)])
o$diff=o[,2]-o[,4]
o=o[country%in%strsplit("AT BE BG HR CY CZ DK EE FI FR DE GR HU IE IT LV LT LU MT NL PL PT RO SK SI ES SE IS LI NO"," ")[[1]]]
o[order(-Pfizer_N)][,c(1,2,4,6,3,5)]|>mutate_if(is.double,round,1)|>print(r=F)
country Moderna_age Pfizer_age diff Moderna_N Pfizer_N
AT 46.8 46.1 0.7 7179 85854
DE 49.7 49.4 0.3 9772 54896
FR 51.3 50.2 1.2 6781 41596
IT 49.2 48.0 1.2 5161 19949
NL 51.3 47.5 3.8 1743 11278
SE 49.9 49.7 0.3 2190 9454
ES 45.3 48.0 -2.7 2232 8977
BE 42.1 42.6 -0.5 2871 8285
PT 44.6 44.9 -0.3 885 6337
NO 44.5 46.3 -1.8 1563 5934
FI 45.6 47.1 -1.5 1059 5555
GR 49.6 47.8 1.8 676 5125
CZ 54.3 48.7 5.6 441 4548
IE 47.5 42.1 5.4 794 4408
PL 47.5 44.2 3.3 500 4300
DK 47.6 51.5 -3.9 861 4250
RO 46.5 41.2 5.3 180 1978
EE 55.7 49.4 6.3 104 1143
HR 47.7 44.7 3.0 178 1118
HU 61.8 49.2 12.6 93 1064
IS 41.2 47.6 -6.4 76 1048
SK 52.8 46.1 6.8 120 857
LT 53.1 47.8 5.4 63 830
SI 62.8 50.6 12.2 58 518
LU 45.7 43.4 2.3 117 457
LV 53.0 40.8 12.2 149 417
CY 43.2 45.5 -2.3 26 214
BG 46.2 46.5 -0.3 33 185
MT 43.2 50.9 -7.7 23 137
LI 44.0 72.0 -28.0 2 5
country Moderna_age Pfizer_age diff Moderna_N Pfizer_N
In February 2024 the Dutch Central Bureau of Statistics published a dataset for COVID and non-COVID ASMR by vaccination status: https://www.cbs.nl/nl-nl/longread/rapportages/2024/covid-vaccinatiestatus-en-sterfte/3-resultaten. It included one set of data where people were classified as vaccinated immediately after the first dose, and another set of data where people were generally only classified as vaccinated 2 weeks after the second dose (or 4 weeks after the first dose in the case of vaccine brands like Johnson & Johnson which had only one primary course dose).
The classification delay reduced unvaccinated ASMR during the main rollout of the primary course doses from March to July 2021, so it resulted in a big advantage to unvaccinated people. In January and February 2021 the classification delay gave an advantage to vaccinated people because then there was a delay before the vulnerable groups of people were priorized during early rollout were counted as vaccinated, even though there were so few people who had been vaccinated during the first two months of 2021 that they don't have much impact on the total ASMR of vaccinated people:
Similarly in the Czech record-level data when I classified people as unvaccinated until 3 weeks from the first dose, it gave an advantage to unvaccinated people because it reduced the ASMR of unvaccinated people and it increased the ASMR of vaccinated people:
library(data.table);library(ggplot2)
b=fread("http://sars2.net/f/czbucketskeep.csv.gz")[dose<=1]
std=fread("http://sars2.net/f/czcensus2021pop.csv")[,.(std=sum(pop)),.(age=pmin(age,95)%/%5*5)][,std:=std/sum(std)]
p=copy(b)[,dose:=ifelse(dose==1,"Vaccinated","Unvaccinated")]
p=rbind(p,copy(b)[week<3,dose:=0][,dose:=ifelse(dose==1,"Vaccinated at least 21 days ago","Unvaccinated or vaccinated less than 21 days ago")])
p=p[,.(cmr=sum(dead)/sum(alive)*365e5),.(dose,age=pmin(age,95)%/%5*5,month)]
p=merge(std,p)[,.(asmr=sum(cmr*std)),.(dose,month)]
p[,dose:=factor(dose,unique(dose))]
p[,month:=as.Date(paste0(month,"-1"))]
xstart=as.Date("2020-12-1");xend=as.Date("2023-1-1")
p=p[month>=xstart&month<=xend]
xbreak=seq(as.Date("2021-1-1"),xend,"6 month");xlab=c(rbind("",2021:2022),"")
ystep=500;ystart=0;yend=max(p$asmr,na.rm=T);ybreak=seq(ystart,yend,ystep)
color=c(hcl(c(210,0)+15,130,60),hcl(c(210,0)+15,80,40))
label=data.frame(x=xstart+.975*(xend-xstart),y=seq(.92*yend,,yend/-12,nlevels(p$dose)),label=levels(p$dose))
sub="The resident population estimates from the 2021 census by 5-year age groups were used as the standard population."|>stringr::str_wrap(78)
ggplot(p,aes(x=month+15,y=asmr))+
geom_vline(xintercept=seq(as.Date("2021-1-1"),as.Date("2023-1-1"),"year"),linewidth=.3)+
geom_vline(xintercept=c(xstart,xend),linewidth=.3,lineend="square")+
geom_hline(yintercept=c(ystart,yend),linewidth=.3,lineend="square")+
geom_line(aes(color=dose,linetype=dose),linewidth=.4)+
geom_point(aes(color=dose),size=.6)+
geom_label(data=label,aes(x=x,y=y,label=label),fill=alpha("white",.85),label.r=unit(0,"lines"),label.padding=unit(.04,"lines"),label.size=0,size=2.75,hjust=1,color=color)+
labs(x=NULL,y=NULL,title="Czech record-level data: All-cause ASMR per 100,000 person-years",subtitle=sub)+
scale_x_date(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak,labels=\(x)formatC(x,digits=0,format='f',big.mark=','))+
scale_color_manual(values=color)+
scale_linetype_manual(values=c(1,1,2,2))+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=7,color="black"),
axis.text.x=element_text(margin=margin(3)),
axis.ticks=element_line(linewidth=.3),
axis.ticks.length=unit(3,"pt"),
axis.ticks.length.x=unit(0,"pt"),
axis.title=element_text(size=8),
legend.position="none",
panel.background=element_blank(),
panel.grid=element_blank(),
plot.subtitle=element_text(size=7,margin=margin(,,4)),
plot.title=element_text(size=7.5,face="bold",margin=margin(1,,4)),
plot.margin=margin(4,5,4,4))
ggsave("1.png",width=3.9,height=2.4,dpi=350*4)
system("mogrify -resize 25% 1.png")
And in the English ONS data when I classified people as unvaccinated until 3 weeks from the first dose, it similarly gave a big advantage to unvaccinated people, because it reduced unvaccinated CMR in early 2021 when the first dose was rolled out (R code):
In an early version of the README file of Kirsch's GitHub repository about the Czech data, Kirsch wrote: "Vaccines were randomly distributed for those wishing to get vaxxed. [...] People were not allowed to select which vaccines they got. [...] The randomization of which vaccine someone got created a perfect real-world randomized clinical trial where we could compute the mortality rates for 1 year after Dose 2 for the two most popular vaccines." [https://github.com/skirsch/Czech/blob/5725ac1b64ede7124e00b72af68892f31736b349/README.md] After I showed him why his claim was likely wrong, he edited the README file to remove the claim.
However in August 2024 when Kirsch was on InfoWars, he resurrected his old claim that people were not allowed to choose if they got a Moderna or Pfizer vaccine. He said: "The only difference in the Czech Republic between getting a Pfizer shot and a Moderna shot is just pure luck. It was... you weren't allowed to choose, you were just assigned a brand. So it's effectively like the world's largest randomized trial." [https://x.com/RealAlexJones/status/1828912865107231096, time 45:51]
When I called Kirsch out for his comments on InfoWars, he asked me to show evidence that people were allowed to choose which vaccine type they got. But I showed him that there was a dropdown menu for choosing your preferred vaccine type in a form for registering a vaccination appointment: [https://ockodoc.mzcr.cz/wp-content/uploads/2023/07/manual_crs.pdf; translated with Yandex]
However I don't know if people were not initially allowed to choose the vaccine type, because an article dated May 2021 said: "You cannot currently choose which vaccine to get. Odds are you will get the Pfizer/BioNTech vaccine, as that one so far has been the most widely available, accounting for about three-fourths of the deliveries." [https://www.expats.cz/czech-news/article/vaccination-in-the-czech-republic-all-you-need-to-know]
Tomáš Fürst posted this plot which shows that among people who were born in 1929 or earlier, vaccinated people had higher mortality than unvaccinated people starting from late 2021: [https://smis-lab.cz/2024/08/17/celkova-umrtnost-ceskych-obyvatel-dle-vakcinacniho-statutu-nejstarsi-rocniky/]
In the following code I took data from one of my buckets files where I had assigned a random birthday for each person so I could account for the aging of people over time with a sub-year precision. From the last column which shows the total mortality rate in 2021-2022, you can see that unvaccinated people had higher mortality than vaccinated people up to age 94, but in ages 95 and above vaccinated people had higher mortality:
b=fread("http://sars2.net/f/czbucketsdaily.csv.gz")
agecut=\(x,y)cut(x,c(y,Inf),paste0(y,c(ifelse(head(y,-1)==y[-1]-1,"",paste0("-",y[-1]-1)),"+")),T,F)
ages=c(0:9*10,91:105)
a=b[date>="2021-1-1",.(alive=sum(alive),dead=sum(dead)),.(dose=pmin(dose,1),age=agecut(age,ages),month=substr(date,1,7))]
a=rbind(a,a[,.(dead=sum(dead),alive=sum(alive),month="Total"),.(dose,age)])
ar=a[,tapply(dead/alive,.(dose,age,month),c)]
m1=ar[1,,];m2=ar[2,,]
m=(m2-m1)/ifelse(m2>m1,m1,m2)*100
disp=round(m2/m1*100)
maxcolor=600;exp=.85
pheatmap::pheatmap(abs(m)^exp*sign(m),filename="0.png",display_numbers=disp,
gaps_col=ncol(m)-1,gaps_row=9,
cluster_rows=F,cluster_cols=F,legend=F,cellwidth=19,cellheight=19,fontsize=9,fontsize_number=8,
border_color=NA,na_col="gray90",
number_color=ifelse((abs(m)^exp>.55*maxcolor^exp)&!is.na(m),"white","black"),
breaks=seq(-maxcolor^exp,maxcolor^exp,,256),
colorRampPalette(hsv(rep(c(7/12,0),5:4),c(.9,.75,.6,.3,0,.3,.6,.75,.9),c(.4,.65,1,1,1,1,1,.65,.4)))(256))
cap="Czech record-level data: Vaccinated CMR as percentage of unvaccinated CMR"
system(paste0("w=`identify -format %w 0.png`;pad=24;convert -pointsize 42 -interline-spacing -3 -font Arial-Bold \\( -size $[w-pad*2]x caption:'",cap,"' -splice $[pad]x20 \\) 0.png -append 1.png"))
Elderly people who live in a senior home have higher ASMR than elderly people who live on their own, and I don't know if for example elderly people who live in a senior home are more likely to be vaccinated than age-matched elderly people who live on their own.
The plot above shows that for some reason in ages 104 and 105+ vaccinated people had more than 10 times as high mortality as unvaccinated people. For people who were born in 1917 so they turned 104 in 2022, people who were born near the end of the year were much more likely to be alive in 2022 than people who were born near the beginning of the year, but my code for assigning a random birthday distributed the birthdays evenly across the year which resulted in a fairly large bias in the case of the very oldest age groups. However it shouldn't matter that much here because I'm simply comparing vaccinated mortality rate against unvaccinated mortality rate. So I don't know how to explain the extremely high percentages in ages 104 and 105+.
But in any case Fürst could've mentioned that there was a sharp increase in the ratio of unvaccinated to vaccinated CMR around ages 94 to 95, but in younger ages unvaccinated people had higher mortality.
If you look at the raw number of deaths by weeks after vaccination but you don't normalize it for the background mortality rate, dose 1 seems to have a bump in deaths around weeks 2 to 10 that's missing for doses 2 and 3. But that's because during the first weeks after dose 1 the background mortality rate was falling rapidly because the COVID wave in early 2021 was passing by. So here when I calculated a baseline in a way where I accounted for the background mortality rate, the baseline for dose 1 has a much bigger drop during the first 15 weeks after vaccination than the baselines of doses 2 and 3:
In the United States first doses were also given at a time when the COVID wave in early 2021 was passing by, which I think explains why the temporal HVE appears to last only about 3 weeks in Kirsch's Medicare data for 2021. But the temporal HVE seems to last longer in the Medicare data for 2022 which is missing the bump around the third week: [moar.html#US_Medicare_data]
library(data.table);library(ggplot2);library(stringr)
b=fread("http://sars2.net/f/czbucketskeep.csv.gz")[,age:=pmin(age,95)%/%5*5]
d=merge(b[dose%in%1:4],b[dose<=1,.(base=sum(dead)/sum(alive)),.(month,age)])[,base:=base*alive]
p=d[,.(y=sum(dead),base=sum(base),pop=sum(alive)),.(x=week,z=paste0("Dose ",dose))]
xstart=0;xend=105;xstep=5;ystart=0;yend=3000;ystep=500
color=hsv(c(12,21,24,30)/36,c(.4,.6,.8,.8),c(.9,1,.7,.6))
sub=paste0(" ","The dashed line is a baseline which was calculated by multiplying the number of person-days for each combination of 5-year age group and ongoing month with the mortality rate among all people who are included in the record-level data for the same combination of age group and ongoing month."|>str_wrap(90))
sub=paste0(sub,"\n ","People are kept under earlier doses after a new dose. Week 0 consists of the day of vaccination and the next 6 days."|>str_wrap(90))
ggplot(p,aes(x,y))+
geom_hline(yintercept=c(ystart,0,yend),linewidth=.3,lineend="square")+
geom_vline(xintercept=c(xstart,xend),linewidth=.3,lineend="square")+
geom_line(aes(color=z),linewidth=.4)+
geom_line(aes(y=base,color=z),linetype=2,linewidth=.4)+
labs(title="Czech record-level data: Deaths by weeks after vaccination",subtitle=sub,x=NULL,y=NULL)+
scale_x_continuous(limits=c(xstart,xend),breaks=seq(xstart,xend,xstep))+
scale_y_continuous(limits=c(ystart,yend),breaks=seq(ystart,yend,ystep))+
scale_color_manual(values=color)+
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"),
legend.background=element_rect(fill="white",color="black",linewidth=.3),
legend.justification=c(1,1),
legend.key.height=unit(11,"pt"),
legend.key.width=unit(15,"pt"),
legend.key=element_rect(fill="white"),
legend.margin=margin(-2,7,3,6),
legend.position=c(1,1),
legend.spacing.x=unit(2,"pt"),
legend.text=element_text(size=7),
legend.title=element_blank(),
panel.grid.major=element_line(linewidth=.3,color="gray90"),
panel.background=element_rect(fill="white"),
plot.margin=margin(4,10,4,4),
plot.subtitle=element_text(size=7,margin=margin(0,0,4,0)),
plot.title=element_text(size=7.5,face="bold",margin=margin(2,0,4,0)))
ggsave("0.png",width=4.5,height=3.2,dpi=400*4)
system("magick 0.png -filter lanczos2 -resize 25% 1.png")
Kirsch picked 9 batches that were given to at last 100 people in ages 70 to 74 in February 2021, and he then calculated a mortality rate for each batch during the first year from vaccination: [https://kirschsubstack.com/p/covid-vax-some-batches-increase-your]
Kirsch wrote that "there was >9X variation between the most deadly Moderna vaccine batch and the safest Pfizer vaccine batch". However actually the ratio of the two batches that are shown in red in his table is only about 8.5 (from .10989/.013001). The ratio becomes over 9 if you compare the Moderna vaccines with a missing batch ID against the Pfizer batch with the highest rate (.12/.013001). But I don't think it makes sense to treat the doses with a missing batch ID as a distinct batch.
But anyway, the sample sizes in his table are way too small because he restricted his analysis to only a single age group and a single month of vaccination.
In the code below I calculated a normalized excess mortality percent for each batch up to the end of 2022, so that I multiplied the number of person-days for each combination of 5-year age group and ongoing month with the mortality rate among all people who are included in the record-level data for the same combination of 5-year age group and ongoing month. I kept people included under a batch even after a dose from another batch.
I got about 19% excess mortality for the Moderna batch 300042721 which was Kirsch's batch with the highest mortality rate, and I got about -16% excess mortality for the Pfizer batch EJ6797 which was Kirsch's batch with the lowest mortality rate. So the Pfizer batch had only about 1.39 times higher normalized mortality (from 119.64/86.28):
> library(data.table)
> download.file("http://sars2.net/f/czbucketsbatchkeep.csv.xz","czbucketsbatchkeep.csv.xz")
> b=fread("xz -dc czbucketsbatchkeep.csv.xz")[,age:=pmin(age,95)%/%5*5]
> a=merge(b[dose>0&batch!=""&type!="Other"],b[dose<=1,.(base=sum(dead)/sum(alive)),.(age,month)])[,base:=base*alive]
> a=a[,.(alive=sum(alive),dead=sum(dead),base=sum(base)),.(batch,age,type)]
> p=a[,.(excess=(sum(dead)/sum(base)-1)*100),.(batch,type)]
> print(p[batch%in%c("300042721","EJ6797")][order(-excess)],nrow=Inf)
batch type excess
1: 300042721 Moderna 19.63627
2: EJ6797 Pfizer -16.28020
I posted code for generating the czbucketsbatchkeep.csv file here: czech2.html#Statistics_by_batch. In the code I assigned a random birthday to each person so that I could model the aging of people with a daily precision.
Kirsch wrote that for people in ages 60 to 64 who got a second dose in June 2021, the Pfizer batch FC2473 got more than 3 times higher mortality rate within a year from vaccination than the Pfizer batch FD0168:
However using the code above my total normalized excess mortality was about -32% for the FC batch and about -39% for the FD batch, so there was not much difference. I included all doses instead of only second doses like Kirsch. I included all age groups and all vaccination months, and I calculated the excess mortality up to the end of 2022 and not during the first year from vaccination.
Kirsch got 48 deaths out of 6238 doses for the FC batch. The upper bound of the 95% confidence interval for the ratio is about twice as the high as the lower bound:
> as.numeric(prop.test(48,6238)$conf) [1] 0.005739933 0.010277757
In the following block of code I picked people who got a second dose from the FC or FD batch in June 2021 and I divided the number of deaths in the 365 days starting from the day of vaccination with the number of doses from the batch. I treated the age of each person as their year of birth subtracted from 2021. I kept people included under a batch even after they got a dose from another batch.
I converted the batch IDs to uppercase and stripped trailing spaces, which is what Kirsch also did. [https://github.com/skirsch/Czech/blob/main/code/vax.py#L59]
When I calculated the ratio of deaths per doses for each 5-year age group, the FD batch actually got a higher mortality rate than the FC batch in each age group from 65-69 to 85-89 (which all had a fairly large sample size of deaths):
rec=fread("CR_records.csv")
d=rec[,.(date=as.IDate(unlist(rec[,.SD,.SDcols=patterns("Datum_")],,F)),dead=DatumUmrti,age=2021-Rok_narozeni)]
d[,batch:=unlist(rec[,.SD,.SDcols=patterns("Sarze_")],,F)][,batch:=toupper(trimws(batch)),batch]
d$dose=rep(1:7,each=nrow(rec))
a=d[dose==2&date>="2021-06-01"&date<"2021-07-01"]
a=a[,.(dead=sum(!is.na(dead)&dead-date<365),doses=.N),.(batch,age=pmin(age,95)%/%5*5)]
a1=a[batch=="FC2473",.(age,dead1=dead,doses1=doses,rate1=dead/doses*1e5)]
a2=a[batch=="FD0168",.(age,dead2=dead,doses2=doses,rate2=dead/doses*1e5)]
o=merge(a1,a2)
o[,rate1vs2:=rate1/rate2]
o[,low:=exp(log(rate1vs2)-qnorm(.975)*sqrt(1/dead1+1/dead2))]
o[,high:=exp(log(rate1vs2)+qnorm(.975)*sqrt(1/dead1+1/dead2))]
print(round(o,2),row.names=F)
# FC2473 # FD0168 # FC vs FD with 95% CI age dead1 doses1 rate1 dead2 doses2 rate2 rate1vs2 low high 15 1 179 558.66 1 709 141.04 3.96 0.25 63.32 20 1 970 103.09 2 2280 87.72 1.18 0.11 12.96 25 0 1706 0.00 1 3469 28.83 0.00 0.00 NaN 30 2 2069 96.67 5 4877 102.52 0.94 0.18 4.86 35 1 2682 37.29 1 6452 15.50 2.41 0.15 38.46 40 8 4111 194.60 11 11006 99.95 1.95 0.78 4.84 45 16 4993 320.45 19 19186 99.03 3.24 1.66 6.29 50 29 4482 647.03 66 28524 231.38 2.80 1.81 4.33 55 38 5404 703.18 148 65597 225.62 3.12 2.18 4.45 60 84 8560 981.31 318 75389 421.81 2.33 1.83 2.96 65 240 31491 762.12 474 58382 811.89 0.94 0.80 1.10 # ratio below 1 70 203 16682 1216.88 301 17613 1708.96 0.71 0.60 0.85 # ratio below 1 75 109 3595 3031.99 213 6880 3095.93 0.98 0.78 1.23 # ratio below 1 80 54 934 5781.58 156 2365 6596.19 0.88 0.64 1.19 # ratio below 1 85 44 440 10000.00 115 1017 11307.77 0.88 0.62 1.25 # ratio below 1 90 36 160 22500.00 59 335 17611.94 1.28 0.84 1.93 95 18 46 39130.43 14 50 28000.00 1.40 0.70 2.81
My ratio in ages 60 to 64 was only about 2.33 and not 3.15 like in Kirsch's table. But I eventually figured out it was because of these three reasons:
df['date_1_'] regardless of the dose number: df['days_until_death'] = (df[dod_col] - df['date_1_']).dt.days.round(0). [https://github.com/skirsch/Czech/blob/eb8406c7c160f3951d798f5cd72a7c84901a60c4/code/vax.py#L117] Kirsch fixed the bug after I pointed it out to him.Kirsch later added this section to his Substack post:
The plot above is a histogram where the y-axis values are always integers, so it's confusing how the y-axis labels include values like 0.5 and 1.5.
Kirsch wrote that "if the vaccines were safe, the mortality rates of the batches I looked at (given to a given group in a specific time period) must have a normal distribution". But that's not necessarily the case if there's confounding factors that vary across batches even after age and month of vaccination are controlled for.
Kirsch's plot had so few batches that it was difficult to even tell if his plot followed a normal distribution or not. And his sample size of total deaths was also fairly small because he only looked at vaccines given in March 2021, but he could've increased his sample size if he would've included all months of vaccination in his plot.
In the histogram plot below I calculated a normalized mortality rate for each batch the same way as earlier, except I plotted the mortality rate as a percentage of the baseline and not as excess mortality relative to the baseline (so now the values are the same as excess mortality values plus 100). I included only batches with at least 10,000 total person-years. But now in my histogram the distribution of mortality rates looks closer to a normal distribution than in Kirsch's plot (even though in plots like this which show excess mortality as a percentage of a baseline, the distribution is typically skewed because the values below 100% typically have a more narrow spread than the values above 100%):
b=fread("xz -dc czbucketsbatchkeep.csv.xz")[,age:=pmin(age,95)%/%5*5]
a=merge(b[dose>0&batch!=""&type!="Other"],b[dose<=1,.(base=sum(dead)/sum(alive)),.(age,month)])[,base:=base*alive]
a=a[,.(alive=sum(alive),dead=sum(dead),base=sum(base)),.(batch,age,type)]
keep=a[,sum(alive),batch][V1/365>=1e4,batch]
a=a[batch%in%keep,.(mort=sum(dead)/sum(base)*100),batch]
step=10;lev=seq(0,max(a$mort),step);yend=30
p=data.table(x=paste0("[",lev,",",lev+step,")"),y=as.integer(table(factor(a$mort%/%step*step,lev))))
p[,x:=factor(x,unique(x))]
png("1.png",1300,800,res=200)
par(mar=c(5,2.3,1.9,0),mgp=c(0,0.6,0))
tit="Histogram of batch mortality percentage relative to baseline"
barx=barplot(p$y,names.arg=p$x,main=tit,cex.main=1.1,xlab=NA,ylab=NA,lwd=1.5,las=2,space=0,ylim=c(0,yend))
text(x=barx,y=p$y-.01*yend,labels=p$y,pos=3,cex=0.8)
dev.off()
In the next plot I calculated ASMR using 5-year age groups but I didn't normalize it in any way for ongoing month or month of vaccination. I again included only batches with at least 10,000 person-years. However now the distribution looked less similar to a normal distribution:
b=fread("xz -dc czbucketsbatchkeep.csv.xz")
b[,age:=ifelse(age<25,0,pmin(age,95)%/%5*5)]
std=fread("czcensus2021pop.csv")[,.(std=sum(pop)),.(age=ifelse(age<25,0,pmin(age,95)%/%5*5))][,std:=std/sum(std)]
a=b[dose>0&batch!=""&type!="Other",.(dead=sum(dead),alive=sum(alive)),.(batch,age)]
a=merge(std,a)[alive/365>=1e4,.(asmr=sum(dead/alive*std*365e5)),batch]
step=50;max=max(a$asmr);lev=seq(0,max,step);yend=10
p=data.table(x=paste0("[",lev,",",lev+step,")"),y=as.integer(table(factor(a$asmr%/%step*step,lev))))
p[,x:=factor(x,unique(x))]
png("1.png",1300,800,res=200)
par(mar=c(5,2.3,1.9,0),mgp=c(0,.6,0))
tit="Histogram of batch ASMR"
barx=barplot(p$y,names.arg=p$x,main=tit,cex.main=1.1,xlab=NA,ylab=NA,lwd=1.5,las=2,space=0,ylim=c(0,yend))
text(x=barx,y=p$y-.01*yend,labels=p$y,pos=3,cex=.8)
dev.off()
Kirsch said that in my triangle plot which showed an age-normalized Moderna-Pfizer ratio by month of vaccination and ongoing month, the results were not statistically significant from July 2021 onwards: [https://kirschsubstack.com/p/a-closer-look-at-the-czech-data-confirms]
Kirsch seems to have determined statistical significance on each row of my heatmap by calculating the mean and standard deviation of each column and then looking at whether the mean was at least 2 standard deviations above 1, so for example in July 2021 he calculated the mean and standard deviation for a vector of monthly values like 0.64, 1.07, 1.72, 1.52, 0.88, and so on:
I don't think his approach was correct, because he should've taken into account the number of deaths and baseline deaths for Moderna and Pfizer vaccinees. I'm not sure if the following is the correct approach either, but when I asked ChatGPT how to calculate the confidence interval for the Moderna-Pfizer ratio for people vaccinated in July 2021, it gave me code like this, where the 95% CI was about 1.20 to 1.53 so the lower bound of the CI was still clearly above 1:
dead1=284;dead2=2266 base1=252.52;base2=2725.01 rate1=dead1/base1;rate2=dead2/base2 ratio=rate1/rate2 se1=sqrt(dead1)/base1;se2=sqrt(dead2)/base2 se=sqrt((se1/rate1)^2+(se2/rate2)^2) exp(log(ratio)+c(-1,1)*qnorm(.975)*se) # [1] 1.195500 1.530073
When I asked ChatGPT how to calculate whether the result was statistically significant, it said "To assess if the result is statistically significant, we check if the 95% confidence interval for the ratio excludes 1. If the interval excludes 1, it indicates that there is a statistically significant difference between the Moderna and Pfizer rates."
When I asked how to calculate a p-value for the result, ChatGPT answered: "To calculate the p-value, we need to test the null hypothesis that the ratio of the two rates is equal to 1. This can be done using a Wald test, which compares the observed value of the ratio to the null value (1), using the standard error of the ratio to compute a Z-score. The p-value is then obtained from the Z-score." And it gave me this code which gave me a p-value of about 2e-6:
z_score=log(ratio)/se 2*(1-pnorm(abs(z_score))) # [1] 1.6133e-06
Here when I applied the same method of calculating confidence intervals to all months of vaccination, I got fairly narrow confidence intervals up to January 2022:
b=fread("http://sars2.net/f/czbucketskeep.csv.gz")
b=b[month>="2020-12"&dose<=1][,age:=pmin(age,100)]
t=b[,.(dead=sum(dead),alive=as.numeric(sum(alive))),.(age,month,vaxmonth,type)]
t=merge(t,t[,.(base=sum(dead)/sum(alive)),age])[,base:=base*alive]
a=t[,.(dead=sum(dead),base=sum(base)),.(vaxmonth,type)]
a1=a[type=="Moderna",.(vaxmonth,dead1=dead,base1=base)]
a2=a[type=="Pfizer",.(vaxmonth,dead2=dead,base2=base)]
o=merge(a1,a2)
rate1=o[,dead1/base1];rate2=o[,dead2/base2]
o[,ratio:=rate1/rate2]
se1=o[,sqrt(dead1)/base1];se2=o[,sqrt(dead2)/base2]
se=sqrt((se1/rate1)^2+(se2/rate2)^2)
o[,low:=exp(log(ratio)-qnorm(.975)*se)]
o[,high:=exp(log(ratio)+qnorm(.975)*se)]
mutate_if(o,is.double,round,2)|>print(r=F)
vaxmonth dead1 base1 dead2 base2 ratio low high 2021-01 1515 744.34 18513 21218.40 2.33 2.21 2.46 2021-02 3501 3732.26 13979 19466.04 1.31 1.26 1.36 2021-03 4961 6404.36 19616 32300.32 1.28 1.24 1.32 2021-04 2887 3748.78 16321 26352.89 1.24 1.20 1.29 2021-05 2186 2433.25 12130 21496.87 1.59 1.52 1.67 2021-06 1674 1510.94 3761 5030.32 1.48 1.40 1.57 2021-07 284 252.52 2266 2725.01 1.35 1.20 1.53 2021-08 136 125.80 1080 1118.41 1.12 0.94 1.34 2021-09 157 93.80 744 586.62 1.32 1.11 1.57 2021-10 99 78.51 754 557.74 0.93 0.76 1.15 2021-11 282 222.76 2275 1838.82 1.02 0.90 1.16 2021-12 249 163.23 1306 871.66 1.02 0.89 1.17 2022-01 81 47.04 585 342.17 1.01 0.80 1.27 2022-02 13 8.14 125 55.09 0.70 0.40 1.25 2022-03 2 3.08 46 17.62 0.25 0.06 1.03 2022-04 0 0.81 19 10.48 0.00 NaN NaN 2022-05 0 0.45 20 8.44 0.00 NaN NaN 2022-06 0 0.25 12 4.52 0.00 NaN NaN 2022-07 2 0.77 10 4.49 1.17 0.26 5.34 2022-08 1 0.87 21 8.76 0.48 0.06 3.55 2022-09 1 0.53 20 7.98 0.76 0.10 5.66 2022-10 0 0.12 14 3.93 0.00 NaN NaN 2022-11 1 0.07 3 1.45 6.85 0.71 65.81 2022-12 0 0.00 0 0.26 NaN NaN NaN vaxmonth dead1 base1 dead2 base2 ratio low high
The record-level dataset contains an entry for a man whose year of birth is 1912 but who got vaccinated in 2023 when he would've been 110 or 111 years old:
> t=fread("curl https://github.com/skirsch/Czech/blob/main/data/CR_records.csv.xz|xz -dc")
> t[Rok_narozeni<1913,1:4]
Pohlavikod Rok_narozeni DatumUmrti Datum_1
1: M 1912 <NA> 2021-09-24
2: F 1912 2021-07-22 2021-02-11
3: F 1911 <NA> 2021-07-27
4: F 1911 2020-12-06 <NA>
5: F 1910 <NA> 2021-03-28
6: F 1912 2020-05-12 <NA>
7: F 1912 2021-10-22 <NA>
8: M 1912 <NA> 2023-11-29
9: F 1912 <NA> <NA>
10: F 1912 <NA> <NA>
11: F 1912 <NA> <NA>
12: F 1912 <NA> <NA>
However according to the Czech-language Wikipedia, there have not been male Czech citizens older than 107 years since at least 2010. [https://cs.wikipedia.org/wiki/Seznam_nejstar%C5%A1%C3%ADch_obyvatel_%C4%8Ceska]
In November 2024 ÚZIS (the Czech institute of health information and statistics) published a new dataset which contains columns for infections and hospitalizations. However the year of birth is only given in 5-year ranges and dates of vaccination and death are given on a weekly and not daily precision. SMIS published the following blog post about the data (translated from Czech): [https://smis-lab.cz/2024/11/07/dalsi-velka-datova-sada-uzis-zverejnena/]
After a long time, we received another large data set from ÚZIS. At first glance, this seems like a holy grail, as the field descriptions indicate that the data contains individual-level information on birth quintile, gender, calendar weeks of all covid vaccines, calendar week of death, and cause of death (covid/non-covid). In addition, the data includes information on tests and covid hospitalizations. No country in the world has yet (completely incomprehensible to us) published such a data set.
However, a closer look at the data raises many doubts.
1. The file contains 12,597,668 lines , which is suspiciously large. Even if the second, third, and other infections that are marked as duplicates in the data description are excluded, that leaves 12,125,969 rows, which is still a lot. There are not that many people in the Czech Republic, even if all the dead since 2020 are included.
2. The file does not contain an (anonymized) person identifier , so it is not possible to match lines that belong to the same person. It is therefore not possible to link the data on vaccination with the data on death and its cause. Therefore, it is also impossible to find out what duplications are the reason for the previous point.
3. When examining the "DatumUmrtiLPZ" column, we find that the numbers listed there approximately correspond to the course of death according to the CZSO. Until the end of 2022, the difference in the weekly number of deaths compared to the CZSO will not exceed one percent, but in the second half of 2024 it will already grow to tens of percent. Data on deaths are therefore usable until the end of 2023 at most. However, even the differences compared to the CZSO in the statistics for the weeks of 2020-2023 are suspicious and should not be there.
4. A far bigger problem, however, is that a total of 1,556,198 lines of the file contain no data (except for DCCI=0), i.e. neither birth nor gender data. Similarly, we have 63,606 rows that contain only the date of death, again without data on birth and sex. A substantial part of these deaths occurred in 2020, and when compared with the data from CSU, it is clear that these deaths actually occurred. Not only can these deaths not be matched to a cause of death, but they cannot even be assigned to a given age group and gender. It is therefore not possible to calculate any mortality rate, as it is always dramatically dependent on age and sex.
Thus, this data set is worthless for analysis. At the same time, its investigation already cost several statisticians hours of time. We will inform the ÚZIS about the problems described above through the applicant and request the publication of data in the same format as here, with the fact that the relevant period will be extended at least until the end of 2023 and a column will be added as to whether the death was caused by covid or not.
The file contains these columns:
| column | translation |
|---|---|
| ID | ID |
| Infekce | order of infection |
| Pohlavi | sex |
| RokNarozeni | year of birth |
| DatumPozitivity | date of positivity |
| DatumVysledku | date of result |
| Vylecen | week of recovery |
| Umrti | week of death |
| Symptom | was symptomatic at time of test |
| TypTestu | type of test |
| Datum_Prvni_davka | first dose date |
| Datum_Druha_davka | second dose date |
| Datum_Treti_davka | third dose date |
| Datum_Ctvrta_davka | fourth dose date |
| Datum_Pata_davka | fifth dose date |
| Datum_Sesta_davka | sixth dose date |
| Datum_Sedma_davka | seventh dose date |
| OckovaciLatkaKod_Prvni_davka | first dose vaccine type code |
| OckovaciLatkaKod_Druha_davka | second dose vaccine type code |
| OckovaciLatkaKod_Treti_davka | third dose vaccine type code |
| OckovaciLatkaKod_Ctvrta_davka | fourth dose vaccine type code |
| OckovaciLatkaKod_Pata_davka | fifth dose vaccine type code |
| OckovaciLatkaKod_Sesta_davka | sixth dose vaccine type code |
| OckovaciLatkaKod_Sedma_davka | seventh dose vaccine type code |
| PrimPricinaHospCOVID | primary reason for COVID hospitalization |
| bin_Hospitalizace | was hospitalized |
| min_Hospitalizace | earliest week of hospitalization |
| dni_Hospitalizace | days of hospitalization |
| max_Hospitalizace | last week of hospitalization |
| bin_JIP | admitted to ICU |
| min_JIP | earliest week of ICU stay |
| dni_JIP | days in ICU |
| max_JIP | last week of ICU stay |
| bin_STAN | hospitalized at regular hospital ward |
| min_STAN | earliest week of regular ward stay |
| dni_STAN | days in regular ward |
| max_STAN | last week of regular ward stay |
| bin_Kyslik | received oxygen therapy |
| min_Kyslik | earliest week of oxygen therapy |
| dni_Kyslik | days of oxygen therapy |
| max_Kyslik | last week of oxygen therapy |
| bin_HFNO | received high-flow nasal oxygen therapy |
| min_HFNO | earliest week of HFNO |
| dni_HFNO | days of HFNO |
| max_HFNO | last week of HFNO |
| bin_UPV_ECMO | received ultra-protective ventilation/ECMO |
| min_UPV_ECMO | earliest week of UPV/ECMO |
| dni_UPV_ECMO | days of UPV/ECMO |
| max_UPV_ECMO | last week of UPV/ECMO |
| Mutace | variant |
| DatumUmrtiLPZ | date of death |
| Long_COVID | long COVID diagnosis date |
| DCCI | Deyo-Charlson Comorbidity Index |
The Infekce column indicates whether the entry is for a person with no infection listed, or for the first, second, or third infection of the person. So part of the reason why there are more rows than people is probably because there's multiple rows for people with 2 or 3 infections. But I don't know why the maximum infection number is only 3 though:
> t[,.N,Infekce] Infekce N 1: NA 7800263 2: 1 4325706 3: 2 450511 4: 3 21188
There's about 1.6 million rows where the year of birth is missing. All except one of them are for people with no vaccine doses listed:
> t[RokNarozeni=="-",.N,OckovaciLatkaKod_Prvni_davka][N>0][order(-N)] OckovaciLatkaKod_Prvni_davka N 1: 1630991 # about 1.6 million rows for unvaccinated people 2: CO01 1 # one row for Pfizer
On most rows with missing year of birth, all other fields including the date of death are also missing, but there's 63695 rows where the year of birth is missing but the date of death is not missing (which are mostly in 2020 like the blog post by SMIS said):
> t=fread("Otevrena-data-NR-26-30-COVID-19-prehled-populace-2024-01.csv")
> t[RokNarozeni=="-"&DatumUmrtiLPZ!="",.N,.(year=substr(DatumUmrtiLPZ,1,4))][order(year)]
year N
1: 2020 37471
2: 2021 9658
3: 2022 5522
4: 2023 6981
5: 2024 4063
If rows with a missing date of birth are excluded and rows for second and third infections are excluded, then the new dataset has only about 10.5 million rows, which is about half a million less than the record-level dataset:
> t[RokNarozeni!="-"&!(!is.na(Infekce)&Infekce>1),.N] [1] 10494977
The new UZIS dataset is missing about 15% of all deaths in 2020 even if you don't exclude the rows with a date of birth missing: [https://ec.europa.eu/eurostat/data/database, czech2.html#Deaths_by_ICD_code_region_age_group_and_year]
> rec=fread("CR_records.csv")
> reclev=rec[,.(reclev=.N),.(year=year(DatumUmrti))]
> eurostat=fread("http://sars2.net/f/czpopdead.csv")[,.(eurostat=sum(dead)),year]
> icd=fread("http://sars2.net/f/czicd.csv.gz")[,.(byicd=sum(dead)),year]
> me=merge(merge(reclev,eurostat),icd)
> ua=\(x,y,...){u=unique(x);y(u,...)[match(x,u)]}
> new=t[!(!is.na(Infekce)&Infekce>1),.(newuzis=.N),.(year=as.integer(substr(DatumUmrtiLPZ,1,4)))]
> merge(me,new,all=T)[year>=2020]|>print(r=F)
year reclev eurostat byicd newuzis
2020 129289 129289 129289 110177
2021 139890 139891 139891 139687
2022 120219 120219 120219 121464
2023 NA NA NA 113807
2024 NA NA NA 70941
However it's because the new dataset is missing all deaths in the first 9 weeks of 2020:
> t[,.N,DatumUmrtiLPZ][order(DatumUmrtiLPZ)][1:6] DatumUmrtiLPZ N 1: 12034894 2: 2020-10 2277 3: 2020-11 2289 4: 2020-12 2318 5: 2020-13 2290 6: 2020-14 2322
The column for the Deyo-Charlson Comorbidity Index shows a comorbidity level expressed as an integer between 0 and 5, where the value is 0 for about 87% of rows:
> t[,.N,DCCI]|>print(r=F)
DCCI N
0 10961996
1 876443
2 347656
3 169019
4 92154
5 150400
The average comorbidity index was about 0.31 for people with first dose type Comirnaty but about 0.38 for people with first dose type SPIKEVAX:
> rec=fread("Czech/data/CR_records.csv")
> k=grep("OckovaciLatka_",names(rec));k2=k-1
> codes=rec[,.(type=unlist(rec[,..k],,F),code=unlist(rec[,..k2],,F))][!duplicated(code)]
> a=t[,.(DCCI=mean(DCCI),.N),.(code=OckovaciLatkaKod_Prvni_davka)]
> merge(codes,a)[order(-N)][N>=1e3]
code type DCCI N
1: CO01 Comirnaty 0.3055174 5820268 # Pfizer
2: 0.1620386 5330971 # unvaccinated
3: CO02 SPIKEVAX 0.3800446 546257 # Moderna
4: CO03 VAXZEVRIA 0.5026497 460052 # AstraZeneca
5: CO04 COVID-19 Vaccine Janssen 0.2669109 429960
6: CO07 Nuvaxovid 0.3323160 5898 # Novavax
7: CO20 Comirnaty Omicron XBB.1.5. 0.5976789 2068
The comorbidity index peaks in people who were born in 1925-1929, so this new dataset might reflect the actual level of comorbidities better than the earlier file ockovani-profese.csv (where the proportion of people who were indicated to have a chronic illness peaked in ages 65-69 but it was very low in the oldest age groups, which might be if information about chronic illness was not listed for people who qualified to be vaccinated based on their age):
> t[,.(DCCI=mean(DCCI),.N),RokNarozeni][order(RokNarozeni)]|>print(r=F)
RokNarozeni DCCI N
- 0.0002066227 1630992
1860-1864 0.0000000000 1
1865-1869 0.0000000000 2
1870-1874 0.0000000000 1
1875-1879 0.0000000000 2
1880-1884 0.0000000000 2
1890-1894 0.0000000000 1
1895-1899 0.0000000000 1
1900-1904 0.0000000000 69
1905-1909 0.0000000000 45
1910-1914 0.0579710145 69
1915-1919 0.7040000000 375
1920-1924 0.9064161207 8089
1925-1929 0.9732974675 45726 # highest comorbidity index
1930-1934 0.9454279068 132412
1935-1939 0.8325555678 226075
1940-1944 0.7322019785 401112
1945-1949 0.6178468076 596779
1950-1954 0.5035088189 672876
1955-1959 0.4384776532 648236
1960-1964 0.4118506011 665654
1965-1969 0.3435418951 706825
1970-1974 0.2727757004 869150
1975-1979 0.2238727484 956530
1980-1984 0.1820616173 781988
1985-1989 0.1483075026 739558
1990-1994 0.1286849167 693842
1995-1999 0.1152187225 525826
2000-2004 0.1038220172 511510
2005-2009 0.0908365436 587704
2010-2014 0.0850617850 560087
2015-2019 0.0456251469 472393
2020-2024 0.0141874725 163736
RokNarozeni DCCI N
Next I calculated an age-normalized excess comorbidity percentage for each vaccine type. I excluded people with a year of birth missing but I didn't exclude rows for second and third infections. I derived the expected total DCCI for each vaccine type by multiplying the number of people in each age group with the average DCCI index for the age group among all rows where the year of birth was not missing, and then I divided the actual total DCCI with the expected total DCCI, subtracted 1, and multiplied the result by 100. The excess comorbidity percentage was now about -4.4% for Pfizer but -2.9% for Moderna, so the difference is not big enough to explain why Moderna would have about 30% higher mortality rate than Pfizer:
> dcci=t[RokNarozeni!="-",.(dcci=sum(as.numeric(DCCI)),.N),.(code=OckovaciLatkaKod_Prvni_davka,born=RokNarozeni)]
> dcci=merge(dcci[,.(base=sum(dcci)/sum(N)),born],dcci)[,base:=base*N]
> dcci=dcci[,.(excesspct=round((sum(dcci)/sum(base)-1)*100,1),N=sum(N)),code]
> merge(codes,dcci)[N>=1e3][order(-N)]|>print(r=F)
code type excesspct N
CO01 Comirnaty -4.4 5820267 # Pfizer
16.0 3699980 # unvaccinated
CO02 SPIKEVAX -2.9 546257 # Moderna
CO03 VAXZEVRIA -8.7 460052 # AstraZeneca
CO04 COVID-19 Vaccine Janssen -7.4 429960
CO07 Nuvaxovid 42.3 5898 # Novavax
CO20 Comirnaty Omicron XBB.1.5. 35.9 2068
When I used a similar method to calculate an age-normalized excess hospitalization percentage by the type of the first vaccine dose, Moderna was about 6 percentage points higher than Pfizer (but my calculation was not adjusted for exposure time or varying hospitalization rates over time):
> ho=t[RokNarozeni!="-",.(hosp=sum(bin_Hospitalizace,na.rm=T),.N),.(code=OckovaciLatkaKod_Prvni_davka,born=RokNarozeni)] > ho=merge(ho[,.(base=sum(hosp)/sum(N)),born],ho)[,base:=base*N] > merge(codes,ho[,.(excesspct=round((sum(hosp)/sum(base)-1)*100,1),N=sum(N)),code])[order(-N)][N>=1e3] code type excesspct N 1: CO01 Comirnaty -27.4 5820268 # Pfizer 2: 91.6 3699980 # unvaccinated 3: CO02 SPIKEVAX -21.7 546257 # Moderna 4: CO03 VAXZEVRIA -28.0 460052 # AstraZeneca 5: CO04 COVID-19 Vaccine Janssen 21.9 429960 6: CO07 Nuvaxovid 222.9 5898 # Novavax 7: CO20 Comirnaty Omicron XBB.1.5. 145.8 2068
Based on the new CSV file published by UZIS in November 2024, it's now possible to calculate hospitalization rate by vaccine type.
In the following code in order to calculate an age-standardized hospitalization rate, I subtracted the start of the 5-year ranges of year of birth from 2020 to assign people to 5-year age groups based on their age on December 31st 2020 (because for example all people born in 1950-1954 belong to the age group 70-74 on December 31st 2020). I didn't bother accounting for the aging of people over time.
The new UZIS data is missing the year of birth of many unvaccinated people, even though it's not missing for any unvaccinated people with a date of hospitalization included. But I think the new data doesn't make it possible to see how many unvaccinated people there are for each combination of month and age, so I took the monthly number of people by vaccine type and age from one of my old bucket files for the record-level data, where I treated the number of vaccinated people in each age group at the end of 2022 as the number of vaccinated people for the rest of the dataset. (In my bucket file I accounted for the aging of people over time, but it shouldn't make much difference here.)
I assigned each hospitalization week to a hospitalization month based on the month of the Thursday of the week.
But anyway, the results below show that people who got a Janssen or AstraZeneca vaccine for the first dose had a high hospitalization rate during the COVID wave in November to December 2021. So it's consistent with my plot for all-cause ASMR by vaccination status where the spike in deaths during the last two months of 2021 was more pronounced for Janssen and AstraZeneca but more flat for Pfizer and Moderna.
In January 2021 vaccinated people had almost as high hospitalization rate as unvaccinated people, which might be because there had not been enough time for the immune response to build up after vaccination, or because vulnerable groups of people were priorized during the early rollout.
library(data.table);library(ggplot2)
t=fread("Otevrena-data-NR-26-30-COVID-19-prehled-populace-2024-01.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}
getweek=\(x)isoweek(as.integer(substr(x,1,4)),as.integer(substr(x,6,7)))
ua=\(x,y,...){u=unique(x);y(u,...)[match(x,u)]}
yemo=\(x){u=unique(x);format(u,"%Y-%m")[match(x,u)]}
t[min_Hospitalizace!="",hospweek:=ua(min_Hospitalizace,getweek)-3]
t[,date1:=ua(Datum_Prvni_davka,getweek)-3]
t[,type:=OckovaciLatkaKod_Prvni_davka]
t[type%in%c("CO01","CO08","CO09","CO16","CO20","CO21"),type:="Pfizer"]
t[type%in%c("CO02","CO12","CO15"),type:="Moderna"]
t[type=="CO03",type:="AstraZeneca"]
t[type=="CO04",type:="Janssen"]
t[type=="",type:="Unvaccinated"]
t[,type:=ifelse(!is.na(date1)&date1>hospweek,"Unvaccinated",type)]
t[RokNarozeni!="-",age:=pmin(100,(2020-as.integer(substr(RokNarozeni,1,4))))]
hosp=t[!type%like%"^CO"&!is.na(hospweek),.(hosp=.N),.(month=yemo(hospweek),type,age=ifelse(age<20,0,age))]
b=fread("http://sars2.net/f/czbucketskeep.csv.gz")
a=b[dose<=1,.(pop=sum(alive)),.(month,type=ifelse(type=="","Unvaccinated",type),age=ifelse(age<20,0,age%/%5*5))]
a=rbind(a,a[month=="2022-12",.(month=yemo(seq(as.Date("2023-1-1"),as.Date("2024-10-1"),"month")),pop),.(type,age)])
me=merge(a,hosp,all=T)[!type%in%c("Novavax","Other")]
me=merge(fread("http://sars2.net/f/czcensus2021pop.csv")[,.(std=sum(pop)),.(age=age%/%5*5)][,std:=std/sum(std)],me)
p=me[,.(hosp=sum(hosp/pop*std*365e5,na.rm=T)),.(month,type)][order(-type)]
p[,month:=as.Date(paste0(month,"-1"))]
p[,type:=factor(type,c("Unvaccinated","Pfizer","Moderna","AstraZeneca","Janssen"))]
xstart=as.Date("2020-1-1");xend=as.Date("2025-1-1")
xbreak=seq(xstart,xend,"6 month");xlab=c(rbind("",2020:2024),"")
ybreak=pretty(p$hosp,7);ystart=0;yend=max(ybreak)
color=c("black",hsv(c(240,355,330,204,36)/360,c(.94,.75,.8,.67,.64),c(.76,.82,.5,.87,.95)))
ggplot(p,aes(x=month+15,y=hosp))+
geom_vline(xintercept=seq(xstart,xend,"3 month"),color="gray85",linewidth=.3)+
geom_vline(xintercept=seq(xstart,xend,"year"),color="gray65",linewidth=.3)+
geom_vline(xintercept=c(xstart,xend),linewidth=.3,lineend="square")+
geom_hline(yintercept=c(ystart,0,yend),linewidth=.3,lineend="square")+
geom_line(aes(color=type),linewidth=.4)+
geom_point(aes(color=type),size=.5)+
labs(title="Age-standardized hospitalization rate per 100,000 person-years",subtitle=paste0("The 2021 census population by 5-year age groups was used as the standard population, where ages below 20 were aggregated to a single age group.")|>stringr::str_wrap(80),x=NULL,y=NULL,color="First dose type")+
scale_x_date(limits=c(xstart,xend),breaks=xbreak,labels=xlab)+
scale_y_continuous(limits=c(ystart,yend),breaks=ybreak)+
scale_color_manual(values=color)+
coord_cartesian(clip="off",expand=F)+
theme(axis.text=element_text(size=7,color="black"),
axis.ticks=element_line(linewidth=.3),
axis.ticks.x=element_line(color=alpha("black",c(1,0))),
axis.ticks.length=unit(3,"pt"),
legend.background=element_rect(fill="white",color="black",linewidth=.3),
legend.direction="vertical",
legend.justification=c(1,1),
legend.key.height=unit(10,"pt"),
legend.key.width=unit(15,"pt"),
legend.key=element_blank(),
legend.margin=margin(3,5,3,4),
legend.position=c(1,1),
legend.spacing.x=unit(1.5,"pt"),
legend.spacing.y=unit(0,"pt"),
legend.text=element_text(size=7),
legend.title=element_text(size=7,face=2,margin=margin(2,,3)),
panel.background=element_blank(),
plot.margin=margin(4,4,4,4),
plot.subtitle=element_text(size=6.7,margin=margin(,,4)),
plot.title=element_text(size=7.5,face="bold",margin=margin(1,,3)))
ggsave("1.png",width=4,height=2.2,dpi=380*4)
system("magick 1.png -resize 25% 1.png")
Here's my old plot for ASMR by vaccine type for comparison, where you can see that the spike in ASMR in the last two months of 2021 is more pronounced for Janssen and AstraZeneca which are viral vector vaccines but the spike is more flat for Moderna and Pfizer which are mRNA vaccines (even though the difference might not necessarily be due to efficacy, but it might also be if people with a higher COVID mortality risk were more likely to get certain types of vaccines): [czech2.html#ASMR_by_month_and_vaccine_type]
