Causal, Bayesian, & non-parametric modeling of the SARS-CoV-2 viral load distribution vs. patient’s age

The viral load of patients infected with SARS-CoV-2 varies on logarithmic scales and possibly with age. Controversial claims have been made in the literature regarding whether the viral load distribution actually depends on the age of the patients. Such a dependence would have implications for the COVID-19 spreading mechanism, the age-dependent immune system reaction, and thus for policymaking. We hereby develop a method to analyze viral-load distribution data as a function of the patients’ age within a flexible, non-parametric, hierarchical, Bayesian, and causal model. The causal nature of the developed reconstruction additionally allows to test for bias in the data. This could be due to, e.g., bias in patient-testing and data collection or systematic errors in the measurement of the viral load. We perform these tests by calculating the Bayesian evidence for each implied possible causal direction. The possibility of testing for bias in data collection and identifying causal directions can be very useful in other contexts as well. For this reason we make our model freely available. When applied to publicly available age and SARS-CoV-2 viral load data, we find a statistically significant increase in the viral load with age, but only for one of the two analyzed datasets. If we consider this dataset, and based on the current understanding of viral load’s impact on patients’ infectivity, we expect a non-negligible difference in the infectivity of different age groups. This difference is nonetheless too small to justify considering any age group as noninfectious.

learned (probability) density's correlation structure and are thus prone to reconstructing 13 unphysical densities. Other methods make use of neural networks (see for example Liu 14 et al. [4]) or restrict the density to specific functional forms (see e.g. Dirichlet Process 15 Mixture Model, or DPMM [5][6][7]). Another approach is to use smooth priors and infer 16 the level of smoothness of the reconstructed density from data via maximum entropy 17 (Density Estimation using Field Theory, or DEFT [8,9]). Chen, Tareen, and Kinney [1] 18 propose an interesting information-theoretical-based modification of DEFT, although it 19 only effectively works in one dimension. Finally, another very commonly used and 20 effective solution to the density estimation problem is the one given by deep neural 21 networks. In a recent work, Liu et al. [4] propose a generative adversarial networks 22 (GAN) which is particularly effective in high dimensions. We refer to their work for 23 comparison with similar neural-network-based approaches.

24
Most of these approaches lack a robust estimate for uncertainties or specify none at 25 all. In the hope of addressing these shortcomings, we propose our novel general-purpose 26 density reconstruction method, which we will refer to as Matérn Kernel Density 27 Estimator (MKDE). This method is very general, works for a generic n-dimensional 28 space, and therefore applies to many different contexts and fields. We presented one 29 paradigmatic example of its many possible applications in Subsec. Causal Structure. In 30 this example, we used MKDE to reconstruct the continuous distributions of the ages 31 and viral loads of patients infected with Covid-19 from age-and-viral-load data samples 32 collected within the general population.  MKDE is capable of reconstructing a smooth density distribution underlying an -34 even limited -discrete dataset. We achieve this result under the hypotheses that the 35 data points are drawn from the underlying density through a Poisson process and that 36 the reconstructed density is a sufficiently smooth function. Since we expect the inferred 37 (probability) density to be strictly positive and to vary on logarithmic scales, we choose 38 the log-normal model where s(x) is the natural logarithm of the density. In the case of multi-dimensional and statistically homogeneous process, and in particular for the Gaussian process s.

57
In order to draw prior samples from the Gaussian field s, we choose a standardized 58 coordinate systemξ. We then transform the standard normally distributed parameters 59 ξ = (ξ k ) k , where k ∈ N is the Fourier-space index according to the mapping 60 s = F −1 A ξ with A := diag( P s ).
September 21, 2022 2/4 Here, F represents the Fourier transform operator and A the amplitude operator in 61 Fourier space. The amplitude operator encodes the Matérn-kernel correlation structure, 62 parametrized with where a s is a scale factor which accounts for the standard deviation in position space, k 0 64 is the magnitude of the characteristic correlation-length wavevector, and γ s is the 65 spectral index of the power spectrum. We assume a s and k 0 to be a priori log-normally 66 distributed since we expect strictly positive variations of the possible power spectra on a 67 logarithmic scale. Similarly, we choose γ s to be normally distributed, since the spectral 68 index could in principle also be negative. We additionally introduce volume factors to 69 ensure that the model parameters are intensive with respect to volume, i.e. they do not 70 depend on the volume in position space.

71
For higher-dimensional data, we expect the correlation structure along each axis (or 72 dimension) to be a priori independent from the others. These different axes could in uniformly-distributed parameter α. For more details on the zero mode degeneracy and 83 factorizing power spectra, we refer to Arras et al. [11].

84
At this stage, we can summarize all the parameters that we have introduced for each 85 independent axis with the scalar-valued parameters α, a s , k 0 , and γ s and the 86 vector-valued ξ k . We set broad priors on these parameters and learn them using MGVI. 87 We set the following priors on the signal parameters: α = [10 −15 , 5.0], a s = (0.3 ± 0.2), 88 k 0 = (4.0 ± 3.0), and γ s = (−6.0 ± 3.0), where the mean and standard deviation specify 89 a Gaussian prior distribution for γ s and log-normal distributions with the given mean 90 and standard deviation for a s and k 0 . For details on the inference of posterior estimates 91 for the MKDE parameters through MGVI and their uncertainty quantification, we refer 92 to Subsec. Inference Fig. 1 illustrates the performance of MKDE in a two dimensional 93 setting.

94
In conclusion, we described MKDE, a Matérn-kernel-based, Bayesian, and 95 non-parametric density estimator that can construct a smooth (probability) density 96 function from an -even limited -set of data samples. The broad priors on the learned 97 parameters, combined with the log-normal model and the Matérn kernel covariance 98 structure, make MKDE very flexible and robust. Furthermore, the Bayesian inference 99 framework allows for posterior uncertainty quantification for the reconstructed density. 100 A software implementation of MKDE is available in NIFTy 7 and is also released as an 101 open-source Python package (DENSe), which can be found at: