#This script is made freely available under license CC-BY-4.0 
#by Prof. Patrick E Meyer. Belgium
#
#
#Data of 2022 analysis - R script. Patrick E. Meyer. October 2023.
#
#reference:
#All-cause Mortality During Covid-19 Vaccinations in European Active Populations (Meyer et al. 2023)
#
#It produces a paired t.test based table and also a spike concordance table

############################################################################
##################### Analysis  ############################################
############################################################################

folder<-"~/covid/"
datafile<-c("data0to14.csv","data15to44.csv","data45to64.csv")
duration<-11 #number of 4-weeks periods considered for the statistics post-injection

####################### Functions ##########################################

#function to get injection starts in one country (when prop of the population already got a shot)
InjectionStart<-function(country,prop=0.05,dataset) { 
  indexwhere<-which(dataset$Where==country)
  indexwhere<-rev(indexwhere)
  start<-which(cumsum(dataset$DoseRate[indexwhere])>=prop)[1]
  dataset$When[indexwhere[start]]
}

#function to get the rows of the dataset corresponding to all the periods investigated in one country
EvaluationPeriod<-function(country,startingpoints,studylength,dataset) {
  start<-which(levels(dataset$When)==startingpoints[country])
  end<-start+studylength
  index<-which(dataset$Where==country & as.numeric(dataset$When)>=start & as.numeric(dataset$When)<end)
  index
}

#function to get the top k maximum values reached during a period  
#starting at startingpoints and for duration of studylength periods (i.e. 11 periods of 4 weeks)
kMaxVar<-function(country,startingpoints,studylength,variable=1,k=2,dataset){
  index<-EvaluationPeriod(country,startingpoints,studylength,dataset)
  ordermax<-order(dataset[index,variable],decreasing=T)
  dataset$When[index[ordermax[1:k]]]
}

#function to get p-values associated to the difference between two variables over a given set of periods
#starting at startingpoints and for a duration of studylength periods
signifgreater<-function(country,startingpoints,studylength=10,variable1=1,variable2=6,dataset){
  index<-EvaluationPeriod(country,startingpoints,studylength,dataset)
  signif<-t.test(dataset[index,variable1],dataset[index,variable2],paired=T) 
  signif$p.value
}


#function to get the average value of a variable over a period
#starting at startingpoints and for a duration of studylength periods 
#in practice the average excess deaths zscore at max injections period (if studylength = 1)
varAvg<-function(country,startingpoints,studylength,variable1=1,dataset){
  index<-EvaluationPeriod(country,startingpoints,studylength,dataset)
  avg<-mean(dataset[index,variable1])
  avg
}

#function to get the average difference between two variables over a set of periods
#starting at startingpoints and for a duration of studylength periods 
#in practice the average difference between zscores of two periods
varDiff<-function(country,startingpoints,studylength,variable1=1,variable2,dataset){
  index<-EvaluationPeriod(country,startingpoints,studylength,dataset)
  avg<-mean(dataset[index,variable1]-dataset[index,variable2])
  avg
}

#Maximal set of countries given by intersection of data sources
ISOcountries=c("AT","BE","CY","DK","EE", "EL","ES","FI","FR","HU","IE","IT","IS","LU","MT","NO","PT","SI","SE")
names(ISOcountries)<-c("Austria","Belgium","Cyprus","Denmark","Estonia","Greece","Spain","Finland","France","Hungary",
                       "Ireland","Italy","Israel","Luxembourg","Malta","Norway","Portugal","Slovenia","Sweden")

############################ ANALYSIS ######################################

read.csv(paste(folder,datafile[1],sep=""))->dataset1
read.csv(paste(folder,datafile[2],sep=""))->dataset2
read.csv(paste(folder,datafile[3],sep=""))->dataset3

countries<-levels(dataset2$Where)

############################## extracting key dates by countries #####################################
injectionStarts1<-sapply(countries,InjectionStart,0.05,dataset1)
injectionStarts2<-sapply(countries,InjectionStart,0.05,dataset2)
injectionStarts3<-sapply(countries,InjectionStart,0.05,dataset3)

maxinjections1<-sapply(countries,kMaxVar,injectionStarts1,duration,7,1,dataset1)
maxinjections2<-sapply(countries,kMaxVar,injectionStarts2,duration,7,1,dataset2)
maxinjections3<-sapply(countries,kMaxVar,injectionStarts3,duration,7,1,dataset3)

index1<-sapply(countries,EvaluationPeriod,injectionStarts1,duration,dataset1)
index2<-sapply(countries,EvaluationPeriod,injectionStarts2,duration,dataset2)
index3<-sapply(countries,EvaluationPeriod,injectionStarts3,duration,dataset3)

############################# Table of p-values of zscores of 2021 vs 2020 #############################
pval<-c(t.test(dataset1$ZscoreCurrent[index1],dataset1$ZscorePast1Y[index1],paired=T)[[3]],
        t.test(dataset2$ZscoreCurrent[index2],dataset2$ZscorePast1Y[index2],paired=T)[[3]],
        t.test(dataset3$ZscoreCurrent[index3],dataset3$ZscorePast1Y[index3],paired=T)[[3]])

avg1<-sapply(countries,varDiff,injectionStarts1,duration,1,6,dataset1)
avg2<-sapply(countries,varDiff,injectionStarts2,duration,1,6,dataset2)
avg3<-sapply(countries,varDiff,injectionStarts3,duration,1,6,dataset3)
avgmat<-rbind(avg1,avg2,avg3)

signifs1<-sapply(countries,signifgreater,injectionStarts1,duration,1,6,dataset1)
signifs2<-sapply(countries,signifgreater,injectionStarts2,duration,1,6,dataset2)
signifs3<-sapply(countries,signifgreater,injectionStarts3,duration,1,6,dataset3)
matsignif<-rbind(signifs1,signifs2,signifs3)

matsignif[avgmat<=0]<-0
matsignif[avgmat>0 & matsignif<=0.05]<-2 
matsignif[avgmat>0 & matsignif<2]<-1

############################# Table of p-values of zscores of 2021 vs 2019 ################################

pvalb<-c(t.test(dataset1$ZscoreCurrent[index1],dataset1$ZscorePast2Y[index1],paired=T)[[3]],
        t.test(dataset2$ZscoreCurrent[index2],dataset2$ZscorePast2Y[index2],paired=T)[[3]],
        t.test(dataset3$ZscoreCurrent[index3],dataset3$ZscorePast2Y[index3],paired=T)[[3]])

avg1b<-sapply(countries,varDiff,injectionStarts1,duration,1,5,dataset1)
avg2b<-sapply(countries,varDiff,injectionStarts2,duration,1,5,dataset2)
avg3b<-sapply(countries,varDiff,injectionStarts3,duration,1,5,dataset3)
avgmatb<-rbind(avg1b,avg2b,avg3b)

signifs1b<-sapply(countries,signifgreater,injectionStarts1,duration,1,5,dataset1)
signifs2b<-sapply(countries,signifgreater,injectionStarts2,duration,1,5,dataset2)
signifs3b<-sapply(countries,signifgreater,injectionStarts3,duration,1,5,dataset3)
matsignifb<-rbind(signifs1b,signifs2b,signifs3b)

matsignifb[avgmat<=0]<-0
matsignifb[avgmat>0 & matsignifb<=0.05]<-2 
matsignifb[avgmat>0 & matsignifb<2]<-1


########################### Table of spike concordance (significant zscores at max injections) #################


nbperiods<-1 #to have the zscore at max injection
zscoreavg1<-sapply(countries,varAvg,maxinjections1,1,nbperiods,dataset1)
zscoreavg2<-sapply(countries,varAvg,maxinjections2,1,nbperiods,dataset2)
zscoreavg3<-sapply(countries,varAvg,maxinjections3,1,nbperiods,dataset3)
matzscores<-rbind(zscoreavg1,zscoreavg2,zscoreavg3)
threshold<-0.675 #given by qnorm(0.75,0,1) 
matzscores[matzscores>=threshold]<-2
matzscores[matzscores<=0]<-0
matzscores[matzscores<2 & matzscores>0]<-1

pval2<-c(binom.test(sum(matzscores[1,]==2),length(countries),0.25)[[3]],
         binom.test(sum(matzscores[2,]==2),length(countries),0.25)[[3]],
         binom.test(sum(matzscores[3,]==2),length(countries),0.25)[[3]])


################### Display Results ###############################################
colnames(matsignif)<-ISOcountries[countries]
rownames(matsignif)<-c("0-14","15-44","45-64")
colnames(matsignifb)<-colnames(matsignif)
rownames(matsignifb)<-rownames(matsignif)
colnames(matzscores)<-colnames(matsignif)
rownames(matzscores)<-rownames(matsignif)
matzscores<-cbind(matzscores,total.pval=pval2)
matsignif<-cbind(matsignif,total.pval= pval)
matsignifb<-cbind(matsignifb,total.pval= pvalb)
ord<-order(colnames(matsignif))

#print("Statistically significant excess death (with symbol 2) in the ten months post injection compared to same period two years before")
#print(matsignifb[,ord])

print("Statistically significant excess death (with symbol 2) in the ten months post injection compared to same period one year before")
print(matsignif[,ord])
print("Statistically significant concordance (with symbol 2) between unusually high excess deaths and maximal injection rate")
print(matzscores[,ord])