Causal, Bayesian, & Non-parametric Modeling of the SARS-CoV-2 Viral Load
Distribution vs. Patient’s Age

Matteo Guardiani1,2,*, Philipp Frank1,2, Andrija Kostić1,2, Gordian
Edenhofer1,2, Jakob Roth1,2, Berit Uhlmann4, and Torsten Enßlin1,2,3

1 Max Planck Institute for Astrophysics, Garching, Germany

2 Fakultät für Physik, Ludwig-Maximilians-Universität München, Munich, Germany

3 Excellence Cluster ORIGINS, Garching, Germany

4 Süddeutsche Zeitung, Munich, Germany

*<EMAIL_ADDRESS>
## Abstract

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.

## Introduction

Children do not seem to be major drivers in the transmission of Severe Acute
Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) in the general population [1].
However, the exact degree to which children and adolescents get infected by,
and are able to transmit the virus is not yet well known. Their role in the
community spread depends on their susceptibility, symptoms, viral load, social
contact patterns, behavior, and existing mitigation strategies as schools and
daycares closings. Among all these variables, the viral load plays a
fundamental role. The viral load might help to predict disease severity [2]
and mortality [3, 4, 5] and can serve as a proxy for the infectivity of the
patient [6, 7, 8]. The severity of a disease, its infectivity, and its
mortality are certainly fundamental parameters that must be considered when
deciding on best-practice preventative measures to fight the pandemic spread.
Research in this direction can enable truly data-driven policymaking like, for
example, school openings and focused lockdowns. For this scope, it is
important to understand how the viral load depends on the patients’ age.

In this work, we examine viral load as a proxy for infectivity. We reanalyze
the age-stratified viral load data from Jones et al. [9] in order to better
understand the actual relation between these variables. We do this in the hope
of gaining insight into fundamental differences reported in the literature
regarding the relationship between viral load and age. We achieve this goal by
developing a flexible, non-parametric, causal, and Bayesian model to
reconstruct the conditional probability density function (PDF) of the viral
load given the patient’s age. The developed method is a second central result
of this work: it can be applied to future studies on SARS-CoV-2, to similar
data from other diseases, and also to many causally connected quantities in
very different contexts.

The non-parametric PDF reconstruction is regularized by mild assumptions on
the smoothness of the underlying statistical processes. In particular, we
assume that the log-densities are Gaussian-process realizations, drawn with an
a priori flexible correlation kernel parametrized by a Matérn family
correlation function. The parameters of this correlation function are then
inferred along with the PDF through a variational inference algorithm. To
achieve this, we adapt methods developed for information field theory, the
information theory for fields [10, 11]. In this context, fields are understood
as spatially varying (physical) quantities. The reconstructed PDFs are
regarded as scalar fields whose values are defined at each point of the two-
dimensional space spanned by age and viral load and represent the probability
of observing a given combination of age and viral load. The toolkit of
information field theory has proven itself to be successful in a wide range of
applications, ranging from 3D tomography [12], over time-resolved astronomical
imaging [13], to causality inference [14].

#### Outline

The rest of this work is structured as follows: in Sec. Related work, we
discuss the state of the art of the research on the assessment of the viral
load-vs-age dependence. In Sec. Model design, we motivate the need for a
causal description and show how this description is built into the model and
the inference scheme we adopt. In Sec. Data, we describe how the data has been
acquired and processed. We then outline the main results of our study in Sec.
Results, focusing on their impact on the infectivity of SARS-CoV-2. Finally,
in Sec. Conclusions we summarize the benefits of our approach while
highlighting its potential limitations and identifying possible future work
directions.

## Related work

Since the outbreak of SARS-CoV-2, efforts have been made in order to
understand whether certain age groups are more susceptible than others. This
could either mean that people from such age groups are more likely to get
infected or that they show more severe symptoms compared to older or younger
individuals. In addition to this, patients from specific age groups could be
more infectious than others, hence more likely vehiculating the disease. To
shed light on these problems, viral load – which is a proxy for infectivity –
can be a useful tool. It is measured by reverse transcription PCR (RT-qPCR)
assays from nasopharyngeal and oropharyngeal swabs via the so-called cycle
threshold (Ct) value. The viral load is the virus concentration in the upper
respiratory tract and it is usually expressed as the number of viral RNA
copies per mL of sample or entire swab specimen or simply by the Ct value.

Several works have analyzed whether viral loads differ between children and
adults [15, 16, 1, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27], and between
young children and adolescents [28, 29, 30], and have led to conflicting
results. At least eight studies from different countries have concluded that
SARS-CoV-2 viral RNA loads among children and adults were comparable [15, 16,
1, 19, 21, 22, 23, 24]. In these studies, dependence between age and viral
load has been tested using one and two-way analysis of variance (ANOVA) [15]
and median of the viral load. A further study [17] found that mean and median
viral load values did not vary conspicuously by age but noted that the highest
values were measured in patients born from $1995$ to $2009$. In contrast, at
least five studies reported significant differences in viral loads of young
children and adults [18, 19, 20, 25, 26, 27]. Euser et al. [18] suggested an
approximately $16$-fold higher viral load in the oldest age group ($>79$
years), compared to the youngest age group ($<12$ years). Here, age-group
differences in the viral load distribution are assessed making use of the
Kruskal-Wallis test and linear regression. Another work [20] estimated the
amount of SARS-CoV-2 in the upper respiratory tract of young children ($<5$
years) to be $10$-fold to $100$-fold greater than in adults whereas a work
from Zachariah et al. [30] showed that mean viral load was significantly
higher in infants ($<1$ year) as compared to older children and adolescents.

Finally, one of the largest and most widely followed studies on the subject –
even though the number of children and adolescents included is fairly small -
was carried out by Jones et al. [9] in early 2020. Dependence between viral
load and age has been tested for different age groups both as categorical data
and treating age as a continuous variable. In order to compare the viral load
of different age categories, the categorical data has been analyzed in a
parametric (Welch’s T-test), non-parametric (Mann-Whitney rank test), and
Bayesian fashion (modeling viral loads as a mixture of gamma distributions).
When considering age as a continuous variable, viral loads have been predicted
from age, type of PCR system, and age-PCR system interaction. This study from
Jones et al. [9] did not reveal large differences in the viral loads of
different age groups, a result that was publicly debated in Germany for its
possible implications for school opening policies. In the following, we
reanalyze this data. In a very recent publication, Jones et al. [31] extended
their initial version of the study. In this newer version, they make use of
thin-plate spline regression to conclude that children and adolescents have a
slightly lower viral load than adults, but that this difference is unlikely to
be clinically relevant.

None of the studies in the literature utilizes the causal framework. In this
work, we choose to adopt this framework since it naturally allows to answer
the central question of whether the age of a patient causes its viral load and
infectivity. It also leads to additional advantages that we will present in
the following sections. The problem of deducing causal directions (in
particular for the bivariate $x\rightarrow y$, $y\rightarrow x$ and $x\perp y$
case) from observational data coming from a joint distribution has been
introduced by the works of J. Pearl [32] and P. Spirtes et al. [33], further
developed by Mooij, J. et al. [34] and is now a central and non-trivial
problem in data analysis. For additional details and motivations behind causal
inference theory and techniques we refer to their works.

## Model design

In Sec. Related work we described the state of the art of statistical analyses
performed in order to investigate the age dependence of the viral load. These
analyses mostly rely on variance tests or correlation assessments. In the
following, we motivate the need for a causal model of the viral load
distribution. In fact, in order to explain the age and viral load data
collected by Jones et al. [9], a well-behaved model should incorporate basic
knowledge about the causal relation between age and viral load. It should
furthermore allow questioning whether the viral load distribution depends on
the patients’ age and quantify the strength of such dependence, if it exists.
We consequently develop a non-parametric causal model and apply it to data.

### Motivation

In order to make statements about the factors that contribute to the spread of
a disease, we need to rely on data to ground these claims. More importantly,
we need to find a model that is capable of identifying and explaining the
relationships that underlie the data. Indeed, this should be a minimal
requirement for any data-analysis task. The choice of an incorrect model can
lead to wrong or misleading conclusions. This is clearly an issue regardless
the nature of the data at hand. However, in the case of data describing the
pandemic spread the choice of an incorrect model can lead to ineffective or
even potentially harmful decisions. It is therefore of vital importance to be
able to capture - at least up to a certain level of uncertainty - the
interdependences underlying the data.

As discussed in the previous section Sec. Related work, the problem of
identifying the relation between age and viral load has been tackled making
use of many different techniques. While categorical data analysis and linear
correlation analysis struggle to pick up all possible dependences, non-
parametric approaches have the drawback of possibly being too (or too little)
complicated to account for the actual functional dependence that is inherent
to the data. For this reason we want to build a model that is flexible in the
sense that it can simply reconstruct independent densities with few degrees of
freedom while also being capable of inferring more complex distributions, when
needed to describe the data. Moreover, we would like the model to
automatically be able to adjust its complexity. This way, the choice of the
model is independent from the data analyst’s choices and the results are more
consistent and reproducible.

This is another benefit of adopting the causal framework. The concepts of
“simple” independence and “more complicated” dependence between variables are
natural to this framework and are described by the causal graph and the
structural causal model of choice [34]. Furthermore, claiming causal
dependence is much stronger than simple linear correlation, which also can
arise from a confounder or by mere coincidence. Causal dependence is
directional. This is a particularly interesting feature as it allows to test
for bias in the data-collection process when a unnatural causal direction is
detected from the data analysis, as we will discuss in Sec. Results.

### Causal structure

The analyzed dataset $d=\\{(i,x_{i},y_{i})\\}_{i=1}^{N}$ consists of indexed
pairs of age $x_{i}=\text{age}_{i}/\text{years}$ and log-viral load
$y_{i}=\log_{10}(\text{viral RNA copies}_{i}/\text{ml)}$ values for each of
the $N$ infected patients, extracted from Fig. 6 in [9]. The dataset is shown
in Fig. 1. For simplicity, we refer to $y$ as the “viral load”, where $10^{y}$
are the number of viral RNA copies per milliliter of sample or entire swab
specimen. We assume that the data points are Poisson process counts drawn from
an underlying stationary density distribution $\varrho(x,y)$. Using the causal
model shown in Fig. 2 we reconstruct the density $\varrho(x,y)$. Our final
density reconstructions are displayed in Fig. 3 and Fig. 4. To model this
density, we express it in terms of the underlying age distribution
$\varrho(x)$ of the patients times the conditional PDF $p(y|x)$ of the viral
load $y$ given the patient’s age $x$,

$\varrho(x,y)=\varrho(x)\ p(y|x).$ (1)

In Eq. 1 the causal direction $x\rightarrow y$ (age influences viral load) is
implicitly introduced. Even though $x\rightarrow y$ might appear to be the
most intuitive causal direction, since the immune system reaction depends on
age and thus age should affect the viral load, we note that selection effects
could introduce different apparent causal structures in the data. If the data
had been collected in such a way that the viral load ($y$) was the deciding
factor for whether a patient would enter the data sample, with an age ($x$)
dependent threshold, the viral load would impact the age distribution in the
sample, leading to an apparent $y\rightarrow x$ causal structure.

As a practical example, one could imagine physicians to request a PCR test for
a child only if their symptoms were worse compared to those of an adult for
which they would have requested a PCR test. Since this and similar selection
effects cannot be fully excluded for the analyzed data (see discussion in
[9]), we will calculate the Bayesian evidence for the possible causal
relations $x\rightarrow y$, $y\rightarrow x$, and $x\perp y$ (i.e. $x$ and $y$
are independent, $p(y|x)=p(y)$). As a final remark, we note that other effects
such as, e.g., delays in data reporting, could in principle also be
confounding the simple $x\rightleftarrows y$ causal direction. Li and White
[41] have for instance shown how reporting delays can have an effect on the
pandemic-spread models. Effectively, a confounding variable $z$ would
influence both $x$ and $y$, namely $x\leftarrow z\rightarrow y$. Although we
cannot completely exclude the presence of such a confounding variable, the
randomization tests that we describe in Subsec. Causal directions and bias
support our implicit assumption of absence of any leading confounding
variable.

To model the $x\rightarrow y$ causal direction, we need to model the age
distribution $\varrho(x)$ of the infected patients according to the causal
structure introduced in Eq. 1. Lacking knowledge on the exact details of the
patient selection process, we assume $\varrho(x)$ to be a log-normally
distributed random variable. The log-normal distribution is a natural choice
since an age density is by definition a strictly positive and continuous
quantity. Another natural assumption is the absence of abrupt changes, since
no sharp age-selecting processes are expected to have shaped it. We fulfill
these assumptions with the choice

$\varrho(x)=\varrho_{0}e^{f(x)},$ (2)

where ${\varrho_{0}=N/100}$ is a reference density and
${f:[0,100]\mapsto\mathbb{R}}$ a smooth function centered around zero. We
accordingly assume $f$ to be drawn from a zero centered Gaussian process with
a prior covariance $F$

$\mathcal{P}(f)=\mathcal{G}(f,F):=\frac{1}{\sqrt{2\pi
F}}\exp\left(-\frac{1}{2}f^{\dagger}F^{-1}f\right).$ (3)

The covariance

$F_{xx^{\prime}}=\expectationvalue{f(x)\,f(x^{\prime})}_{(f)}:=\int\mathcal{D}f\
\mathcal{P}(f)\,f(x)\,f(x^{\prime})$ (4)

determines the degree of smoothness of the logarithmic distribution function,
as well as the characteristic length scale and the amplitude of its
variations. We assume this correlation structure to be translation invariant
$F_{xx^{\prime}}=F(x-x^{\prime})$, since we only expect it to depend on age
differences - and not on a particular age value - and parametrize it with a
Matérn kernel. Invoking the Wiener-Khinchin theorem, we can represent such a
translation-invariant correlation function in Fourier space with a spectral
density of

$P_{f}(k)=\frac{a_{f}^{2}}{\left[1+(k/k_{f})^{2}\right]^{\nicefrac{{\gamma_{f}}}{{2}}}},$
(5)

with $a_{f}$ specifying the amplitude of the variations in $f$, $1/k_{f}$ the
characteristic length-scale above which the variations become uncorrelated,
and $\gamma_{f}$ the spectral index, which determines the smoothness of the
variations. We infer all three covariance parameters
$p_{f}\coloneqq(a_{f},k_{f},\gamma_{f})$ from the data. In order to ensure
that the model is flexible enough to fit the data, we set mildly informative
priors on the covariance parameters. We denote by $\mathcal{P}(f|p_{f})$ the
probability of a specific realization of $f$ given the Matérn kernel
parameters $p_{f}$, as described by Eqs. 3 to 5.

Next, we have to specify the distribution of the viral load given the age,
$p(y|x)$. To do so, we note that we can directly model an independent
distribution for which $p(y|x)=p(y)$ in the same way as we modeled $\rho(x)$,
i.e. by choosing $p(y|x)\propto e^{g(y)}$. Of course, $p(y|x)$ can in general
exhibit an arbitrary complicated dependence on $x$. We model any additional
complicated entanglement between the age $x$ and the viral load $y$ with a new
function $h(x,y)$. By doing so, we introduce a degeneracy between $g(y)$ and
$h(x,y)$, since they can both model $y$-only dependent structures. In
principle, the function $h(x,y)$ can in fact model any function $p(y|x)$
without the necessity of introducing $g(y)$. To solve this problem, we choose

$p(y|x)\propto\frac{e^{g(y)+h(x,y)}}{\int
e^{h(\tilde{x},y)}\differential{\tilde{x}}}.$ (6)

To ensure that only $g(y)$ models strictly $y$-dependent features and that all
possibly complicated $x-$dependent features are captured by $h(x,y)$, we have
to prevent $h(x,y)$ from modeling any strictly $y$-dependent structure that
has been already captured by $g(y)$. We do this with the denominator $\int
e^{h(\tilde{x},y)}\differential{\tilde{x}}$ in Eq. 6, that removes from
$p(y|x)$ any $y$-only structure contained in $h(x,y)$, hence eliminating the
degeneracy between $g(y)$ and $h(x,y)$.

We can verify that $h(x,y)$ indeed satisfies the desired property of not
encoding any structure that could be represented by $g(y)$ by simply
substituting $h(x,y)\mapsto h^{\prime}(x,y)\coloneqq h(x,y)+g^{\prime}(y)$ in
Eq. 6, where $g^{\prime}(y)$ depends only on $y$. This results in $h$ and
$h^{\prime}$ leading to the same conditional PDF

$\displaystyle\begin{split}p^{\prime}(y|x)&\propto\frac{e^{g(y)+h^{\prime}(x,y)}}{\int
e^{h^{\prime}(\tilde{x},y)}\differential{\tilde{x}}}=\frac{e^{g(y)+h(x,y)+\cancel{g^{\prime}(y)}}}{\cancel{e^{g^{\prime}(y)}}\int
e^{h(\tilde{x},y)}\differential{\tilde{x}}}\\\ &=\frac{e^{g(y)+h(x,y)}}{\int
e^{h(\tilde{x},y)}\differential{\tilde{x}}}\propto p(y|x)\end{split}$ (7)

and we can therefore conclude that only $g(y)$ can model strictly
$y$-dependent features.

As desired, for vanishing $h$ it still holds

$\displaystyle\evaluated{p(y|x)}_{h(x,y)=0}$ $\displaystyle\propto$
$\displaystyle e^{g(y)}\propto p(y),$ (8)

which implies independence, $y\perp x$. Thus, a non-trivial $h(x,y)$ models
the $x\rightarrow y$ causal influence while a trivial $h(x,y)\equiv 0$
represents causal independence. We now want to make sure that the more
complicated $x-$dependence modeled by $h(x,y)$ is only introduced if it is
strictly needed to explain the data. This way, we can clearly distinguish the
causal scenarios $x\rightarrow y$ or $y\rightarrow x$ from the independent
$x\perp y$ scenario. In fact, the distinction between $g$ and $h$ would be
meaningless without a prior choice that favors independence between $x$ and
$y$. Hence, we assume $g$ and $h$ to be drawn from zero-centered Gaussian
processes. In this way, without any information coming from the data, the most
likely realizations of both of these functions are the identically-zero
functions $g(y)=h(x,y)\equiv 0$. However, if the data exhibits strictly
$y$-dependent features, meaning that the marginal distribution $p(y)$ is non-
trivial, these features can only be represented by a non-zero $g(y)$.

Strictly $y$-dependent features are indeed clearly visible in the data. For
example, we notice that in Fig. 3 higher viral loads are by far more rare than
lower viral loads. In this case, the data-inferred $g(y)$ is non-zero and
shows this decreasing feature. Following the same reasoning, the most likely
distribution for $h(x,y)$ in absence of data is identically-zero everywhere.
Again, $h(x,y)$ will only be non-zero in case that the data triggers some
coupling between $x$ and $y$. Thus, the model favors independence of $x$ and
$y$ (by favoring $h\equiv 0$) and the inferred density will be entangled in
$x$ and $y$ only in the case the data exhibits this feature. This shows how
the level of complexity of the model is adapted to the data automatically,
without having to make any additional explicit model choice.

We again assume a Matérn-kernel-shaped correlation structure for the Gaussian
process $g$, with covariance parameters
${p_{g}\coloneqq(a_{g},k_{g},\gamma_{g})}$. We set the priors on $p_{g}$ as
for $p_{f}$ and learn these parameters from the data as well. Since the
typical length scales and amplitudes of the variations of $h(x,y)$ in $x$ and
$y$ directions are not in principle a priori similar, we assume the covariance
for $h$ to be shaped by a direct product of individual Matérn kernels in the
$x$ and $y$ directions. For their corresponding prior parameters,
$p_{h}=(p_{h}^{(x)},p_{h}^{(y)})$ with
$p_{h}^{(i)}=(a_{h}^{(i)},k_{h}^{(i)},\gamma_{h}^{(i)})$ and $i\in\\{x,y\\}$,
we use similar hyper-priors as before, i.e. as for $p_{f}$ and $p_{g}$,
respectively. For more details on the prior choices we refer to the Appendix
S1 Appendix. Matérn-kernel density reconstruction.. We call the ensemble of
all these kernel parameters $p\coloneqq(p_{f},p_{g},p_{h}).$ The details on
how the Gaussian process for $h$ with a Matérn kernel product covariance
structure is set up is described in the Appendix S1 Appendix. Matérn-kernel
density reconstruction., where we describe a multi-dimensional density
estimator that is agnostic to causal directions.

Lastly, we normalize the conditional PDF

$p(y|x)=\frac{e^{g(y)+h(x,y)}}{\int
e^{h(\tilde{x},y)}\differential{\tilde{x}}}\
\left(\int\frac{e^{g(\tilde{y})+h(x,\tilde{y})}}{\int
e^{h(\tilde{x},\tilde{y})}\differential{\tilde{x}}}\differential{\tilde{y}}\right)^{-1}$
(9)

such that the full model density reads

$\varrho(x,y)=\varrho_{0}e^{f(x)}\ \frac{e^{g(y)+h(x,y)}}{\int
e^{h(\tilde{x},y)}\differential{\tilde{x}}}\
\left(\int\frac{e^{g(\tilde{y})+h(x,\tilde{y})}}{\int
e^{h(\tilde{x},\tilde{y})}\differential{\tilde{x}}}\differential{\tilde{y}}\right)^{-1}.$
(10)

A schematic representation of the forward causal model is shown in Fig. 2. The
assumed causal structure ${x\rightarrow y}$ is introduced in the model by the
asymmetry between the roles of the $x$ and $y$ coordinates and the zero-
centered Gaussian process priors on $f$, $g$, and $h$. Interchanging $x$ and
$y$ leads to a model that follows the opposite causal direction ${y\rightarrow
x}$. This allows to empirically distinguish these causal directions by
calculating the model evidences for the two opposite scenarios, namely
${x\rightarrow y}$ and ${y\rightarrow x}$, as well as to test for $x\perp y$
by enforcing $h=0$.

Fig 1: The cobas (left) and the LC 480 (right) datasets. The lower thresholds
$y_{\text{min}}$ and $y^{\prime}_{\text{min}}$ with which the data has bene
filtered are shown in red and orange, respectively. The number of data counts
is color coded with a logarithmic color scheme (see colorbar). Fig 2: Graph
structure of the causal model of age $\mathbf{x}$ and viral load $\mathbf{y}$.
$\mathcal{F}$ denotes the Fourier transform operator. Starting from the top,
standard-normally distributed excitations $\xi$ are drawn from the latent
priors for $f,g,$ and $h$. They are then transformed into the non-parametric
signal $\rho(x,y)$ (Eq. 9) by taking the inverse Fourier transform of the
Matérn kernel parametrization described in Eq. 5. Finally, the signal
$\rho(x,y)$ is compared to the data $d$ through a Poissonian likelihood (Eq.
12).

The evidence already inherenty penalizes more complicated models (see [44]),
i.e. models that require a higher number of degrees of freedom. Using the ELBO
as a criterion for model selection has therefore the additional benefit of
regularizing the solutions. In other words, the evidence further helps the
data analyst to univocally pick the lowest complexity model which can best
explain the data, in addition to the already self-adjusting level of
complexity provided within the model itself.

Finally, we note that the choice of $x$ and $y$ is up to now completely
arbitrary. Therefore, the method described in this work can be used to assess
the causal direction of any continuous two-dimensional data distributions,
given the data and some loose (prior) information about the correlation
structure. the proposed method is thus very general and can be used for
different datasets and analysis schemes with respect to those described in
this paper. For example, it could be used to test the causal relation between
the concentration of greenhouse gasses in the atmosphere and the average
global temperature rise in Celsius degrees.

### Likelihood

In order to construct the likelihood $\mathcal{P}(d|\varrho(\cdot,\cdot))$, we
bin the data into a fine two dimensional grid over the $x$ and $y$ coordinates
with $90\times 128$ pixels, such that

$\displaystyle n_{ij}(d)=\sum_{m=1}^{N}\int_{i\Delta x}^{(i+1)\Delta
x}\differential{x}\int_{j\Delta y}^{(j+1)\Delta
y}\delta(x-x_{m})\delta(y-y_{m})\differential{y}$

contains the number of cases within the $(i,j)^{\text{th}}$ pixel of size
$\Delta x=1\,\text{yr}$ and $\Delta y\simeq 0.04\log_{10}(\text{viral RNA
copies}_{i})/\text{ml}$. These counts $n_{ij}$ are then compared with the
model’s expectations

$\displaystyle\begin{split}\lambda_{ij}&\coloneqq\lambda_{ij}(\varrho)=\int_{i\Delta
x}^{(i+1)\Delta x}\differential{x}\int_{j\Delta y}^{(j+1)\Delta
y}\varrho(x,y)\differential{y}\\\ &\approx\Delta x\Delta y\
\varrho\quantity(\quantity(i+\frac{1}{2})\Delta
x,\quantity(j+\frac{1}{2})\Delta y)\end{split}$ (11)

via a Poisson likelihood

$\mathcal{P}(d|\varrho)=\prod_{i,j}\frac{\lambda_{ij}^{n_{ij}}}{n_{ij}!}e^{-\lambda_{ij}}.$
(12)

### Inference

The full model involving the data $d$ as well as all the unknown quantities,
which compose the signal vector ${s\coloneqq(f,g,h,p)},$ reads

$\displaystyle\begin{split}\mathcal{P}(d,s)&=\mathcal{P}(d|s)\,\mathcal{P}(s),\quad\text{where}\\\
\mathcal{P}(d|s)&=\mathcal{P}(d|\varrho[f,g,h])\quad\text{and}\\\
\mathcal{P}(s)&=\mathcal{P}(f|p_{f})\mathcal{P}(g|p_{g})\mathcal{P}(h|p_{h})\,\mathcal{P}(p_{f})\,\mathcal{P}(p_{g})\,\mathcal{P}(p_{h}).\end{split}$
(13)

At this stage, we need to convert our causal model into an inference machine
for the signal vector $s$. We do this conversion by reformulating the model in
the language of information field theory [10, 11], transforming the
coordinates of the signal vector $s=s(\xi)$ such that the prior on the new
$\xi$ coordinates becomes an uncorrelated Gaussian
$\mathcal{P}(\xi)=\mathcal{G}(\xi,\mathbb{1})$ as described by Knollmüller and
Enßlin [35]. We then implement the resulting model using the Python package
Numerical Information Field Theory (NIFTy)[36, 37, 38] and finally use NIFTy’s
implementation of Metric Gaussian Variational Inference (MGVI) [39] to
approximate the posterior distribution in the new coordinates

$\mathcal{P}(\xi|d)=\frac{\mathcal{P}(d|\xi)\mathcal{P}(\xi)}{\mathcal{P}(d)}\approx\mathcal{G}(\xi-\overline{\xi}_{d},\Xi_{d})$

with a Gaussian which has posterior mean $\overline{\xi}_{d}$ and covariance
$\Xi_{d}$, where the $d$ suffix indicates the dataset used in the inference.
This Gaussian posterior encodes the approximate result of the inference in the
new coordinates. In order to translate this into the signal coordinates, we
have to transform $\mathcal{P}(\xi|d)$ to $\mathcal{P}(s|d)$ using the
relation $s=s(\xi)$. This relation is non-linear and the resulting PDF is
neither Gaussian nor practical to obtain analytically. In order to evaluate
moments from the posterior distributions of the desired quantities, MGVI
provides $\xi$-samples drawn from the approximate Gaussian posterior
$\xi\hookleftarrow\mathcal{G}(\xi-\overline{\xi}_{d},\Xi_{d})$. These
$\xi$-samples can be converted via the coordinate transformation $s=s(\xi)$
into the signal space, where they represent (approximate) signal posterior
samples. Making use of these posterior samples it is possible to calculate the
posterior expectation values and model uncertainties of any desired quantity
$q(s)$:

$\displaystyle\begin{split}\overline{q}&\coloneqq\expectationvalue{q(s)}_{(s|d)}\approx\frac{1}{N_{\text{s}}}\sum_{i=1}^{N_{\text{s}}}q(s(\xi_{i}))\\\
\sigma_{q}^{2}&=\expectationvalue{(q(s)-\overline{q})^{2}}_{(s|d)}.\end{split}$
(14)

Here, $\xi_{i}$ denotes the $i^{\text{th}}$ of the $N_{\text{s}}$ drawn
samples. In particular, the posterior mean of the conditional PDF $p(y|x)$ and
of any quantity which can be calculated from $p(y|x)$ can thereby be obtained,
as well as the resulting uncertainties characterized via their uncertainty
dispersion. In general, these posterior signal samples will not follow
Gaussian statistics because the transformation is typically non-linear.
Furthermore, since MGVI is a variational inference approach, the calculated
uncertainties will be slightly smaller compared to the ones given from the
accurate posterior. However given the complexity and size of the model, we
need to use an approximate inference method as MGVI. For details about this
methodology, as well as extensive performance and accuracy tests, we refer to
Knollmüller and Enßlin [39].

Fig 3: Cobas dataset analysis. Left: The cobas dataset. Middle: the
reconstructed density distribution $\varrho(x,y)$ as a function of the age
($x$) and the viral load ($y$) in a logarithmic coloring scheme. Right: The
$2$D conditional probability distribution $p(y|x)$ of the viral load (Eq. 9)
obtained by fitting the model. The data and the results of the analysis are
shown for two different data-filtering thresholds $y_{\text{{min}}}$ (top) and
$y^{\prime}_{\text{{min}}}$ (bottom). Fig 4: LC 480 dataset analysis. The same
as in Fig. 3, but for the LC $480$ dataset.

### Identifying causal directions with the Evidence Lower Bound

We now focus on understanding the causal relations given by the interplay
between the variables. The proposed causal model should allow to discriminate
between all possible causal structures, namely $x\rightarrow y$, $y\rightarrow
x,$ and $x\perp y$. Since we have a different model for each of these causal
directions, we can use the Bayesian evidence to select the model that is
better suited to the data. We do this once again by exploiting our variational
inference scheme. It has been shown that the evidence is indeed a consistent
and robust criterion for model selection [43, 42].

First, we need to define the independent model, i.e. a model for which age and
viral load are regarded as statistically independent variables. Such a model
is built by setting a very tight zero-centered prior on $h$, hence removing it
from the model for all practical purposes. We can then estimate the Bayesian
evidence both for the causal – hence dependent – models ($x\rightarrow y$ and
$y\rightarrow x$) and compare it with the evidence of the independent model
($x\perp y$). More precisely, we calculate the so called Evidence Lower Bound
(ELBO) as a proxy for an evidence. Using the posterior uncertainty covariance
$\Xi$ as well as the posterior samples provided by MGVI, we can compute the
ELBO [44] for each model. This is lower than the exact logarithm of the
evidence by the information difference (as measured in nits) between the exact
and approximated posterior of a model. If too much information is not lost
during MGVI, the ELBO should be a good approximation of the exact log
evidence. Furthermore, we can assume that deviations from the exact log
evidence should be similar among different models, thereby reducing the effect
of the MGVI approximations on differences between log-evidences. Thus, we can
use the ELBO as a good proxy for the log model evidence ratios of similar
models. The stochastic sampling steps performed in order to estimate the ELBO
introduce a sampling uncertainty. This uncertainty can be in principle reduced
by taking more samples, at the expense of larger computational costs. We state
this numerical one-sigma uncertainty for all MGVI and ELBO based log-
evidences. The $y\rightarrow x$ model is obtained by swapping the coordinates
of the $x\rightarrow y$ model. The quantity of interest is then the logarithm
of the evidence ratio of each causal model with the independent one,

$\Delta E_{x\rightleftarrows y}=\log\frac{p(d|x\rightleftarrows y)}{p(d|x\perp
y)},$ (15)

where $x\rightleftarrows y$ denotes either $x\rightarrow y$ or $y\rightarrow
x$. $\Delta E_{x\rightarrow y}$ indicates the log evidence in favor of the
causal model $x\rightarrow y$ with respect to the independent one and
similarly $\Delta E_{y\rightarrow x}$ the one for $y\rightarrow x$. Comparing
$\Delta E_{x\rightarrow y}$ with $\Delta E_{y\rightarrow x}$ also allows to
discriminate between the two possible causal directions in the dataset. We
identify the preferred causal direction underlying the data in subsec. Causal
directions and bias, where we also discuss its implications.

## Data

We make use of RT-PCR viral load data collected from the Charité Institute of
Virology and Labor in Berlin, Fig. 6 in Jones et al. [9]. The data was
acquired with two different PCR instruments, Roche cobas $6800$/$8800$ (cobas
dataset, which we denote by $d_{\text{C}}$ and is comprised of $\approx 2200$
data points) and Roche LightCycler $480$ II (LC $480$ dataset, which we denote
by $d_{\text{L}}$ comprised of $\approx 1350$ data points). In the following,
we will show the difference between the two datasets. As can be seen from the
count difference in the raw data plots as of Fig. 6 of [9] and in the
histogramed data in Figs. 3 and 4, for low viral loads ($y\lesssim 5$ in units
of $\log_{10}(\text{RNA copies/ml})$), the LC $480$ dataset shows a roughly
uniform count distribution in the whole viral load domain. In contrast, the
cobas dataset exhibits an increasing number of counts in the $y\in[2.0,3.8]$
viral load domain followed by a descending trend in counts in the
$y\in[3.8,5.0]$ region. For this reason and in order to better understand the
possible shortcomings of both instruments, we analyzed the data in two
different ways.

Since the major differences between the datasets arise for viral load values
$y\in[3.8,5.0]$ we define two lower thresholds for the viral load
($y_{\text{min}}\coloneqq 3.8$ and $y^{\prime}_{\text{min}}\coloneqq
5.4\simeq\log_{10}(250000)$ in units of $\log_{10}(\text{viral RNA
copies/ml})$) and discard any data point for which the viral load is lower
than $y_{\text{min}}$ and $y^{\prime}_{\text{min}}$, respectively. We set the
lower threshold $y_{\text{{min}}}$ at the value for which the number of counts
of the cobas dataset is maximum (see Fig. 1). Below this threshold, the
counts’ density lowers dramatically. As for $y^{\prime}_{\text{{min}}}$, the
value of $250000$ in units of viral RNA copies/ml indicates the threshold for
the isolation of infectious virus in cell cultures at more than $5\%$
probability as described by Wörfel et al. [40]. We then analyze the cobas and
LC $480$ datasets first neglecting the data below $y_{\text{min}}$ threshold
and then below $y^{\prime}_{\text{min}}$. For the sake of simplicity, we will
denote the cobas and LC $480$ datasets with $y_{\text{min}}$ and
$y^{\prime}_{\text{min}}$ as a lower threshold as $d_{\text{C}}$,
$d^{\prime}_{\text{C}}$, $d_{\text{L}}$, and $d^{\prime}_{\text{L}}$,
respectively. This way we can highlight the differences between the two
datasets and investigate possible sources of systematic errors in the viral
load measurement instruments or different selection effects in the data
collection process.

The datasets are not explicitly provided by the authors. Therefore, we
acquired the data by means of a plot digitizer algorithm from Fig. 6 of [9].
Since the age coordinates are not labeled precisely, but only a rough interval
$\Delta_{\text{age}}\sim 0-100\,\text{yr}$ is provided, we do not expect the
acquired data points to be accurate – especially in the age domain. For this
reason the obtained age axis could be affected by a global shift up to $5-10$
years in any direction. But, since our model is translation invariant, this
kind of systematic bias in the data extraction process (an overall shift in
the age values for all data points) does not affect the causal inference
machinery. Nonetheless, the reader should keep this in mind when interpreting
the results concerning the conditional PDF of the viral load given the age.
Concerning the viral load coordinates, for which more precise units were
given, the data should be regarded as more reliable. For our aim of building a
method to analyze age and viral load data in a non-parametric and causal
fashion providing uncertainty estimates, this level of accuracy is sufficient.
The results we present from here on are given for the measured coordinate
values without considering any uncertainty with respect to the real quantities
(age, viral load). Nevertheless, we do not believe systematic or random error
contributions in the data extraction procedure to significantly affect our
results since the shapes of the learned distributions are translation
invariant.

## Results

In the following, we discuss the main results and their relevance for the
infectivity of different age groups.

### Age dependence of the viral load

For both lower thresholds and datasets $d_{\text{C}}$, $d_{\text{L}}$, we
estimate the model parameters $s$ by means of the MGVI algorithm implemented
in NIFTy. In Figs. 3 and 4 we show the data (left panel) and the correspondent
reconstructed underlying densities $\varrho(x,y)$ (central panel).

A multi-modal age distribution is clearly visible as well as an overall
decrease for growing viral loads in all of the reconstructed densities. The
conditional PDF $p(y|x)$ of the viral load $y$ given the age $x$ is shown in
Fig. 5. In the conditional probability, the multi-modal age structure
displayed by the density is not visible since it has been absorbed by
$\varrho(x)$. What this effectively means is that age-only selection effects –
e.g., testing one or many specific age groups more than others or demographics
in general – have been modeled, and the resulting conditional probability
distribution does not depend on such effects.

For all datasets and ages, a general descending trend in the viral load
probability distribution is clearly visible. The reconstruction based on the
$d_{\text{C}}$ dataset exhibits significant differences in the viral load for
different ages (Fig. 5, first panel). For infected patients approximately
above the age of $60$, the distribution exhibits a distinct maximum for viral
loads of $y\simeq 8$ (in units of $\log_{10}(\text{RNA copies/ml})$).
Furthermore, we show that this feature is indeed triggered by the data, and is
not just the result of an over-fit of sample noise. To do so, we apply a
random permutation $r$ to the viral load values in the $d_{\text{C}}$ dataset,
the only one that exhibits a possible ${x\rightarrow y}$ causal structure. We
then analyze the randomized dataset
${d_{\text{C}}^{r}=\\{(i,x_{i},y_{r(i)})\\}_{i=1}^{N}}$ in the same way as
seen for ${d_{\text{C}}=\\{(i,x_{i},y_{i})\\}_{i=1}^{N}}$.The resulting
conditional PDF $p^{r}(y|x)$ reconstructed from the randomized dataset (Fig.
5) does not exhibit any clear age-dependent structure in the viral load,
indicating that the differences seen in the real data are not just a shot
noise effect.

Fig 5: Viral-load probabilities for specific ages. Viral load conditional PDF
$p(y|x)$ for specific ages $x\in\\{10,50,80\\}$. Panel one and two (top)
display the results for the cobas dataset $d_{\text{C}}$ and for the
randomized cobas dataset $d_{\text{C}}$ respectively. The latter serves as a
null test, since the randomization erases any (causal) relation between the
age ($x$) and the viral load ($y$) other than shot noise. Panel three and four
(bottom) show $p(y|x)$ for the cobas dataset $d^{\prime}_{\text{C}}$ and for
the LC $480$ dataset $d^{\prime}_{\text{L}}$ respectively, both with the
higher viral load threshold $y^{\prime}_{\text{{min}}}$. The shaded regions
represent $1\sigma$ and $2\sigma$ uncertainty contours of the approximate
posterior.

### Causal directions and bias

As discussed in Subsec. Causal structure, we can use ELBO ratios to identify
causal relations. In our analysis, this corresponds to choosing between the
$x\rightarrow y$, $y\rightarrow x,$ and $x\perp y$ models. Identifying the
causal direction underlying the data will allow to show that selection effects
distorting the expected causal relation $x\rightarrow y$ are subdominant and
that for the cobas dataset $d_{\text{C}}$ we indeed see evidence for an age
dependence of the viral load distribution, $p(y|x)\neq p(y)$.

As a matter of fact, the log-evidence ratio for the $d_{\text{C}}$ dataset is
$\Delta E_{\text{C},x\rightarrow y}=4.6\pm 1.0$ for $d_{\text{C}}$, which
clearly favors the dependent model, but this value decreases to $\Delta
E^{\prime}_{\text{C},x\rightarrow y}=-1.5\pm 1.0$ when considering
$d^{\prime}_{\text{C}}$, hence $y^{\prime}_{\text{min}}$ as a lower threshold.
We highlight that a log-evidence difference between the two compared models of
$1$ unit corresponds to a factor of $e\approx 2.7$ for the Bayesian odds ratio
between the two. Thus, $\Delta E_{\text{C},x\rightarrow y}=4.6\pm 1.0$ implies
given equal model priors, $p(x\rightarrow y)=p(x\perp y)$, a posterior model
odds ratio of $p(x\rightarrow y|d):p(x\perp y|d)=e^{3.5\pm 0.7}\approx
99.5_{[36.6]}^{[270.4]}$ in favor of a causal dependence between viral load
and age for the cobas dataset with the lower threshold $y_{\text{min}}$.

For the opposite causal direction we get $\Delta E_{\text{C},y\rightarrow
x}=-0.2\pm 1.1$, which shows that there is no strong $y\rightarrow x$
structure in the data. For the LC dataset $d_{\text{L}}$ these evidence
differences with respect to the independent model become $\Delta
E_{\text{L},x\rightarrow y}=-4.6\pm 1.0$ and $\Delta
E^{\prime}_{\text{L},x\rightarrow y}=-3.4\pm 1.0$ respectively for the two
thresholds. Therefore the independent model is favored in both cases. This
shows that the cobas ($d^{\prime}_{\text{C}}$) and LC $480$
($d^{\prime}_{\text{L}}$) datasets are in agreement for viral loads which are
higher than $y^{\prime}_{\text{min}}$, but in case we include the data lying
in the viral load region $y\in[y_{\text{{min}}},y^{\prime}_{\text{{min}}}]$,
the causal age-dependent structure becomes visible in $d_{\text{C}}$ and is
not negligible anymore.

It is known that model evidences can vary strongly with different data
realizations. Moreover, in order to calculate the evidences, we invoked
approximations and stochastic calculation steps. Thus, proper null-tests are
required in order to validate and calibrate the evidence ratio calculation. We
provide such null-tests by repeating the data randomization step described
above several times, thereby producing many randomized datasets. By
construction, these dataset should not exhibit any causal structure. Indeed,
for the $10$ randomized-dataset realizations performed, we find much lower
log-evidence ratios between dependent (causal) and independent models with
respect to the ones found for the original cobas dataset $d_{\text{C}}$. The
average difference between these tests is $\expectationvalue{\Delta
E_{\text{random},x\rightarrow y}}_{\text{randomizations}}=-6.0\pm 1.0$. Since
none of the randomized dataset realizations reaches comparably high log-
evidence ratios with respect to the independent model, all these findings
support robustly the argument that the dependent model is a more suited
description of the $d_{\text{C}}$ dataset. As previously discussed, this
result also supports the assumption that no additional variable $z$ is
confounding $x\leftarrow z\rightarrow y$ the bivariate causal direction
$x\leftrightarrows y$. Had there been a strong confounder $z$, we would have
expected a more complicated data distribution that could have possibly been
detected when comparing the actual data with the randomized datasets. Even if
this did not appear to be the case, we cannot completely rule out the presence
of a weak confounder.

The results of the evidence calculations are displayed in Fig. 6. This figure
also indicates that the independent model interpretation of the data is
favored for all other datasets and threshold combinations (except for
$d_{\text{C}}$), since the evidence for the independent model is always
higher.

These contradicting results for the two datasets (or thresholds)
$d_{\text{C}}$ and $d^{\prime}_{\text{C}}$ might have several possible
explanations. First, they could be caused by a potential accuracy loss of the
PCR devices below certain viral load values, as suggested for the cobas
dataset $d_{\text{C}}$ below $y^{\prime}_{\text{{min}}}$ in [9], hence for the
Roche cobas 6600/8800 PCR system. It could also mean that the opposite is
happening and the Roche LC $480$ PCR device is less sensitive than the Roche
cobas 6600/8800 in the $y\in[3.8,5.0]$ region. This possibility would be
supported by the fact that the cobas dataset exhibits a clear age dependent
structure in such viral load region, but the (swab) data processed by the
cobas PCR device contains no information on the patients’ ages. Hence it would
be surprising that a systematic effect in the measurements could introduce an
age dependence on the viral load distribution. Furthermore, this pattern –
that older patients exhibit higher viral loads – is plausible from a medical
perspective.

Nevertheless, we cannot exclude the possibility that selection effects have
been introduced in the data. We have already shown that “viral load causing
age” effects ($y\rightarrow x$) are subdominant. Nevertheless, selection
effects could still have been introduced for instance by collecting age and
viral load subsamples from a “viral-load biased” population sample. This could
happen, e.g., if symptomatic patients had been predominantly tested. Since
children are less likely to show symptoms than adults, the sample would then
include mainly those children with higher viral loads.

And finally, a combination of such competing effects could have affected the
results and thereby imprinted a spurious age dependence to the viral-load
distribution. Given only the statistical data in our possession, this
possibility cannot be ruled out completely.

Fig 6: Evidence results. Natural logarithm of the evidence ratio with respect
to the correspondent independent model for the $x\rightarrow y$ cobas datasets
$d_{\text{C}}$ and $d^{\prime}_{\text{C}}$, for the $y\rightarrow x$ causal
model $d_{\text{C}}^{\text{inv}}$ and
$\left.d^{\prime}\right._{\text{C}}^{\text{inv}}$(left panel). The log-
evidence ratios labeled with $d^{r}_{C}$ display the ($y$-)randomized datasets
(for 10 different realizations). The same is also shown for the LC $480$
dataset (right). The error bars represent the numerical uncertainty associated
to the stochastic estimate of the ELBO. A positive logarithm of the evidence
ratio denotes a preference for the given causal model with respect to the
correspondent independent model and vice versa. The strength of this
preference can be determined by taking the exponential of the log-evidence
ratio between the causal model and the correspondent independent one. For
example, a log-evidence ratio of $5$ corresponds to a $e^{5}\approx 141$-fold
preference for the examined causal model.  Fig 7: Infectivity proxy vs. age.
$I(y)$ obtained from Fig. 1g in [40] and fit with a probit function (top left,
Original), then projected for different patients’ ages for the cobas dataset
$d_{\text{C}}$ (top right), $d^{\prime}_{\text{C}}$ (bottom left) and for the
LC $480$ dataset $d_{\text{L}}$ (bottom right). The upper and lower curves are
obtained by translating $I(y)$ of one unit in viral load as in $I(y-1)$ and
$I(y+1)$ in order to give upper and lower bounds similar as those shown in
Fig. 1g in [40]. The shaded regions show $1\sigma$ and $2\sigma$ uncertainty
contours of the approximate posterior.

### Impact on infectivity

After having established a potential age difference in the viral load
distribution for $d_{\text{C}}$, we investigate whether this difference \- if
real - would be relevant for the infection dynamics. For this purpose, we need
to link the viral load to the infectivity $I(y)$ of the virus, i.e. the
probability of transmitting the infection. Infectivity can be measured in
different ways. In our analysis, we choose the projected virus isolation
success based on probit distribution described in [40] as a proxy for
infectivity. This represents the infection success rate for cell cultures
exposed to saliva with different viral load $y$ and can be seen as the blue
curve in Fig. 7, labeled with “Original”. We will from now on refer to this
parameter as “infectivity proxy” and indicate it with $I(y)$.

The projected infectivity as a function of the age is then given by the
expectation value of the infectivity parameter over the conditional PDF
$p(y|x)$

$I(x)\coloneqq\expectationvalue{I(y)}_{p(y|x)}=\int
I(\tilde{y})p(\tilde{y}|x)\differential{\tilde{y}}.$ (16)

The result, together with the uncertainties resulting from our PDF modeling
are shown in Fig. 7. No relative differences larger than $0.18$ for the
(projected) infectivity of the different age groups is found, with typical
values of $I(x)\approx 0.3$ for all datasets. This means that at most a $50\%$
difference in infectivity due to different viral load between different age
groups \- but more likely a smaller one - should be expected.

In order to characterize the uncertainty resulting from our $I(y)$ modeling,
due to the uncertainty of the original determination of this function and due
to the uncertainty in the identification of viral loads with different
instruments, we repeat the analysis while shifting the original $I(y)$ curve
by one order of magnitude upwards and downwards in $y$. The resulting maximal
relative difference in the infectivity of the different age groups is $\simeq
0.3$. Thus, even though our model allows us to show that infectivity exhibits
an age dependence if the $d_{\text{C}}$ dataset with $y_{\text{{min}}}$
provides a valid picture, the viral load differences between different age
groups though are not strong enough to impact on the infection dynamics at a
level that justifies regarding any age group as noninfectious or even
significantly less infectious.

## Conclusions

In order to investigate the controversial results reported in the literature,
we developed a causal model to assess the dependence of viral loads of
patients infected with COVID-19 on age. The developed model is capable of
reconstructing two-dimensional density distributions from data counts and to
learn causal directions. The model complexity is set in a user-independent
fashion making the results more robust and consistent. Furthermore, its causal
nature allows to make predictions about non-directly measured quantities (e.g.
the infectivity assessment in Sec. Results and Fig. 7) and to additionally
test for bias in the data-collection process.

Although the benefits of a causal analysis were already discussed in Pearl’s
original work [32] and are usually well recognized, causal inference is not so
commonly applied to real-world data science, often because it requires the
implementation of complicated ad-hoc models. With the hope of making it useful
for a wide range of data-science applications we make our causal model freely
available. This model is very flexible and generic (as it only requires a set
of $x$ and $y$ data pairs and mild priors on their correlation structures) and
it can therefore be used in future epidemiological studies as well as in
completely different fields. We provide the source code under an open source
license for usage in further studies and applications at
https://gitlab.mpcdf.mpg.de/ift/public/causal_age_viral_load_model. As a side
product, we also developed a causal-direction-agnostic density estimator,
which is described in more detail in the S1 Appendix. Matérn-kernel density
reconstruction..

Using our novel method to model causal relations non-parametrically, we have
reanalyzed the SARS-CoV-2 age and viral load data presented in [9].In doing
so, we have found statistically significant differences in the viral load
distribution of different age groups when regarding the cobas dataset
$d_{\text{C}}$ for viral loads within the interval of $10^{3.8}$ to $10^{5.0}$
in units of viral RNA copies/ml of sample or entire swab specimen as reliable.
These differences become irrelevant if this region is ignored in the analysis.

We cannot completely exclude that selection effects in the data-collection
process may have introduced an apparent causal relation between viral load and
age, but the observed trend – a statistically-significant increase in the
viral load with age – fits with the generally accepted notion that the immune
system response gets weaker with age. Assuming this trend to be real, we
showed, however, that its expected impact on the infectivity of different age
groups is at most moderate. For this reason we cannot exclude any age group
from being considered as a potentially significant source of infection.

The region of the cobas dataset relevant for this trend is described in [9] as
containing an artifact, suggesting that the correct interpretation of the data
is that viral load, hence infectivity, is predominantly age independent. Here,
we want to point out that other studies on the age dependence of the viral
load present in the literature [18, 15] make the opposite claim. Moreover, in
their most recent publication, Jones et al. [31] acknowledge an age dependence
of the viral load. This dependence is quantitatively similar to the one we
have detected with our method. Furthermore, the causal evidence tests
presented in Sec. Results favor considering the age dependence of the viral
load as a real effect and not just as an artifact. These tests also disfavor
the reverse-causal-direction model (here: that the viral load of a patient
“causes” its age), which would indicate that strong selection effects have
affected the data-collection process. We introduced these tests as a new tool
to detect potential systematic effects in similar datasets.

In conclusion, the results of our analysis ultimately confirm some of the
findings in the literature - i.e. that the viral load is only modestly
dependent on the age - but with a much higher sensitivity and robustness. Of
central importance are the methods here developed. While being tailored to
describe the Covid-19 pandemic data, they can be easily adapted for more
general purposes and can prove very useful also for future pandemics or for
new – and possibly more infective – mutations of SARS-CoV-2.

## Supporting information

#### S1 Appendix. Matérn-kernel density reconstruction.

Matérn-kernel density reconstruction. We aim to infer an unknown distribution
from a random realization of discrete data. In the following, we present a
general-purpose density estimator that serves this scope. Using a fully
Bayesian framework, we can propagate uncertainties through each inference step
and extract information from the correlation structure inherent to the data.
We provide our method as open-source software.

According to Chen, Tareen, and Kinney [45], the problem of extracting a smooth
density function from a limited set of data samples is a challenging and well-
known problem in statistical learning and data analysis. The most common ad-
hoc methods to empirically derive (probability) densities from data usually
involve histograms or Kernel Density Estimation (KDE) [46, 47]. These methods
do not infer the smoothness of the learned (probability) density’s correlation
structure and are thus prone to reconstructing unphysical densities. Other
methods make use of neural networks (see for example Liu et al. [48]) or
restrict the density to specific functional forms (see e.g. Dirichlet Process
Mixture Model, or DPMM [49, 50, 51]). Another approach is to use smooth priors
and infer the level of smoothness of the reconstructed density from data via
maximum entropy (Density Estimation using Field Theory, or DEFT [52, 53]).
Chen, Tareen, and Kinney [45] propose an interesting information-theoretical-
based modification of DEFT, although it only effectively works in one
dimension. Finally, another very commonly used and effective solution to the
density estimation problem is the one given by deep neural networks. In a
recent work, Liu et al. [48] propose a generative adversarial networks (GAN)
which is particularly effective in high dimensions. We refer to their work for
comparison with similar neural-network-based approaches.

Most of these approaches lack a robust estimate for uncertainties or specify
none at all. In the hope of addressing these shortcomings, we propose our
novel general-purpose density reconstruction method, which we will refer to as
Matérn Kernel Density Estimator (MKDE). This method is very general, works for
a generic $n$-dimensional space, and therefore applies to many different
contexts and fields. We presented one paradigmatic example of its many
possible applications in Sec. Causal structure. In this example, we used MKDE
to reconstruct the continuous distributions of the ages and viral loads of
patients infected with Covid-19 from age-and-viral-load data samples collected
within the general population.

MKDE is capable of reconstructing a smooth density distribution underlying an
– even limited – discrete dataset. We achieve this result under the hypotheses
that the data points are drawn from the underlying density through a Poisson
process and that the reconstructed density is a sufficiently smooth function.
Since we expect the inferred (probability) density to be strictly positive and
to vary on logarithmic scales, we choose the log-normal model

$\varrho(x)=e^{s(x)},$ (17)

where $s(x)$ is the natural logarithm of the density. In the case of multi-
dimensional data, both the density $\varrho(x)$ and its logarithm $s(x)$ in
Eq. 17 are functions defined over $n$-dimensional vector spaces, with $n$
being the dimension of the data space (e.g. space, time, age, viral load, …).
For the sake of simplicity, we will initially show how a one-dimensional
density is reconstructed, emphasizing how to generalize to the multi-
dimensional case only when this generalization is non-trivial.

Fig 8:  Example of MKDE performance. Left: random realization of a smooth
probability density with independent covariance structure in the $x$ and $y$
direction (ground truth, GT). Middle: Poissonian counts drawn from the GT.
Right: MKDE reconstruction from the counts. In the same panel, we also show
the marginals of the reconstructed density (blue lines), together with their
one and two-sigma uncertainty estimates (shaded areas in blue) against the
marginals of the GT (orange lines).

Smoothness is a common and ubiquitous assumption when dealing with physical
data. To fulfill this assumption in the MKDE, we parametrize the two-point
correlation structure of the Gaussian process that determines the value of the
log-density $s(x)$ at each point and for each dimension with a Matérn kernel.
This kernel is a very flexible choice [54]. Moreover, it is very well suited
to represent a priori homogeneous covariance structures like the ones we want
to model. We define the Gaussian process
$s:=s(x)=(s_{x})_{x\in\mathbb{R}}\hookleftarrow\mathcal{P}(s)=\mathcal{G}(s,S)$
on the one-dimensional position space $\mathbb{R}$. This process has a
homogeneous covariance structure $S_{xx^{\prime}}=C_{s}(x-x^{\prime})$ that
can be efficiently represented in Fourier space thanks to the Wiener-Khinchin
theorem. Invoking this theorem, we can use the power spectrum $P_{s}(k)$ to
fully determine the Fourier transform of the two-point correlation function
$C_{s}(x-x^{\prime})$ for any stationary and statistically homogeneous
process, and in particular for the Gaussian process $s$.

In order to draw prior samples from the Gaussian field $s$, we choose a
standardized coordinate system $\hat{\xi}$. We then transform the standard
normally distributed parameters $\xi=(\xi_{k})_{k}$, where $k\in\mathbb{N}$ is
the Fourier-space index according to the mapping

$s=\mathcal{F}^{-1}\,A\,\xi\quad\text{with
}A\coloneqq\text{diag}(\sqrt{P_{s}}).$ (18)

Here, $\mathcal{F}$ represents the Fourier transform operator and $A$ the
amplitude operator in Fourier space. The amplitude operator encodes the
Matérn-kernel correlation structure, parametrized with

$A_{kk^{\prime}}=\,2\pi\,\delta(k-k^{\prime})\frac{a_{s}}{\left[1+k^{2}/k_{0}^{2}\right]^{-\gamma_{s}/4}},$
(19)

where $a_{s}$ is a scale factor which accounts for the standard deviation in
position space, $k_{0}$ is the magnitude of the characteristic correlation-
length wavevector, and $\gamma_{s}$ is the spectral index of the power
spectrum. We assume $a_{s}$ and $k_{0}$ to be a priori log-normally
distributed since we expect strictly positive variations of the possible power
spectra on a logarithmic scale. Similarly, we choose $\gamma_{s}$ to be
normally distributed, since the spectral index could in principle also be
negative. We additionally introduce volume factors to ensure that the model
parameters are intensive with respect to volume, i.e. they do not depend on
the volume in position space.

For higher-dimensional data, we expect the correlation structure along each
axis (or dimension) to be a priori independent from the others. These
different axes could in fact have very different meanings (and units), as they
could represent – for instance – space and time, temperature, pressure, and
volume, age and viral load (as seen in Sec. Causal structure), or a different
combination of these and other continuous quantities. Therefore, in the
$n$-dimensional data space we can decompose the amplitudes of the correlation
structure

$A_{k,q,\ldots,n}=\bigotimes_{i\in\\{k,q,\ldots,n\\}}A_{i}$ (20)

along each independent axis, each modeled by an individual amplitude operator
$A_{i}$, for $i\in\\{k,q,\ldots,n\\}$. The zero modes of the individual axes
must be treated separately in order to avoid degeneracy. Thus, in the proposed
model the zero mode is shared among all directions and inferred independently
through an a priori strictly-positive and uniformly-distributed parameter
$\alpha$. For more details on the zero mode degeneracy and factorizing power
spectra, we refer to Arras et al. [55].

At this stage, we can summarize all the parameters that we have introduced for
each independent axis with the scalar-valued parameters $\alpha,a_{s},k_{0}$,
and $\gamma_{s}$ and the vector-valued $\xi_{k}$. We set broad priors on these
parameters and learn them using MGVI. We set the following priors on the
signal parameters: $\alpha=[10^{-15},5.0]$, $a_{s}=(0.3\pm 0.2)$,
$k_{0}=(4.0\pm 3.0)$, and $\gamma_{s}=(-6.0\pm 3.0)$, where the mean and
standard deviation specify a Gaussian prior distribution for $\gamma_{s}$ and
log-normal distributions with the given mean and standard deviation for
$a_{s}$ and $k_{0}$. For details on the inference of posterior estimates for
the MKDE parameters through MGVI and their uncertainty quantification, we
refer to Sec. Inference. Fig. 8 illustrates the performance of MKDE in a two
dimensional setting.

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

#### S2 Appendix.

Priors. Throughout the analysis, we have made the use of the following priors
on the signal parameters: ${a_{f}=(0.3\pm 0.1)}$, ${k_{f}=(5.0\pm 3.0)}{\
\text{yr}^{-1}}$, and ${\gamma_{f}=(-3.0\pm 2.12)}$, where the mean and
standard deviation specify a Gaussian prior distribution for $\gamma_{f}$ and
log-normal distributions with the given mean and standard deviation for
$a_{f}$ and $k_{f}$.

## Acknowledgments

We thank Dr. Matthia Sabatelli (Montefiore Institute, University of Liège) for
feedback on preliminary versions of this work and the Information Field Theory
group at the Max Planck Institute for Astrophysics for the fruitful
discussions.

## References

  * 1. Colson P, Tissot-Dupont H, Morand A, Boschi C, Ninove L, Esteves-Vieira V, et al.  Children account for a small proportion of diagnoses of SARS-CoV-2 infection and do not exhibit greater viral loads than adults.  European Journal of Clinical Microbiology & Infectious Diseases. 2020;39(10):1983–1987.
  * 2. Zheng S, Fan J, Yu F, Feng B, Lou B, Zou Q, et al.  Viral load dynamics and disease severity in patients infected with SARS-CoV-2 in Zhejiang province, China, January-March 2020: retrospective cohort study.  Bmj. 2020;369.
  * 3. Pujadas E, Chaudhry F, McBride R, Richter F, Zhao S, Wajnberg A, et al.  SARS-CoV-2 viral load predicts COVID-19 mortality.  The Lancet Respiratory Medicine. 2020;8(9):831–934.  doi:10.1016/S2213-2600(20)30354-4.
  * 4. Fajnzylber J, Regan J, Coxen K, Corry H, Wong C, Rosenthal A, et al.  SARS-CoV-2 viral load is associated with increased disease severity and mortality.  Nature Communications. 2020;11(1):5493.  doi:10.1038/s41467-020-19057-5.
  * 5. Maltezou HC, Raftopoulos V, Vorou R, Papadima K, Mellou K, Spanakis N, et al.  Association Between Upper Respiratory Tract Viral Load, Comorbidities, Disease Severity, and Outcome of Patients With SARS-CoV-2 Infection.  The Journal of Infectious Diseases. 2021;doi:10.1093/infdis/jiaa804.
  * 6. Yonker LM, Neilan AM, Bartsch Y, Patel AB, Regan J, Arya P, et al.  Pediatric severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2): clinical presentation, infectivity, and immune responses.  The Journal of pediatrics. 2020;227:45–52.
  * 7. Jefferson T, Spencer E, Brassey J, Heneghan C.  Viral cultures for COVID-19 infectivity assessment. Systematic review.  medRxiv. 2020;.
  * 8. Kawasuji H, Takegoshi Y, Kaneda M, Ueno A, Miyajima Y, Kawago K, et al.  Transmissibility of COVID-19 depends on the viral load around onset in adult and symptomatic patients.  PLoS ONE. 2020;15(12).  doi:10.1371/journal.pone.0243597.
  * 9. Jones TC, Mühlemann B, Veith T, Biele G, Zuchowski M, Hoffmann J, et al.  An analysis of SARS-CoV-2 viral load by patient age.  https://virologie-ccm.charite.de/fileadmin/user_upload/microsites/m_cc05/virologie-ccm/dateien_upload/Weitere_Dateien/Charite_SARS-CoV-2_viral_load_2020-06-02.pdf
  * 10. Enßlin TA, Frommert M, Kitaura FS.  Information field theory for cosmological perturbation reconstruction and nonlinear signal analysis.  Physical Review D. 2009;80(10):105005.
  * 11. Enßlin TA.  Information theory for fields.  Annalen der Physik. 2018; p. 1800127.
  * 12. Leike RH, Glatzle M, Enßlin TA.  Resolving nearby dust clouds.  Astronomy & Astrophysics. 2020;639:A138.
  * 13. Arras P, Frank P, Haim P, Knollmüller J, Leike R, Reinecke M, et al.  The variable shadow of M87.  arXiv preprint arXiv:200205218. 2020;.
  * 14. Kurthen M, Enßlin T.  A Bayesian Model for Bivariate Causal Inference.  Entropy. 2020;22(1):46.
  * 15. Mahallawi WH, Alsamiri AD, Dabbour AF, Alsaeedi H, Al-Zalabani AH.  Association of Viral Load in SARS-CoV-2 Patients With Age and Gender.  Frontiers in Medicine. 2021;8:39.  doi:10.3389/fmed.2021.608215.
  * 16. Jacot D, Greub G, Jaton K, Opota O.  Viral load of SARS-CoV-2 across patients and compared to other respiratory viruses.  Microbes and infection. 2020;22(10):617–621.
  * 17. Kleiboeker S, Cowden S, Grantham J, Nutt J, Tyler A, Berg A, et al.  SARS-CoV-2 viral load assessment in respiratory samples.  Journal of Clinical Virology. 2020;129:104439.
  * 18. Euser S, Aronson S, Manders I, van Lelyveld S, Herpers B, Jansen R, et al.  SARS-CoV-2 viral load distribution in different patient populations and age groups reveals that viral loads increase with age.  medRxiv. 2021;doi:10.1101/2021.01.15.21249691.
  * 19. Costa R, Bueno F, Albert E, Torres I, Carbonell-Sahuquillo S, Barres-Fernandez A, et al.  Upper respiratory tract SARS-CoV-2 RNA loads in symptomatic and asymptomatic children and adults.  medRxiv. 2021;.
  * 20. Heald-Sargent T, Muller WJ, Zheng X, Rippe J, Patel AB, Kociolek LK.  Age-related differences in nasopharyngeal severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) levels in patients with mild to moderate coronavirus disease 2019 (COVID-19).  JAMA pediatrics. 2020;174(9):902–903.
  * 21. Chung, E., Chow, E., Wilcox, N., Burstein, R., Brandstetter, E., Han, P., Fay, K., Pfau, B., Adler, A., Lacombe, K. & Others  Comparison of symptoms and RNA levels in children and adults with SARS-CoV-2 infection in the community setting.  JAMA Pediatrics. 175, e212025-e212025 (2021)
  * 22. Baggio, S., Arnaud, G., Yerly, S., Bellon, M., Wagner, N., Rohr, M., Huttner, A., Blanchard-Rohner, G., Loevy, N., Kaiser, L. & Others  SARS-CoV-2 Viral Load in the Upper Respiratory Tract of Children and Adults With Early Acute COVID-19.  Clinical infectious diseases: an official publication of the Infectious Diseases Society of America, 73(1), 148–150.
  * 23. Polese-Bonatto, M., Sartor, I., Varela, F., Giannini, G., Azevedo, T., Kern, L., Fernandes, I., Zavaglia, G., David, C., Santos, A. & Others  Children Have Similar Reverse Transcription Polymerase Chain Reaction Cycle Threshold for Severe Acute Respiratory Syndrome Coronavirus 2 in Comparison With Adults.  The Pediatric Infectious Disease Journal. 40, e413 (2021)
  * 24. Ade, C., Pum, J., Abele, I., Raggub, L., Bockmühl, D. & Zöllner, B.  Analysis of cycle threshold values in SARS-CoV-2-PCR in a long-term study.  Journal Of Clinical Virology. 138 pp. 104791 (2021)
  * 25. Cendejas-Bueno, E., Romero-Gómez, M., Escosa-García, L., Jiménez-Rodríguez, S., Mingorance, J., García-Rodríguez, J., Group, S. & Others  Lower nasopharyngeal viral loads in pediatric population. The missing piece to understand SARS-CoV-2 infection in children?.  Journal Of Infection. 83, e18-e19 (2021)
  * 26. Bullard, J., Funk, D., Dust, K., Garnett, L., Tran, K., Bello, A., Strong, J., Lee, S., Waruk, J., Hedley, A. & Others  Infectivity of severe acute respiratory syndrome coronavirus 2 in children compared with adults.  Cmaj. 193, E601-E606 (2021)
  * 27. Bellon, M., Baggio, S., Jacquerioz Bausch, F., Spechbach, H., Salamun, J., Genecand, C., Tardin, A., Kaiser, L., L’Huillier, A. & Eckerle, I.  Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) viral load kinetics in symptomatic children, adolescents, and adults.  Clinical Infectious Diseases. 73, e1384-e1386 (2021)
  * 28. Maltezou HC, Magaziotou I, Dedoukou X, Eleftheriou E, Raftopoulos V, Michos A, et al.  Children and adolescents with SARS-CoV-2 infection: epidemiology, clinical course and viral loads.  The Pediatric Infectious Disease Journal. 2020;39(12):e388–e392.
  * 29. L’Huillier AG, Torriani G, Pigny F, Kaiser L, Eckerle I.  Culture-competent SARS-CoV-2 in nasopharynx of symptomatic neonates, children, and adolescents.  Emerging infectious diseases. 2020;26(10):2494.
  * 30. Zachariah P, Halabi KC, Johnson CL, Whitter S, Sepulveda J, Green DA.  Symptomatic infants have higher nasopharyngeal SARS-CoV-2 viral loads but less severe disease than older children.  Clinical Infectious Diseases. 2020;71(16):2305–2306.
  * 31. Jones TC, Biele G, Mühlemann B, Veith T, Schneider J, Beheim-Schwarzbach J, et al.  Estimating infectiousness throughout SARS-CoV-2 infection course.  Science. 2021;doi:10.1126/science.abi5273.
  * 32. Pearl, J.  Causality: Models, Reasoning and Inference.  Cambridge University Press, 2009.
  * 33. Spirtes, P., Glymour, C. & Scheines,  R. Causation, Prediction, and Search, 2nd Edition.  The MIT Press, 2001, 12
  * 34. Mooij, J., Peters, J., Janzing, D., Zscheischler, J. & Schölkopf, B.  Distinguishing Cause from Effect Using Observational Data: Methods and Benchmarks.  Journal Of Machine Learning Research. 17, 1-102 (2016), http://jmlr.org/papers/v17/14-518.html
  * 35. Knollmüller J, Enßlin TA.  Encoding prior knowledge in the structure of the likelihood.  arXiv e-prints. 2018; p. arXiv:1812.04403.
  * 36. Selig M, Bell MR, Junklewitz H, Oppermann N, Reinecke M, Greiner M, et al.  NIFTy - Numerical Information Field Theory. A versatile Python library for signal inference.  Astron. & Astrophys.. 2013;554:A26.  doi:10.1051/0004-6361/201321236.
  * 37. Steininger T, Dixit J, Frank P, Greiner M, Hutschenreuter S, Knollmüller J, et al.  NIFTy 3 - Numerical Information Field Theory: A Python Framework for Multicomponent Signal Inference on HPC Clusters.  Annalen der Physik. 2019;531(3):1800290.  doi:10.1002/andp.201800290.
  * 38. Arras P, Baltac M, Ensslin TA, Frank P, Hutschenreuter S, Knollmueller J, et al.. NIFTy5: Numerical Information Field Theory v5; 2019\.
  * 39. Knollmüller J, Enßlin TA.  Metric Gaussian Variational Inference.  arXiv e-prints. 2019; p. arXiv:1901.11033.
  * 40. Wölfel R, Corman VM, Guggemos W, Seilmaier M, Zange S, Müller MA, et al.  Virological assessment of hospitalized patients with COVID-2019.  Nature. 2020;581(7809):465–469.  doi:10.1038/s41586-020-2196-x.
  * 41. Li, L., White L.F. Bayesian back-calculation and nowcasting for line list data during the COVID-19 pandemic.  PLoS Computational Biology 17, 1-22 (2021,7),
  * 42. Cherief-Abdellatif, B. Consistency of ELBO maximization for model selection.  Proceedings Of The 1st Symposium On Advances In Approximate Bayesian Inference.  96 pp. 11-31 (2019,12,2).
  * 43. Knuth, K., Habeck, M., Malakar, N., Mubeen, A. & Placek, B.  Bayesian evidence and model selection.  Digital Signal Processing. 47 pp. 50-67 (2015), Special Issue in Honour of William J. (Bill) Fitzgerald
  * 44. Kostić A, Leike R, Hutschenreuter S, Kurthen M, Enßlin TA.  Bayesian Causal Inference with Information Field Theory.  Manuscript in preparation (2022).
  * 45. Chen WC, Tareen A, Kinney JB.  Density Estimation on Small Data Sets.  Phys. Rev. Lett.. 2018;121(16):160605.  doi:10.1103/PhysRevLett.121.160605.
  * 46. Silverman BW.  Density estimation for statistics and data analysis. vol. 26.  CRC press; 1986.
  * 47. Sheather SJ.  Density Estimation.  Statistical Science. 2004;19(4):588–597.
  * 48. Liu Q, Xu J, Jiang R, Wong WH.  Density estimation using deep generative neural networks.  Proceedings of the National Academy of Sciences. 2021;118(15).
  * 49. Ferguson TS.  A Bayesian Analysis of Some Nonparametric Problems.  The Annals of Statistics. 1973;1(2):209 – 230.  doi:10.1214/aos/1176342360.
  * 50. Müller P.  Bayesian nonparametric data analysis.  Cham, Switzerland: Springer; 2015.
  * 51. Gelman A.  Bayesian data analysis.  Boca Raton: CRC Press; 2014.
  * 52. Kinney JB.  Estimation of probability densities using scale-free field theories.  Phys Rev E. 2014;90:011301.  doi:10.1103/PhysRevE.90.011301.
  * 53. Kinney JB.  Unification of field theory and maximum entropy methods for learning probability densities.  Phys Rev E. 2015;92:032107.  doi:10.1103/PhysRevE.92.032107.
  * 54. Genton MG.  Classes of Kernels for Machine Learning: A Statistics Perspective.  J Mach Learn Res. 2002;2:299–312.
  * 55. Arras P, Frank P, Haim P, Knollmüller J, Leike R, Reinecke M, et al.. M87* in space, time, and frequency; 2021.


* [139] J. Ren, Z.-Z. Xianyu, and H.-J. He, Higgs Gravitational Interaction, Weak Boson Scattering, and Higgs Inflation in Jordan and Einstein Frames, JCAP 06 (2014) 032, [arXiv:1404.4627].
  * [140] L. Wu, H. Zhang, and B. Zhu, Non-minimal Higgs-gravity coupling and naturalness at the LHC and beyond, JCAP 07 (2019) 033, [arXiv:1901.06532].
  * [141] V.-M. Enckell, K. Enqvist, S. Rasanen, and E. Tomberg, Higgs inflation at the hilltop, JCAP 06 (2018) 005, [arXiv:1802.09299].
  * [142] R. Jinno, M. Kubota, K.-y. Oda, and S. C. Park, Higgs inflation in metric and Palatini formalisms: Required suppression of higher dimensional operators, JCAP 03 (2020) 063, [arXiv:1904.05699].
  * [143] M. He, Y. Mikura, and Y. Tada, Hybrid metric-Palatini Higgs inflation, arXiv:2209.11051.
  * [144] S. Capozziello, T. Harko, T. S. Koivisto, F. S. N. Lobo, and G. J. Olmo, Hybrid metric-Palatini gravity, Universe 1 (2015), no. 2 199–238, [arXiv:1508.04641].
  * [145] T. Eadkhong, P. Dam-O, P. Channuie, and D. Momeni, Nonminimally-coupled warm Higgs inflation: Metric vs. Palatini Formulations, arXiv:2303.00572.
  * [146] A. Berera and L.-Z. Fang, Thermally induced density perturbations in the inflation era, Phys. Rev. Lett. 74 (1995) 1912–1915, [astro-ph/9501024].
  * [147] N. Kaewkhao and B. Gumjudpai, Cosmology of non-minimal derivative coupling to gravity in Palatini formalism and its chaotic inflation, Phys. Dark Univ. 20 (2018) 20–27, [arXiv:1608.04014].
  * [148] D. V. Gal’tsov, Conformal and kinetic couplings as two Jordan frames of the same theory: Conformal and kinetic couplings, Eur. Phys. J. C 80 (2020), no. 5 443, [arXiv:2001.03221].
  * [149] C. Germani and A. Kehagias, New Model of Inflation with Non-minimal Derivative Coupling of Standard Model Higgs Boson to Gravity, Phys. Rev. Lett. 105 (2010) 011302, [arXiv:1003.2635].
  * [150] H. Sheikhahmadi, E. N. Saridakis, A. Aghamohammadi, and K. Saaidi, Hamilton-Jacobi formalism for inflation with non-minimal derivative coupling, JCAP 10 (2016) 021, [arXiv:1603.03883].
  * [151] S. Karydas, E. Papantonopoulos, and E. N. Saridakis, Successful Higgs inflation from combined nonminimal and derivative couplings, Phys. Rev. D 104 (2021), no. 2 023530, [arXiv:2102.08450].
  * [152] L. Amendola, Cosmology with nonminimal derivative couplings, Phys. Lett. B 301 (1993) 175–182, [gr-qc/9302010].
  * [153] S. Capozziello, G. Lambiase, and H. J. Schmidt, Nonminimal derivative couplings and inflation in generalized theories of gravity, Annalen Phys. 9 (2000) 39–48, [gr-qc/9906051].
  * [154] S. V. Sushkov, Exact cosmological solutions with nonminimal derivative coupling, Phys. Rev. D80 (2009) 103505, [arXiv:0910.0980].
  * [155] S. Capozziello and G. Lambiase, Nonminimal derivative coupling and the recovering of cosmological constant, Gen. Rel. Grav. 31 (1999) 1005–1014, [gr-qc/9901051].
  * [156] M. Ostrogradsky, Mémoires sur les équations différentielles, relatives au problème des isopérimètres, Mem. Acad. St. Petersbourg 6 (1850), no. 4 385–517.
  * [157] R. P. Woodard, Ostrogradsky’s theorem on Hamiltonian instability, Scholarpedia 10 (2015), no. 8 32243, [arXiv:1506.02210].
  * [158] I. Dalianis, S. Karydas, and E. Papantonopoulos, Generalized Non-Minimal Derivative Coupling: Application to Inflation and Primordial Black Hole Production, JCAP 06 (2020) 040, [arXiv:1910.00622].
  * [159] J. D. Bekenstein, The Relation between physical and gravitational geometry, Phys. Rev. D 48 (1993) 3641–3647, [gr-qc/9211017].
  * [160] M. Zumalacárregui and J. García-Bellido, Transforming gravity: from derivative couplings to matter to second-order scalar-tensor theories beyond the Horndeski Lagrangian, Phys. Rev. D 89 (2014) 064046, [arXiv:1308.4685].
  * [161] M. Minamitsuji, Disformal transformation of cosmological perturbations, Phys. Lett. B 737 (2014) 139–150, [arXiv:1409.1566].
  * [162] S. Tsujikawa, Disformal invariance of cosmological perturbations in a generalized class of Horndeski theories, JCAP 04 (2015) 043, [arXiv:1412.6210].
  * [163] G. Domènech, A. Naruko, and M. Sasaki, Cosmological disformal invariance, JCAP 10 (2015) 067, [arXiv:1505.00174].
  * [164] J. Sakstein and S. Verner, Disformal Gravity Theories: A Jordan Frame Analysis, Phys. Rev. D 92 (2015), no. 12 123005, [arXiv:1509.05679].
  * [165] C. van de Bruck, T. Koivisto, and C. Longden, Disformally coupled inflation, JCAP 03 (2016) 006, [arXiv:1510.01650].
  * [166] J. Ben Achour, D. Langlois, and K. Noui, Degenerate higher order scalar-tensor theories beyond Horndeski and disformal transformations, Phys. Rev. D 93 (2016), no. 12 124005, [arXiv:1602.08398].
  * [167] A. Escrivà and C. Germani, Beyond dimensional analysis: Higgs and new Higgs inflations do not violate unitarity, Phys. Rev. D 95 (2017), no. 12 123526, [arXiv:1612.06253].
  * [168] I. P. Lobo and G. G. Carvalho, The geometry of null-like disformal transformations, Int. J. Geom. Meth. Mod. Phys. 16 (2019), no. 11 1950180, [arXiv:1707.01784].
  * [169] S. Sato and K.-i. Maeda, Hybrid Higgs Inflation: The Use of Disformal Transformation, Phys. Rev. D 97 (2018), no. 8 083512, [arXiv:1712.04237].
  * [170] D. Gal’tsov and S. Zhidkova, Ghost-free Palatini derivative scalar–tensor theory: Desingularization and the speed test, Phys. Lett. B 790 (2019) 453–457, [arXiv:1808.00492].
  * [171] K. Karwan and P. Channuie, Generalized Conformal Transformation and Inflationary Attractors, Phys. Rev. D 100 (2019), no. 2 023514, [arXiv:1811.03006].
  * [172] A. Delhom, I. P. Lobo, G. J. Olmo, and C. Romero, A generalized Weyl structure with arbitrary non-metricity, Eur. Phys. J. C 79 (2019), no. 10 878, [arXiv:1906.05393].
  * [173] T. Qiu, Z. Xiao, J. Shi, and M. Aljaf, Potential-driven Inflation with Disformal Coupling to Gravity, Phys. Rev. D 102 (2020), no. 6 063506, [arXiv:2004.14712].
  * [174] S. Rasanen and P. Wahlman, Higgs inflation with loop corrections in the Palatini formulation, JCAP 11 (2017) 047, [arXiv:1709.07853].
  * [175] T. Takahashi and T. Tenkanen, Towards distinguishing variants of non-minimal inflation, JCAP 04 (2019) 035, [arXiv:1812.08492].
  * [176] K. Dimopoulos, A. Karam, S. Sánchez López, and E. Tomberg, Palatini R 2 quintessential inflation, JCAP 10 (2022) 076, [arXiv:2206.14117].
  * [177] S. Panda, A. Rana, and R. Thakur, Constant-Roll Inflation in modified $f(R,\phi)$ gravity model using Palatini Formalism, arXiv:2212.00472.
  * [178] A. Bedroya, R. Brandenberger, M. Loverde, and C. Vafa, Trans-Planckian Censorship and Inflationary Cosmology, Phys. Rev. D 101 (2020), no. 10 103502, [arXiv:1909.11106].
  * [179] M. AlHallak, N. Chamoun, and M. S. Eldaher, Natural Inflation with non minimal coupling to gravity in R 2 gravity under the Palatini formalism, JCAP 10 (2022) 001, [arXiv:2202.01002].
  * [180] J. Beltrán Jiménez and A. Delhom, Ghosts in metric-affine higher order curvature gravity, Eur. Phys. J. C 79 (2019), no. 8 656, [arXiv:1901.08988].
  * [181] J. Beltrán Jiménez and A. Delhom, Instabilities in metric-affine theories of gravity with higher order curvature terms, Eur. Phys. J. C 80 (2020), no. 6 585, [arXiv:2004.11357].
  * [182] J. Annala and S. Rasanen, Stability of non-degenerate Ricci-type Palatini theories, arXiv:2212.09820.
  * [183] D. Lovelock, The Einstein tensor and its generalizations, J. Math. Phys. 12 (1971) 498–501.
  * [184] T. Damour, S. Deser, and J. G. McCarthy, Nonsymmetric gravity theories: Inconsistencies and a cure, Phys. Rev. D 47 (1993) 1541–1556, [gr-qc/9207003].
  * [185] J. Annala, Higgs inflation and higher-order gravity in Palatini formulation, Master’s thesis, Helsinki U., 2020.
  * [186] E. Gildener and S. Weinberg, Symmetry Breaking and Scalar Bosons, Phys.Rev. D13 (1976) 3333.
  * [187] C. Dioguardi, A. Racioppi, and E. Tomberg, Inflation in Palatini quadratic gravity (and beyond), arXiv:2212.11869.
  * [188] E. Tomberg and H. Veermäe, Tachyonic preheating in plateau inflation, JCAP 12 (2021), no. 12 035, [arXiv:2108.10767].
  * [189] N. Koivunen, E. Tomberg, and H. Veermäe, The linear regime of tachyonic preheating, JCAP 07 (2022), no. 07 028, [arXiv:2201.04145].
  * [190] K. D. Lozanov and M. A. Amin, Equation of State and Duration to Radiation Domination after Inflation, Phys. Rev. Lett. 119 (2017), no. 6 061301, [arXiv:1608.01213].
  * [191] K. D. Lozanov and M. A. Amin, Self-resonance after inflation: oscillons, transients and radiation domination, Phys. Rev. D 97 (2018), no. 2 023533, [arXiv:1710.06851].
  * [192] F. Dux, A. Florio, J. Klarić, A. Shkerin, and I. Timiryasov, Preheating in Palatini Higgs inflation on the lattice, JCAP 09 (2022) 015, [arXiv:2203.13286].
  * [193] F.-Y. Zhang, Reheating predictions in non-minimally coupled inflationary models with radiative corrections, Phys. Dark Univ. 39 (2023) 101169\.
  * [194] K. Dimopoulos and T. Markkanen, Non-minimal gravitational reheating during kination, JCAP 06 (2018) 021, [arXiv:1803.07399].
  * [195] E. J. Copeland, A. R. Liddle, and D. Wands, Exponential potentials and cosmological scaling solutions, Phys. Rev. D 57 (1998) 4686–4690, [gr-qc/9711068].
  * [196] L. H. Ford, Gravitational Particle Creation and Inflation, Phys. Rev. D 35 (1987) 2955.
  * [197] I. Antoniadis, A. Guillen, and K. Tamvakis, Late time acceleration in Palatini gravity, JHEP 11 (2022) 144, [arXiv:2207.13732].