Implementation
Overall structure
The general implementation of o2geosocial follows the
structure of outbreaker2
and incorporates C++ functions into a R framework using Rcpp. The
main function of the package, called outbreaker(), uses
Monte Carlo Markov Chains (MCMC) to sample from the posterior
distribution. For each case, it infers the infection date, the infector,
and the number of missing generations between the case and their
infector. It takes five lists as inputs: i) ‘moves’, ii) ‘likelihoods’,
iii) ‘priors’, iv) ‘data’, and v) ‘config’. These five lists can be
generated and modified using the functions custom_moves(),
custom_likelihoods(), custom_priors(),
create_config() and outbreaker_data().
Examples of how these functions are used to customize the model will be
presented in the Tutorial section.
The package o2geosocial was developed for the study of outbreaks where full information on the onset date, location and age group of the cases is available. It builds upon outbreaker2 by adding the effect of the location and the age-stratified contact data on the reconstruction of the transmission clusters. However, unlike outbreaker2, o2geosocial does not take genetic sequences as input. Genetic information (or genotype) is only used to exclude connections between cases, i.e. two cases cannot be from the same cluster if they belong to different genetic groups. Therefore, o2geosocial is adapted to reconstructing transmission clusters from large datasets (up to several thousand cases) where genetic sequencing is not informative, either because the mutation rate of the virus is slow, or because sequencing is scarce. On the other hand, outbreaker2 is more adapted to reconstructing transmission chains if genetic sequences are informative and available for a high proportion of the cases.
In o2geosocial the genetic component of the likelihood is a binary value depending on the genotype of the cases. There can be only one genotype reported per transmission tree, hence the number of trees will be at least equal to the number of unique genotypes reported in the dataset. The genotype of a tree T is the genotype reported for at least one of the cases belonging to T. For each genotyped case \(i_{gen}\) and at every iteration, only cases from trees with the same genotype as \(i_{gen}\), or without reported genotype should be listed as potential infectors. Therefore, ungenotyped cases belonging to a tree with a different genotype will not be potential infectors of \(i_{gen}\) at this iteration. Ungenotyped cases can be placed on any tree.
Routinely collected surveillance data may include thousands of cases
from unrelated outbreaks. Therefore, it was crucial to develop the
algorithm to be computationally efficient. We added a pre-clustering
step to the function outbreaker() to reduce the pool of
potential infectors for each case, which can greatly reduce the
computing time of the MCMC (Figure 1). In the pre-clustering step, the
potential infectors for each case are listed. Cases with overlapping
potential infectors, and their potential infectors, are grouped
together. Cases from different groups cannot infect one another. The
group each case belongs to is listed in the variable ‘is_cluster’, which
is an element of the list returned by the function
outbreaker_data(). We used three criteria to group together
cases that may be connected: For each case \(i\), only cases reported in a radius of
\(\gamma\) km, less than \(\delta\) days before \(i\), and from similar or unreported
genotype can be classified as potential infectors of \(i\). The parameter \(\gamma\) and \(\delta\) are defined as inputs in the
function create_config().
Cases classified in the same group after the pre-clustering stage may
still be unrelated (e.g. if several importations were reported
simultaneously in a nearby region). Because of this, we defined a
likelihood threshold \(\lambda\) which
allows connections to be discarded as unlikely. If the likelihood of
connection between cases \(i\) and
\(j\) is less than \(\lambda\), we consider it to be more likely
that \(i\) was an import and unrelated
to \(j\). In o2geosocial, the
variable \(\lambda\) can either be an
absolute (the log-likelihood threshold will be \(log(\lambda)\)) or relative value (a
quantile of the likelihood of all connections in all samples), and is
defined by the variables ‘outlier_threshold’ and ‘outlier_relative’ in
create_config(). If the initial number of importations is
too small, unlikely connections between cases can alter the inferred
infection dates of cases and bias the generated transmission trees.
Therefore, we first run a short MCMC and compute \(n\), the minimum number of connections with
a likelihood lower than \(\lambda\) in
the samples; then we add \(n\) imports
to the starting tree and run a long MCMC. Finally, we remove the
connections with likelihood lower than \(\lambda\) in the final samples and return
the infector, infection date and probability of being an import for each
case (Figure 1).
Likelihood and priors
The functions custom_likelihoods() and
custom_priors() can be used to edit each component of the
likelihood and priors.
Genotype
The genetic component of the likelihood that a case \(i\) of genotype \(g_i\) was infected by a case \(j\) belonging to the tree \(\tau_j\) is defined as a binary value: \[G(g_{i},g_{\tau_{j}})= \left\{\begin{array}{l} 1\text{ if }g_{i}\text{ unknown}\\ 1\text{ if }g_{\tau_{j}}\text{ unknown }\\ 1\text{ if }g_{i}\text{ and }g_{\tau_{j}}\text{ both known and }g_i=g_{\tau_{j}} \\ 0\text{ otherwise} \end{array} \right. \]
Therefore, ungenotyped cases may not be potential infector if they belong to a tree that contains another genotype than \(g_i\). The algorithm must consider the genotype of each tree containing potential infectors.
Conditional report ratio
As in the packages outbreaker and outbreaker2, we allow for missing cases in transmission chains. The number of generations between linked cases \(i\) and \(j\) \(\kappa_{ji}\) is equal to one if \(j\) infected \(i\). We define \(\rho\) as the conditional report ratio of the trees. The conditional report ratio is not the same as the overall report ratio of an outbreak because entirely unreported clusters, or missing cases before the ancestor of a tree will not change the value of \(\rho\). Only unreported cases within transmission chains will impact the conditional report ratio. By default, the probability of observing \(\kappa_{ji}\) missing generation between \(i\) and \(j\) from the conditional report ratio \(p(\kappa_{ji} | \rho)\) follows an exponential distribution.
The conditional report ratio is estimated during the MCMC runs using
a beta distribution prior. The two parameters of the beta prior can be
changed using the variable prior_pi in
create_config() (default to \(Beta(10, 1)\)).
Time
The time component of the likelihood is similar to what is used in
the default version of outbreaker2.
The probability of \(t_i\) being the
infection date of the case \(i\)
reported at time \(T_i\) depends on the
distribution of the incubation period \(f\). The incubation period should be
defined using the variable f_dens in the function
outbreaker_data().
The probability that \(i\) was
infected by \(j\) given their
respective inferred dates of infection \(t_i\) and \(t_j\) is defined by the generation time of
the disease \(w^{\kappa_{ji}}(t_i-t_j)\) (variable
w_dens in outbreaker_data()), and the number
of generations \(\kappa_{ji}\) between
\(i\) and \(j\). The operator \(w^{\kappa_{ji}}\) was defined as \(w^{\kappa_{ji}}= \Pi_{\kappa_{ji}}w\),
where \(\Pi\) is the convolution
operator.
Space
In o2geosocial, we introduce a spatial component to the likelihood of connection between cases. By default, we implemented a gravity model to illustrate the probability of connection between two regions \(k\) and \(l\). Gravity models depend on the population sizes \(m_k\) and \(m_l\), and the distance between regions \(d_{kl}\). Given spatial parameters \(a\) and \(b\), the probability that a case in the region \(k\) was infected by a case reported in \(l\) is \(s(k,l)\):
\[s(k,l)= \frac{p_{kl}}{\Sigma_hp_{hl}}=\frac{F(d_{kl}, b)* m_k^a*m_l^c}{\Sigma_h(F(d_{hl}, b)*m_h^a *m_l^c)}=\frac{F(d_{kl}, b)* m_k^a}{\Sigma_h(F(d_{hl}, b)*m_h^a)}\]
If we use an exponential gravity model, \(F(d_{kl},a)=e^{-b*d_{kl}}\); for a
power-law gravity model: \(F(d_{kl},a)=(\frac{1}{d_{kl}})^b\). The
exponential gravity model has been shown to perform well to represent
short-distance mobility patterns. As o2geosocial aims at
reconstructing transmission in a community or a region, the default
model in the function outbreaker() is an exponential
gravity model. The type of gravity model can be changed by setting the
parameter spatial_method to “power-law”:
create_config(spatial_method = "power_law"). Other mobility
models can be used by developing modules. We give an example of module
in the tutorial where we replace the exponential gravity model by a
Stouffer’s rank model.
The parameters \(a\) and \(b\) are estimated during the MCMC run. This requires re-computing the probability matrix twice per iteration and can be time-consuming. Therefore, if either \(a\) or \(b\) is not fixed, we allow for a maximum of 1 missing generation between cases (\(max(\kappa_{ji})=2\)) and only compute \(s^1(k,l)\) and \(s^2(k,l)\) for regions that could potentially be connected. By default, the prior distribution of \(a\) and \(b\) are uniform.
Age
Social contact matrices provided by large scale quantitative
investigations such as the POLYMOD
study can be used to calculate the probability of contact
between cases of different age groups in different countries. These
social contact matrices are imported using the R package socialmixr.
In o2geosocial, given the age group of each case \(\alpha_{(1,..,N)}\) and the age-stratified
social contact matrix, we introduced \(a^{\kappa_{ji}} (\alpha_i,\alpha_j)\), the
probability that a case aged \(\alpha_j\) infected a case aged \(\alpha_i\). This corresponds to the
proportion of contacts to \(\alpha_i\)
that came from individuals of age \(\alpha_j\). The contact matrix can be
modified using the variable a_dens in
outbreaker_data().
Overall likelihood
The overall likelihood that a case \(i\) was infected by the case \(j\) is equal to \(L_i(t_i, j, t_j, \theta ) = log(f(t_i - T_i)) + L_{ji}(t_i, t_j, \theta)\) where \(L_ji(t_i,t_j,\theta)\) is the likelihood of connection between \(i\) and \(j\) and is defined as:
\[L_ji(t_i,t_j,\theta)=log(\rho(\kappa_{ji}|\rho)*w^{\kappa_{ji}}(t_i-t_j)*a^{\kappa_{ji}} (\alpha_i,\alpha_j )*G(g_i,g_{\tau_j})*s^{\kappa_{ji}} (r_i,r_j | a,b))\]
Tree proposals
At every iteration of the MCMC, a set of movements is used to propose
an update of the transmission trees. This update is then accepted or
rejected depending on the posterior density of the new trees. In
o2geosocial, eight default movements are tested at each
iteration. Three of them were already part of outbreaker2
and were not modified (cpp_move_t_inf() to change the
infection date of cases; cpp_move_pi() to change the
conditional report ratio; cpp_move_kappa() to change the
number of generations between connected cases); two were edited to scan
of the transmissions trees to prevent from having different genotypes in
the same tree (cpp_move_alpha().
cpp_move_swap_cases()); The last three are new movements
and were not part of outbreaker2:
cpp_move_a()andcpp_move_b()propose new values of the spatial parameters \(a\) and \(b\) and compute the matrix of probability of connection between regions.cpp_move_ancestor(): An ancestor is defined as the first case of a transmission tree. As only non-ancestors are moved incpp_move_alpha(), we added this function to ensure good mixing of the ancestors of the transmission trees. For each ancestor \(i\), an index case is drawn from the poll of potential infectors, while another link is randomly picked and deleted. We then compare the new posterior value, and accept or reject the change.
