key: cord-0918664-3hm0ulno
authors: Finlay, C.; Bassett, B. A.
title: Climate & BCG: Effects on COVID-19 DeathGrowth Rates
date: 2020-07-15
journal: nan
DOI: 10.1101/2020.07.13.20152991
sha: 839d359df94e124d6c39bf73c8f0ea8d3f79ab24
doc_id: 918664
cord_uid: 3hm0ulno

Multiple studies have suggested the spread of COVID-19 is affected by factors such as climate, BCG vaccinations, pollution and blood type. We perform a joint study of these factors using the death growth rates of 40 regions worldwide with both machine learning and Bayesian methods. We find weak, non-significant (< 3 sigma) evidence for temperature and relative humidity as factors in the spread of COVID-19 but little or no evidence for BCG vaccination prevalence or PM2.5 pollution. The only variable detected at a statistically significant level (> 3 sigma) is the rate of positive COVID-19 tests, with higher positive rates correlating with higher daily growth of deaths.

The COVID-19 pandemic has triggered extensive efforts to predict the severity of COVID-19 to aid in decision making around interventions such as lockdown and the closure of schools 1 [1] . Regions hit hard by the pandemic, such as Wuhan, Lombardy and New York, where doubling times of 2-3 days and high crude mortality rates stand in stark contrast to other countries only mildly affected such as Hong Kong, South Korea and New Zealand.

Potential explanations for the apparent differences in the transmissivity (encoded by the time-dependent reproductive number, R t ) and lethality (encoded by the Infection Fatality Rate, IFR), currently fall into four broad categories. The first posits that differences are largely fictitious, driven by the heterogeneity of testing and reporting of cases and deaths; a known issue and one that we are particularly concerned with in this paper.

The next three categories posit that the differences are primarily real and are driven by (1) cultural and policy factors (swift lockdown, efficient contact tracing and quarantining, use of masks, obedient populations or social structures that are naturally distant or isolated), (2) local environmental factors (such as temperature and humidity variations, population density, comorbid factors, vaccinations, vitamin D levels, blood type etc...) or (3) existence of multiple strains with different transmissibility or lethality [4] . Our primary interest lies in disentangling the first category (testing) from a subset of potential factors in the second category.

Finding the relative contributions of each of these four categories is key in understanding and optimally fighting the pandemic. The widely different testing capabilities between countries, particularly between the developed and developing world, imply that if not correctly treated, testing variability will create spurious evidence that can lead to false hope and sub-optimal interventions.

In the wake of the spreading pandemic there have been a host of studies that have examined the possible impact of climate [8]- [20] , blood type [21, 22] , haplogroup [32] , pollution [27] - [31] and BCG vaccination prevalence [23, 24, 25, 26] on the spread of COVID-19. Our main conclusion is that testing issues are significant and the factors above are likely not the main causes of variability in growth rates of deaths worldwide.

We note that both the basic reproductive number, R 0 , and the IFR will be affected by the four categories above in general. However, since we do not currently have access to the true number of infections we cannot address the potential effect of environmental factors on the IFR. Similarly the Case Fatality Rate cannot be used for this purpose due to testing differences around the world. We therefore focus on their potential effect on R 0 , which we quantify through the daily growth rate of deaths of countries around the world for which we have sufficiently good data.

To explore the potential impact of climate (e.g. through temperature, relative humidity and UV Index), BCG vaccination prevalence, blood type and pollution (PM 2.5 ) on COVID-19, one must think carefully about choice of both data and methodology.

We choose to focus on deaths instead of the number of confirmed infections in the belief that these are significantly less affected by sampling and testing issues than infection numbers: the number of tests per 1000 population (Tests/1k) currently varies by more than three orders of magnitude across the world 2 , making infection numbers highly biased and correlated with confounding variables such as GDP and healthcare. In countries with limited testing capability, only the more severe patients are typically tested.

Since these are also the patients most likely to die, it is likely that patients who die from COVID- 19 are more likely to be tested than typical COVID-19 sufferers, who may be mostly asymptomatic. Without detailed knowledge of testing protocols in each country separating out the effects of the factors of interest (climate, BCG vaccination etc...) is extremely challenging. Although deaths are not immune to testing issues, with many missing deaths shown through excess mortality studies, especially in overwhelmed medical systems [48, 49] , we argue that this is likely to be much less of an issue in the first month after initial COVID-19 deaths, which we focus on, since we are interested in R 0 .

The second key choice is whether to focus on absolute death counts or growth rates. Absolute death counts are also highly problematic. First, testing efficiencies may vary systematically from country to country. Second, the widely varying start dates of the pandemic in different countries are hard to deal with rigorously. As a result we focus on the daily growth rate of deaths, G c , for country c during the initial phase of the epidemic. The initial daily growth correlates with R 0c , the base reproductive number in each country; in simple models R 0c = (1 + G c ) τ , where τ is the period in days for which a person is infectious on average 3 . Separating out the true causes of variation in G c or R 0c is still a highly complex problem, and we discuss our approach in detail in the additional material in section (5). Part of the complexity arises from the fact that R 0c depends explicitly both on properties of the virus and on social factors (e.g. average number and nature of contacts), as well as any of the other potential factors we wish to study. Hence, we expect there to be a large intrinsic variation in the growth of deaths from country to country, driven e.g. by the average number of people in a household, population density [5] etc... that may correlate with, and hence lead spurious evidence for, potential factors such as climate, vaccination coverage, blood type and pollution.

We model this complexity by allowing each country, c, to have its own unique base growth rate, denoted G 0c , that is estimated from its own data. This freedom allows the model to account for the myriad unmodelled factors specific to each country (population density, GDP, culture, health care quality etc...). However, we tie these base country growth rates, G 0c , together through a parent distribution, expressing the prior belief that there is a single dominant strain of COVID-19 globally.

Our primary goal in this study is to examine worldwide data to investigate whether climate and the BCG vaccination prevalence are important drivers of COVID-19 spread. Because of the danger of confounding variables we used the base growth rate, G 0c for each country in a machine learning feature extraction algorithm to pick the most important additional variables to include in our computationally intensive hierarchical Bayesian analysis. This lead us to exclude PM 2.5 pollution. In addition we undertook a separate analysis including A+ blood type. The data for blood type came from unpublished online sources and is therefore kept separate since it is less trustworthy.

The remaining four most promising environmental factors were Temperature (T, • C), Relative Humidity (RH, %), BCG vaccination coverage (BCG, %), Ultra-Violet Index (UVIndex); the latter included as a proxy for vitamin D production. To this set we added two testing-related variables to serve as diagnostics for potential contamination from testing issues: (1) the fraction of tests that return positive (Pos-Rate, %) and (2) the number of tests per 1,000 population (Tests/1k), yielding our final set of global parameters, Θ ≡ (T, RH, BCG, UV, Pos-Rate, Tests/1k). Due to concerns over data integrity, we study the impact of A+ blood type separately in section (5.4.1).

We then explicitly fit for Θ using our sample of 40 regions in 37 countries for which there is data for all of the parameters in Θ. We cannot use all the death data for countries since a constant growth rate model fails quickly for most countries due to Non-Pharmaceutical Interventions (NPIs) or nonlinear effects. To address this we extract the initial pure exponential part of the data on a country-by-country basis, as illustrated in Fig. (1) and discussed in detail in section (5), which is Page 2 of 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint the only data we use in our main analysis.

The early phase death data for each country c, D c (t), are then fit to the model:

where β = 10 −2 converts our growths to percentages and the X c are the country-specific data corresponding to Θ. We use a hierarchical Bayesian analysis in which each country's growth rate is drawn from a parent distribution which allows each country's base growth rate, G 0c , to be intrinsically different rather than forcing any differences to be due only to the parameters in the factors encoded by Θ. We use Monte Carlo Markov Chain (MCMC) methods to simultaneously fit the base growth rates, G 0c , for all countries, our parameters of interest, Θ, and the parent distribution hyperpriors, implying a large Bayesian hierarchical model with 90 parameters in total. After marginalising over all country growth rates and parent distribution parameters, we are left with the marginal distributions on the Θ, our parameters of interest, providing our main results.

We assess the importance of each of the Θ parameters both through standard model selection metrics such as the Bayesian evidence and Bayesian Information Criterion (see Table 1 ) and by computing the statistical significance with which the parameters deviate from zero in the marginalised posterior chains. We now discuss these results.

Figure (2) shows the main results of our paper: only the parameter associated to the positive rate or fraction of tests (Pos-Rate) is non-zero at more than 3σ. This result is stable to significant changes in our priors and hyperpriors and to including or removing the other parameters in Θ in the analysis. In addition, the positive rate is selected as the most important variable by all model-selection metrics (see Table 1 ) and by our machine learning analysis, section (5.5). A plot of regional death growth rates versus positive rate is shown in Fig. (3) showing the clear correlation.

The simplest explanation of this result is that regions that experience the most rapid spread of the disease, i.e. those with the largest R 0 , were also the regions that were least able to keep up with testing demands and hence where the rate of Polymerase Chain Reaction (PCR) tests used for COVID-19 returning positive were the highest on average. This is then not a cause, but rather a result of, high growth rate. We find that Pos-Rate correlates positively with Tests per 1000 population (see Fig. 10 ): high death growth rate and growing positive rate likely spurred increased testing. Neither of these observations help us identify the cause of increased death growth rates, however.

Temperature is the only other variable which is nonzero at more than 2σ in our multivariate fits (see section 5.4.3 for more discussion) and relative humidity the only other parameter non-zero at more than 1.5σ, providing some weak evidence for climatic impact on the spread of the disease. Increasing temperature and humidity tends to decrease the spread of COVID-19, in agreement with previous studies.

In our model-selection metrics, where we compare fits with one variable at a time, relative humidity is preferred over temperature by the Akaike Information Criteria (AIC) and Bayesian Information Criteria (BIC) and is the only variable other than Pos-Rate, that has a BIC value more than 2 units lower than the No Factor model (which has Θ = 0) 4 . The remaining parameters have poor or mixed results relative to the No Factor model, and all perform very poorly relative to the model with just the positive rate.

In particular, our regression and machine learning analyses provide no evidence in support of BCG fraction. This is not surprising if we look at figure (4) which plots regional death growth rates versus BCG coverage: there is no discernible correlation. Our results are therefore in agreement with [25, 26] . Further we do not find any correlation between UV index and death growth rate. This is pertinent since UV index is relevant to natural production of vitamin D [39] and vitamin D has been suggested as a protective factor against the spread of the disease [36, 37, 38] . Now we briefly discuss results for the other testing 4 ∆ BIC= 2 − 6 is the standard demarcation of positive evidence [41] .

Page 3 of 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint The results in both cases are consistent: only the coefficient associated with the positive rate of tests is non-zero at more than 3σ, with temperature non-zero at 2σ in the multivariate case. Notice that the significance of both the temperature and relative humidity increase in the multivariate fits relative to the univariate fits where they are fit alone. Both BCG and UV Index are consistent with no effect.

variable, Tests / 1k. We find that this is not a significant correlate of growth rate. This is perhaps not unexpected. If testing efficiency in each country is approximately constant, i.e. the rate of true infections detected changes only slowly with time, then this has little impact on the growth rate of deaths. It is only if the testing efficiency varies rapidly with time that we would expect this to be significant. Our results suggest that this was not the case, at least within the initial phase of the epidemic in each countries. Neither the Bayesian, nor the machine learning analyses found the number of tests per 1000 population to be significant. We present technical details of the data, algorithms and analysis in section (5).

Finally, in our separate analysis of blood types discussed in sections 5.4.1 and 5.5 we analyse the potential correlation of blood type with growth rate. We find marginal evidence for A+ blood type being relevant,

CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint subject to the caveats discussed in section (5.4.1). How should our results be taken in the context of the many claims of climate, blood, BCG etc... being significant factors in COVID-19 spread? First, many studies were based on confirmed COVID-19 test cases which, as we discussed earlier, are affected by differences in testing capability and protocols between countries.

Secondly, many of the studies present regressions that do not allow for unmodelled confounding sources of variation in the growth. Hence, if a country shows high growth the algorithm will try to force one of the potential factors under study to explain it. Instead the hierarchical Bayesian framework allows the base growth rate of each country to be different, and hence potential factors will only be given credit for the difference in the growth rate if they provide a genuinely better fit. Further, many studies do not model the intrinsic uncertainties associated with the data as we have done. We too find that the best-fits are non-zero (as can be seen by looking at the peaks of the posteriors in Figure  ( 2) or at the last rows of Table (6). Hence our results are not in disagreement with regression results, the issue is about the statistical significance of such claims.

Finally, although we do not detect environmental factors at more than 3σ, it is interesting to examine how big the environmental factors would be if our bestfit parameter values describe reality: a 15 • C increase (decrease) in temperature implies a 5.25% decrease (increase) in the base daily growth rate while a increase (decrease) of 20% in relative humidity would mean a 2.24% decrease (increase) in base daily growth rate. The decreased spread at higher temperature and humidity agrees with previous work [9].

This may not seem significant, but for a city such as Johannesburg, where both temperature and humidity drop significantly in winter, the combined effect could add more than 5% to the base daily growth rate. For a daily growth rate of 5%, which was approximately the value in May 2020, this would halve the doubling time of the disease, a significant impact.

Page 5 of 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint

Contrary to previous claims our analysis of growth rates for deaths from countries worldwide are consistent with no effect from climate, pollution or BCG vaccination. The only significant correlation detected is with the positive rate of tests: a country that intrinsically had a high R 0 (due to high population density etc...), would naturally tend to be more overwhelmed and hence run low on testing kits earlier, leading to increased fraction of positive tests. We did find some weak suggestive evidence, at < 3σ, that temperature and relative humidity correlate with death growth rates.

A separate analysis of blood type data shows that A+ type is the most important blood type, though the significance is marginal, both because the data quality is low and the statistical significance is weak. Our combined statistical and machine learning analysis finds no evidence for PM 2.5 pollution, other blood types or UV Index as drivers of COVID-19.

More data could be obtained by dropping the requirement that all countries in the sample have data for all the potential factors, which could potentially allow for some of the effects to be detected at higher statistical significance but at the cost of making model comparison significantly more difficult.

The data and code for our Bayesian analysis is available at https://github.com/chrisfinlay/covid19/.

We thank Michael van Niekerk for extracting data from BCG ATLAS, Niayesh Afshordi, Ewan Cameron, Inger Fabris-Rotelli, Ben Holder and Nadeem Oozeer for discussions and comments. This research has been conducted using resources provided by the United Kingdom Science and Technology Facilities Council (UK STFC) through the Newton Fund and the South African Radio Astronomy Observatory.

We chose to look at the growth rate in deaths as they are less likely than confirmed cases to be affected by the widely varying testing protocols between countries. The cumulative death data for each country comes from the Centre for Systems Science and Engineering (CSSE) at Johns Hopkins University (JHU) [45] .

We cannot, however, take all data on deaths for all times since it is clear that a simple exponential model fails quickly: as soon as interventions occur the simple exponential model fails and we will obtain contaminated estimates of the growth rate due to the flattening of the curve which will skew our analysis; see e.g. Fig.  (1) . As a result we want to know which data should be included in our analysis for each country. This reduces to knowing the time, t 0c , at which the growth of deaths deviates from a simple exponential. We then only use data at t < t 0c for country c.

Doing this cut by hand introduces the possibility of human bias. Instead we compute t 0c for each country by fitting a non-linear sigmoid transition function to the growth rate of the death data of each country over time:

This function allows the growth rate to smoothly transition from an initial value (which we are interested in) to a final value, reflecting the impact of social interventions or nonlinearity in the system. In this first phase of the analysis each country, c, therefore has 5 parameters (D 0c , G 0c , δG c , α c , t 0c ) to describe the trajectory of its deaths, where δG c and α c represent the change in growth rate and suddenness of the transition with time from G 0c to G 0c + δG c for each country, c (the index, c, has been left out in Eq.

(3) below):

To maximise the probability that this is a good fit to the data the non-linear model in Eq. (3) we only consider data up to the time when a country passes 1000 deaths. All five parameters are treated hierarchically with their own parent distribution where the means are Normally distributed and the standard deviation is Half Normal distributed.

We only use data from t < t 0c in our main analysis (see section 5.2), where t 0c is given by the marginalised mean from the chains. Typical values for t 0c were around 15 days. We excluded countries where t 0c < 10 because of concerns about quality of the underlying model fit in such cases, i.e. a simple exponential model Page 6 of 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint was not a good fit even at early times which could lead to spurious growth rates.

The resulting fits to country death data for both the full 5-d model (green) and just the region t < t 0c are given in Figure (9) , showing that the technique performs well in isolating the initial exponential growth from a first-principles approach.

Once we have t 0c for each country we further cut our data by using the following rules:

• We only use countries for which there were a minimum of 20 deaths by 25 April 2020. • We choose day zero as the first day a country passes 3 deaths. This leads to typical starting numbers of deaths of around 5. • Data for our variables of interest (Climate, BCG, etc...) were not available for all countries. To ensure the same amount of data for all models only countries with data for all our parameters in Θ were included.

We did not model the potential correlations and interactions of our parameters with these data cuts. This could potentially alter our conclusions: if a parameter were extremely important in determining growth rates, then countries with very small numbers of deaths would systematically not make it past our data cuts and hence the signal from that data would be lost. However, it is unclear how to model this censorship rigorously and it is left to future work.

After these data cuts we were left with 40 provinces in 37 countries with 613 data points in total, shown in Figure (9) . The only country with more than one province was Canada which include Alberta, Ontario, Quebec and British Columbia (BC).

For each region and country the climate data was gathered from the Dark Sky API, [44] , where the locations sampled are taken from the latitudes and longitudes given for each country or region in the JHU COVID-19 data. For each country, c, the mean temperature, mean relative humidity and mean UV Index are calculated as an average over a N day window starting 28 days prior to day 0 for each region, where N is the number of days of deaths data for that country. The latter is chosen as an estimate of the average time from infection to death 5 . Since mean climate variables change relatively slowly the exact delay is not important. See Fig (9) for the distribution of t 0c times. We used UV Index as a proxy for natural vitamin D production. Data for all regions and countries is shown in Table (5) .

Bacillus Calmette-Guérin (BCG) vaccination policies have been in place in many countries across the world starting from widely differing start dates. We would like to estimate the percentage of the population in each country that has received a BCG vaccine in their lifetime. For this purpose we need the age demographics of each country as well as the dates when BCG vaccinations became/stopped being mandatory.

To estimate the percentage of the population vaccinated by BCG we need to draw on three sources of data. These are the BCG Atlas, World Health Organisation (WHO) BCG vaccination rates amongst 1-year olds and the age demographics for each country from the United Nations (UN). Firstly we looked to the BCG Atlas [46]. This is a heterogeneous dataset which implies that not all information was available for each country. We collected the following fields from the BCG Atlas: We did not use the last field directly due to missing information on this field, specifically if it was for a small age range of the population or its entirety. For some countries data was missing from fields 1-3. In these cases, where appropriate and possible, the missing data was obtained from fields 4 and 5.

Age demographics for each country were used to compute the expected fraction of the population who have had the BCG vaccination. In 15% of cases the BCG Atlas did not have data and we instead used WHO data [6] to perform the estimate. The derived BCG fractions and the origin of the data, are shown in Table  ( 2) while a plot of the BCG fractions versus growth rates are shown in Fig. (4) .

There is one caveat here: our BCG coverage is the estimate of the percentage of the population of each region and country that has had the BCG vaccination. This is arguably not the optimal quantity to use in our analyses however; it might be better to use the fraction of vaccinated population weighted by the probability of infection as a function of age. However, the latter is currently unknown and hard to compute even in the best of situations: how can we know how many people were exposed but never got infected? The large scatter in Figure (4) suggests that this is unlikely to make a significant difference to our conclusions that BCG vaccine is not important in the spread of COVID-19.

Data on prevalence of blood type for each country was taken from [33, 34] while PM 2.5 pollution data came from [35] . We discuss the blood type data and results in section (5.4.1).

Page 7 of 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. 

As described earlier our basic regression model for the deaths in country c at day t is:

which is assumed valid for T ≤ t 0c , as described before. G c is measured in percent and we model it's potential dependence on our variables of interest as:

where G 0c and X c are the country-specific base growth rate and data for climate, BCG etc..., shown in Table  ( 5) and Θ are the global parameters we are interested in. In general the X c could be time-dependent. In this analysis, because of missing data and the fact that our data for each region typically spans a short period (t 0c is less then 3 weeks as discussed in section 5.1.1), we use average values for X c . Our goal is to determine if any of the parameters Θ are rigorously required to be non-zero by the data. One limitation of Eq. (5) is that it is linear in the country data X c . We justify this by noting that the growth rates and the underlying data vary over narrow ranges, so that retaining only the linear terms in the Taylor series expansion of G c (Θ) is a reasonable step. For the temperature variables we have verified that assuming instead a step change in the growth at around Page 8 of 17

. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint 10 • C with two hierarchical growth parameters, one on each side of the step, did not lead to any increase in significance in the detection of temperature as an effect.

We write a hierarchical Bayesian probabilistic model by assuming a Poisson likelihood, suitable for count data, with mean given by the deterministic forward model defined in Eq. (5), together with priors and hyperpriors for all our 90 parameters, as shown in the schematic graphical representation in Fig. (5) ; see e.g. [40] . We do not try to model missing deaths due to testing irregularities. As long as the fraction of missing deaths remains approximately constant in the early phase of the spread within the country that we consider this should have little impact on our results.

Since our data covers 37 countries and 40 provinces worldwide and we know that the growth rate of the disease will depend both on properties of the disease (which are universal), and properties specific to each country (e.g. culture and population density) it is natural to model the data with a hierarchical Bayesian structure which allows growth rates to vary somewhat from country to country but to also be somewhat similar between countries.

The priors and hyperpriors for each variable to be estimated is chosen to be:

• µ D0 ∼ N (7, 3 2 ); the hyperprior for the mean on the prior for initial deaths, D 0c , for each country (D 0c > 3 as part of our data cuts). • σ D0 ∼ Half Normal(3); the standard deviation on the prior for initial deaths for each country. • µ G0 ∼ N (20, 20 2 ); the mean (in percent) on the prior for the growth rate of each country, G 0c . • σ G0 ∼ Half Normal (20) ; the standard deviation on the prior for the growth rate of each country. • D 0c ∼ N (µ D0 , σ 2 D0 ), the prior on the initial deaths parameter for each country, c. All countries together will be denoted by D 0 . • G 0c ∼ N (µ G0 , σ 2 G0 ); the prior on the growth rate parameter for each country, c. All countries together will be denoted by G 0 . • Θ ∼ N (0, I); the prior on the vector of parameters of key interest.

The complete vector of all parameters is therefore:

The joint prior over our full set of parameters is assumed factorizable, i.e. a product of the prior distributions listed above. We assume our data to be statistically independent both between countries and from day-to-day. Our model simultaneously fits the death data from all countries. We use NumPyro [43] for Monte Carlo Markov Chain (MCMC) probabilistic sampling from our posterior using the No U-Turn Sampling (NUTS) algorithm [42] . NumPyro has diagnostic tools built in and allows for easy running on accelerators, such as GPUs, which were key in being able to iterate quickly through our many models including all the multivariate and univariate combinations, as well as looking at the effect of priors.

The MCMC simulations were generally run as 4 independent chains with each chain starting from a random position sampled from the prior for our parameters and run for 2000 steps in order for length scales (and other sampling options) to be determined automatically by NumPyro. After this initial "burn-in" each chain was typically run for 5000 samples. Chains were long enough to ensure that the Gelman-Rubin convergence test was always less than 1.01.

Each chain was then thinned by a factor of two in order to further improve independence of samples. This leads to the final 10k samples collected for each model and led to good, converged, trace plots and posteriors; see Fig. (2). We used ChainConsumer 6 for the analysis of the chains, the model selection metrics in Table ( 1) and some of the plots [52] .

The data for blood type prevalence comes from online heterogeneous collections of published and unpublished data sources covering a wide range of publication dates. As a result, data integrity is an issue and we have chosen to separate the blood analysis from our main Bayesian analysis.

A random forest analysis finds only A+ blood type as potentially relevant amongst the different ABO blood types (see section 5.5). Here we present the results of both univariate and multivariate A+ blood type Bayesian analyses using the formalism described in section (5.2). Fig. (6) shows the result of the Bayesian analysis for the A+ blood type coefficient: it is consistent with zero in both the univariate and multivariate cases. If we look at the AIC, BIC and DIC for A+ blood type we find values of: (10.7, 10.7, 23.4) compared to the no-factor models results of (15.3, 11.6, 22.9). The A+ AIC and BIC are better than the no-factor model, but the DIC value is worse. In addition, the random forest analysis found A+ blood type to be more important than both BCG and and Tests/1k.

As a result we conclude that there is somewhat conflicting evidence regarding the potential effect of A+ blood type, but that none of the evidence is strong. 6 https://samreay.github.io/ChainConsumer/ Page 9 of 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint Figure 6 : The univariate and multivariate results for the coefficient of A+ blood type are consistent with zero, i.e. no effect.

In the results section we focused on the best-fit values of Θ. The reverse question is interesting too: what is the impact of including the Θ parameters on the hyperprior parameters, mean, µ G0 , and standard deviation, σ G0 , of the base growth rates, G 0c , for the 40 regions? If the Θ do explain some of the variation we would expect the standard deviation on the base growth rate (parent distribution) to get smaller, which is captured by the hyperparameter σ G0 . The mean of the base growth rates (parent distribution), given by µ G0 , are nothing more than an offset dictating the value when all factors are zero i.e. Θ = 0. This is what we see: in the No-Factor model σ G0 = 7.9% which changes to 6.3% when we allow all the Θ to vary, a significant shrinkage. On the other hand, if we include all parameters except for the two testing parameters, σ G0 returns to 7.9%, while σ G0 = 6.7%(7.7%) when we add the positive rate (Tests/1k) parameters alone; showing that it is primarily the positive rate parameter that drives the shrinkage in the uncertainty in the parent distribution on the base growth rates.

To assess the stability of our results we fit each of our key parameters in Θ both alone (univariate) and simultaneously with all the other parameters (multivariate). As shown in Table ( 3) and Fig. (2) , the results are quite stable. Positive rate and temperature are the still the most significantly parameters in both cases: the positive rate is non-zero at more than 3σ in both cases while the significance for temperature increases from 1.45σ to 2.19σ when going from the univariate to multivariate fit. Relative Humidity also increases in significance in the multivariate case.

Multivariate Univariate BCG Vaccine 0.01 ± 0.04 0.00 ± 0.04 Temperature -0.29 ± 0.14 -0.19 ± 0.13 Relative Humidity -0.10 ± 0.08 -0.09 ± 0.09 Tests / 1k -0.11 ± 0.12 -0.14 ± 0.11 Positive Rate 0.25 ± 0.08 0.27 ± 0.08 UV Index -0.01 ± 0.95 -0.47 ± 0.89 To access correlations between parameters we show in Fig. (10) the one and two-dimensional marginalised posterior plots. Correlations are typically small, though there are is some small positive correlation between Tests per 1k and Positive Rate, and a small negative correlation between Temperature and UV Index.

Our full results for all the hyperpriors, the base country growth rates, G 0c and the Θ parameters for each of the different univariate and multivariate models are shown in Table (6) .

Since our goal in this analysis is to assess whether there is evidence for external factors such as climate, blood type etc... in addition to known country-specific factors, an important internal check is to look for potential sensitivity to our priors and hyperpriors.

In our hierarchical analysis, data for each country tries to pull the measured growth rates of countries to their own best-fit values, while the hierarchical nature of our model described in section (5.3) attempts to pull them all together. Between this tug of war, the algorithm searches for joint values of Θ that will improve the fits to all the data. One concern might be that if we make the hyperpriors on the parent distributions much stronger or much weaker, we allow the algorithm to given less or more freedom to the base growth rates G 0c which in turn may affect the best-fit parameters Θ and our main conclusions.

To test this we tightened the following hyperpriors and priors listed in section (5.3) by an order of magnitude to:

• µ G0 ∼ N (20, 2 2 ); the mean (in percent) on the prior for the growth rate of each country, G 0c . • σ G0 ∼ Half Normal(2); the standard deviation on the prior for the growth rate of each country. • Θ ∼ N (0, 0.1); the prior on the vector of parameters of key interest.

The results of these changes are shown in Fig. (8) . As expected tightening the prior on the Θ pulled most Page 10 of 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint parameters slightly closer to zero but also tighten the posterior so that none of our conclusions were altered: the statistical significance of all variables was unaltered.

To provide a largely independent test of our Bayesian results we also undertook a machine learning analysis of the deaths data using random forest regression [53] . Random forests are a powerful ensemble method that naturally provide feature selection capability, and are hundreds of times faster to run than our computationally intensive hierarchical Bayesian framework. While random forest does not provide estimates of statistical significance of factors it does allow us to rank additional factors in terms of importance and hence to explore the potential of additional explanatory features for inclusion in the main set of parameters, Θ, for the Bayesian analysis.

To undertake the random forest regression we used the results from our No-Factor Bayesian run (i.e. Θ = 0) to obtain the base growth rates G 0c for each country and used these as the target for the random forest regression with the X c data as features. In the random forest analysis we augmented the data both with PM 2.5 pollution and extra blood types data (namely A-, AB+, O-, AB-, B-).

Feature selection was done using the standard Random Forest impurity method [53] and by ranking variables by the impact they had on the average Root-Mean-Square-Error (RMSE) of the regression: leaving out important explanatory variables is expected to significantly degrade the performance of the algorithm, leaving out irrelevant features does not. All feature importance scores were averages over 500 random 70-30 training-test splits of the data. We found that the only variables that made more than one percent difference to the RMSE value were the Postive Rate (3.98%) and A+ (2.28%). This lead us to select A+ blood type as a variable in our separate full Bayesian analysis shown in section (5.4.1).

The results from the impurity-based feature selection are shown in Table (4). Again Positive Rate and A+ were the most significant features, followed by Tests per 1k and BCG vaccine coverage, PM 2.5 and the other blood types at significantly less importance. As a result blood types, other than A+, were not included in the Bayesian analysis since additional parameters significantly increased the computational complexity of the analysis, both because of increased time to convergence and an increase in the number of models to run (since we fit both multivariate and univariate models in all cases). The distribution of a selection of feature importances is shown in Fig. (7) In summary the machine learning analysis confirms, to the extent that they overlap, the results of our much more intensive Bayesian analysis: Positive Rate is the most important feature, followed by A+ (subject to the . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint the priors are tightened by a factor of ten. We see that the significance of the best-fit changes by less than 0.5σ although the means are typically shifted towards zero, other than BCG and Positive Rate. Our main conclusions are unchanged and stable to changing the priors dramatically.

[6] https://www.who.int/data/gho/data/indicators/indicatordetails/GHO/bcg-immunization-coverage-among-1-year-olds-(-) . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

The copyright holder for this preprint this version posted July 15, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 15, 2020. . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review)

The copyright holder for this preprint this version posted July 15, 2020. . https://doi.org/10.1101/2020.07.13.20152991 doi: medRxiv preprint

Estimating the effects of nonpharmaceutical interventions on COVID-19 in Europe

Proceedings of the National Academy of Sciences

Climate & BCG: Effects on COVID-19 Death Growth Rates

Data analysis using regression and multilevel/hierarchical models

The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo