Decoding the Fundamental Drivers of Phylodynamic Inference

Abstract Despite its increasing role in the understanding of infectious disease transmission at the applied and theoretical levels, phylodynamics lacks a well-defined notion of ideal data and optimal sampling. We introduce a method to visualize and quantify the relative impact of pathogen genome sequence and sampling times—two fundamental sources of data for phylodynamics under birth–death-sampling models—to understand how each drives phylodynamic inference. Applying our method to simulated data and real-world SARS-CoV-2 and H1N1 Influenza data, we use this insight to elucidate fundamental trade-offs and guidelines for phylodynamic analyses to draw the most from sequence data. Phylodynamics promises to be a staple of future responses to infectious disease threats globally. Continuing research into the inherent requirements and trade-offs of phylodynamic data and inference will help ensure phylodynamic tools are wielded in ever more targeted and efficient ways.


Introduction
Phylodynamics combines phylogenetic and epidemiological modeling to infer epidemiological dynamics from pathogen genome data (Volz et al. 2013;du Plessis and Stadler 2015;Baele et al. 2018). Analyses are usually conducted within a Bayesian Markov chain Monte Carlo (MCMC) framework, meaning that the output consists of posterior samples for parameters of interest, such as the basic reproductive number, R 0 (i.e., the average number of secondary infections from the index case in an otherwise fully susceptible population). Input data usually consist of partial-or whole-genome sequences with associated sample collection dates. In the case of birth-death-sampling models (Stadler 2010), both sequence and date data inform the branching of inferred trees by either temporally clustering lineages or via sequence similarity. Internal nodes are assumed to coincide with transmission events, such that they provide information about patterns of transmission that sampling time data alone cannot (Featherstone et al. 2022). Sampling times, or date data, are similar to standard epidemiological time series data whereas sequence data introduce the evolutionary aspect. The widely used birth-death model uses sampling times to infer a sampling rate which is also informative about transmission rates (Stadler et al. 2012;Boskova et al. 2018).
Application of phylodynamic methods has intensified since the onset of the SARS-CoV-2 pandemic. Moreover, larger and more densely sequenced outbreaks are being studied. Although the value of pathogen genome data is now well established, an increasingly pertinent question is whether the inclusion of more sequence data after a point is of diminishing returns for some densely sequenced outbreaks (Hill et al. 2021;Porter et al. 2022). Since the answer to this question will naturally vary with each data set and pathogen considered, we believe it is most appropriately addressed by a method to quantify the individual effects of date and sequence data. This would help broaden our understanding of the phylodynamic tools that now feature in infectious disease surveillance. Such a method also has the potential to help target sampling efforts of future outbreaks for optimization of knowledge gain against resource expenditure.
Earlier work showed that sequence sampling times, referred to here as "date data," can drive epidemiological inference under the birth-death model (Volz and Frost 2014;Boskova et al. 2018;Featherstone et al. 2021). However, each stopped short of proposing a formal method to measure this effect in regular application. The birthdeath model is most applicable to the question at hand since it includes a rate of sampling. The coalescent is another key phylodynamic model, but it typically conditions on sampling dates which therefore precludes a comparison of date and sequence effects. Some coalescent formulations include a sampling rate (Volz and Frost 2014); MBE however, these are used less often than the birth-death or standard Kingman coalescent (Kingman 1977). The Kingman coalescent also assumes a low sampling proportion relative to population size such that its typical formulation would be inappropriate for many densely sequenced outbreaks (Boskova et al. 2018), where the question of the effect of large amounts of sequence data is most relevant.
Building upon these earlier results, we introduce a theoretical framework and a new method to quantify and visualize the effect of sequence and dates for any parameter under the birth-death with continuous sampling. We focus on continuous sampling because it is most relevant to how emerging outbreak data are collected. It also classifies which data source is driving the inference but crucially also indicates whether a binary classification is meaningful. We believe these observations will form a critical addition to the phylodynamic toolkit used to inform public health decisions because they clearly quantify the added-knowledge acquired from genomes in a given analysis.

Isolating Date and Sequence Data
We conduct four analyses for a given data set to contrast the effects of complete data, date data, sequence data, and the absence of both ( fig. 1A). We focus on inferring R 0 , with all other parameters fixed (the sampling proportion, becoming-uninfectious rate, molecular clock rate, and substitution model parameters), but this new approach is applicable to any parameter under the birth-death with any combination of priors. First, we use complete data to fit a birth-death model and infer the posterior distribution of R 0 . This represents the combined effects of dates and sequences. Second, to isolate the effect of date data, we remove sequence information and retain dates, thus integrating over the prior on tree topology. This is traditionally referred to as "sampling from the prior," but this term should be avoided in the context of models where the sampling times are treated as data, such as the birth-death. Third, to isolate the effect of sequence data, we keep sequence data and remove dates. This requires estimation of all sampling dates, analogously to how removing sequence data causes integration over topology. We use a novel Markov chain Monte Carlo (MCMC) operator to estimate dates which are implemented in the feast v17 package for BEAST 2 (Bouckaert et al. 2019). Briefly, the operator adjusts the time between the final sample date and end of the birth-death process such that dates can rescale relative to each other and in absolute time. See supplementary figure S1, Supplementary   FIG. 1. Graphical summary of the process to quantify signal and classify signal drivers. (A) Colored boxes give examples of data under each of the four treatments with letters in brackets giving shorthand notation for each. From left to right: Marginal Prior results from the removal of both date and sequence data. Dates Only includes date data while ignoring sequence data through a constant phylogenetic likelihood. This can be represented as converting all sequence characters to "N" (i.e., alignments of entirely missing data). Full data represents the usual combination of both sequence and date data. This produces a reference distribution from which the Wasserstein metric to other posteriors is calculated. Sequence Only corresponds to the removal and re-estimation of dates while sequence data are retained. (B) Example posterior output for R 0 with color corresponding to each treatment in A. The Wasserstein metric is calculated as the difference in the inverse distribution function of each posterior from Full Data integrated over 0 to 1. Example values for the Wasserstein metric are given in white boxes. (C ) The plane with x and y axes W D and W S and shaded classification regions. r SD is the Euclidean distance from the origin to a point (W D , W S ), with higher values indicating that one or both of data and sequence data drive differing signals from the reference posterior. d SD is the vertical distance from a point (W D , W S ) to the line y = x, with points closer to this line corresponding to more similar data and sequence data effects such that classification is less meaningful. In the example, the distance from the posterior under only sequence data to the full data posterior (W S ) is the smallest, leading to classification as "Seq-Driven."

MBE
Material online for a visual representation. Lastly, and for completeness, we conduct the analysis with both date and sequence data removed. This formally corresponds to the marginal prior conditioned on the number of samples. The resulting Wasserstein metric, W N , is useful for quantifying whether full data offer information in addition to the prior.

Quantifying Data Signal
We employ the Wasserstein metric in one dimension to measure a "distance" between each of the sequence posterior, date posterior, or marginal prior, and the posterior derived from the complete data. We write these distances as W † , with † being D, S, or N for the date, sequence, and marginal-prior distributions, respectively. For example, the Wasserstein distance W D from the date data to complete data posterior is: where F D and F F are cumulative distribution functions for posterior R 0 under date and complete data, respectively. F −1 † then maps from cumulative probability to the corresponding value on the domain of the parameter of interest (e.g., [0, 1] → R 0 ). The units of W † are equivalent to the units of the parameter of interest (e.g., years −1 for the birth rate, becoming uninfectious rate, sampling rate, and unitless for R 0 ). See supplementary figure S2, Supplementary Material online for a visual interpretation of the expression. We used the transport R package to calculate the Wasserstein metric (Schuhmacher et al. 2020).
As in figure 1C, we can now consider a plane where the axes are W D and W S . We classify the data source with the lowest Wasserstein distance from the complete data posterior as contributing most to the posterior from full data. We refer to this as a classifier, although we emphasize that this is not a classifier as in machine learning literature. There is no prediction or statistical modeling of the driving data source beyond min (W D , W S ). In this case, the lines y = x mark the classification boundary as in the shading in figure 1C.
Finally, we can quantify the disagreement in signal between each data source. We define disagreement with respect to the full data posterior, r SD as the magnitude of the vector (W D , W S ���� �→ ) leading to each point in the plane. This is also the radius from the origin to the point. Values near zero indicate that the posteriors under date data, sequence, and complete data are all near identical and classification as date-or sequence-driven is less meaningful. Larger values signify that one or both data sources drive differing posteriors and classification is more meaningful. We also define disagreement without respect to full data d SD = |W S − W D | as a quantification of disagreement between date and sequence posteriors without reference to the full data (supplementary fig. S3C, Supplementary Material online). Visually, this corresponds to the vertical distance to the classification boundary (y = x) such that smaller values correspond to less meaningful classification. r SD and d SD are similar in that when r SD is near-zero, d SD necessarily is too. d SD also accounts for the case where r SD is high, but both date and sequence data have similarly sized effects. In this case, r SD is higher whereas d SD is lower and classification of one or another as driving analysis is inaccurate.

Results
We simulated 600 alignments to explore the differing signals in date and sequence data using the Wasserstein metric. These derive from 100 simulated outbreaks of 500 cases, sampled with proportion 1, 0.5, or 0.05 (n = 25, 250, or 500 accordingly) and used to simulate sequences with an evolutionary rate of 10 −3 or 10 −5 (subs/site/ time). Higher evolutionary rates induce more site patterns and therefore more informative sequence data. We estimated R 0 under each data treatment with all other parameters fixed using a birth-death tree prior. In all analyses, simulated data provided information in addition to the prior (W N > 0.06, supplementary fig. S5, Supplementary Material online). Among the 600 data sets, we observe a mixture of cases where date and sequence data infer similar or dissimilar posterior R 0 . This supports the core assumption that date and sequence data can have differing signals concealed in their combination ( fig. 2). Classification using the Wasserstein metric results in a mix of date-and sequence-driven classifications, supporting that our proposed method is sensitive to differences between data sets ( fig. 3). Most data sets were classified as date driven (372/600), which is consistent with earlier work showing that dates are highly influential under the birth-death (Volz and Frost 2014).

Reliability of Wasserstein Metric
We conducted a subsampling analysis to ensure that Wasserstein metric values reflected differences between date and sequence data, rather than noise alone. For each of the 600 simulated data sets, we subsampled posterior R 0 distributions corresponding to each of the 4 data treatments 100 times, taking a number of samples equal to half the length of the postburnin posterior. This yielded 60,000 subsampled posteriors for which we recalculated Wasserstein values and reclassified each as in the original simulation study. Of the 60,000, subsampled posteriors, only 99 were misclassified across 99 of the 600 original data sets. In other words, misclassification occurred only once out of 100 replicates for 99 of the 600 simulated data sets.
Of the 600 simulated data sets, those where misclassification occurred had substantially smaller differences between Wasserstein distance to the date-only and sequence-only posteriors (d SD ) (supplementary fig. S3, Supplementary Material online). Misclassification occurred where d SD ≤ 8.27 × 10 −2 in the complete R 0 posteriors, with no error above this level. Differences below MBE 8.27 × 10 −2 correspond to a level of difference between data-and sequence-only posteriors where classification is of little significance (i.e., R 0 ± 0.08 assuming both posteriors have similar uncertainty). Classification is wholly reliable above this threshold, validating the trends observed in the simulation study and that the Wasserstein metric is sensitive to differing signals concealed in date and sequence data.

Observations about the Effects of Sequence and Dates
The distribution of W S is more diffuse than W D , meaning the sequence data posterior tends to differ more from full data than date data ( fig. 3A). This again aligns with previous results showing that date data usually drive inference under the birth-death.
Low sequence diversity, measured here in the number of site patterns, seems to preclude sequence data from driving inference ( fig. 3C, supplementary fig. S4, Supplementary Material online). This matches the expectation that fewer site patterns result in less sequence information to inform the posterior. It appears as a necessary but not sufficient condition for a data set to have at least one site pattern per sample in order for analysis to be sequence driven (supplementary fig. S6, Supplementary Material online). On the other hand, the date span does not follow an equivalent trend with lower diversity coinciding with analyses being sequence driven. Here, the relative date span is the time between the first and last sample, divided by the height of the outbreak tree and thus should be indicative of the information content of the dates through the proportion of the outbreak's duration that they capture. The distribution of relative date span appears random across classifications, unlike the distribution of site patterns ( Comparison of full data, date data, sequence data, and marginal prior R 0 for four of the 600 simulated data sets where the true R 0 value is 2.5.

MBE
lower sampling still yields a usable sample size. This matches expectation because lower sampling proportion increases the likelihood of sampling more divergent lineages, resulting in more site patterns to drive inference. One of the core assumptions of phylodynamics is that inferred phylogenetic trees bear a resemblance in their shape and timing to underlying transmission networks (Featherstone et al. 2022). Based on this, we would expect that analyses with less uncertainty in posterior epidemiological estimates have posterior tree distributions that are more divergent. To test this expectation, we calculated tree distance metrics pairwise for each of the 2400 posterior tree distributions, originating from the 600 simulated outbreaks and four data treatments. We sampled 100 trees from each posterior and considered the pairwise topological distance between each. We focus on the mutual clustering information (MCI) introduced by Smith (2020) (supplementary fig. S7, Supplementary Material online), but all observed patterns were repeated when using other metrics including Robinson-Foulds and Informationcorrected Robinson-Foulds. The full data treatment consistently had the lowest pairwise MCI since these analyses incorporated the most information to constrain tree space. The marginal prior, reflecting the least information, reported the highest MCI values. Between these two extremes, sequence-only analyses consistently had a lower

MBE
MCI than date-only treatments (based on visual inspection). This can be explained by sequence data providing information about favored topologies via the phylogenetic likelihood. Conversely, date-only data incorporate no topological information and consistently have higher MCI values. Given that date data tend to drive inference, especially where the evolutionary rate is low and to an extent where sampling proportion is high, this presents strong evidence that the chronology of tips can be equally, if not more informative than the topology of trees. This highlights the importance of curating sampling time data as carefully as sequence data is curated. Moreover, this applies at the phylogenetic level of phylodynamic inference, as distinct from the choice of the phylodynamic model and its parameters.

Australian SARS-CoV-2 Clusters
We analyzed data from two SARS-CoV-2 transmission clusters from 2020 in Australia to demonstrate that datedriven and sequence-driven analyses can arise in practice (table 1). These clusters were chosen out of the 595 identified by Lane et al. (2021) because they possessed the most complete sampling-time data to compare with sequence data. They are also ideal examples of the size and type of genome data set that public health practitioners may seek phylodynamic insight from as each sample most likely reflects transmission stemming from a single, unknown source as supported by contact tracing.
Analysis of the first cluster is classified as sequence driven, with r SD = 0.223 indicating an appreciable difference between the sequence posterior and complete data. d SD = 0.078 adds that date data also drive an effect of similar size, offering the interpretation that both date and sequence data are influential in this analysis. The second cluster is classified as date driven but with r SD = 0.009 and d SD = 0.008. The low r SD value indicates a near-negligible difference between date, sequence, and complete data posteriors. Since r SD is low, d SD is necessarily also low. Due to this, classification is effectively meaningless, and it can be concluded that both date and sequence data drive a highly congruent signal. Moreover, W N for both analyses is more than double each of W S and W D , which affirms that both sources of data contribute to the posterior deviating from the prior and are therefore informative with respect to the prior in both analyses.

H1N1
The above simulation study and SARS-CoV-2 empirical data consist of analyses where all parameters are fixed except R 0 . However, this extent of prior certainty is rare in practice. In particular, there exists an inherent trade-off between prior certainty in the evolutionary rate and the strength of signal due to sequence data. Moreover, this affects the inference of the sampling proportion, which relies on sequence information to inform the evolutionary distance between samples and tMRCA, in turn informing the proportion of lineages captured in sampling. To explore the effects of these trade-offs, we consider North-American H1N1 data from the 2009 swine flu pandemic originally studied by Hedge et al. (2013). This data set was chosen because it represents a small sample from a much larger outbreak, which helps to ensure to ensure the presence of sequence information with which to test the effects of prior uncertainty in the evolutionary rate.
We again fit a birth-death process, estimating two R e values before and after early July 2009 (R e 1 and R e 2 , respectively capturing the exponential and postexponential phases of the outbreak) and the sampling proportion (p). Analyses were conducted with either a fixed evolutionary rate, 4 × 10 −3 subs/site/year following Hedge et al. (2013), or a uniform prior U(10 −4 , 10 −2 ), and under each of the four data treatments defined above. We placed a noninformative β(1, 1) prior to the sampling proportion and fixed the becoming-uninfectious rate (δ = 365.25(days)/4), such that comparing analyses with a fixed or uniform prior on the evolutionary rate reveals the effects of uncertainty in evolutionary rate on R e and sampling proportion. Posterior distributions for R e and p are presented in figure 4 and Wasserstein values in table 2.
The patterns in posterior R e under each treatment verify the trade-off between estimating tip dates and evolutionary rate. This is reflected in higher W S values when using a prior on the evolutionary rate than when fixing it (table 2), which corresponds to posterior R e under Full Data overlapping less with the posterior from sequence data. Moreover, the sequence data posterior is wider where the uniform clock prior is used, reflecting that more sequence information is devoted to ascertaining the evolutionary rate where the prior is more diffuse, which would otherwise further ascertain R e (fig. 4A). There is also a disparity between the posterior tMRCA between each clock prior under Full Data (Only the fixed prior tMRCA is shown in fig. 4A), which is expected given the well-known nonidentifiability between the evolutionary rate and tMRCA.
Estimates of sampling proportion appear to be strongly sequence driven regardless of the evolutionary rate prior ( fig. 4B). The posterior sampling proportion is two orders of magnitude lower for the sequence data and full date treatments, both containing sequence information, than for date data alone. These estimates are far more in line with expectations too, given the disparity between the

MBE
H1N1 sample size (n = 100), and scale of the 2009 H1N1 pandemic. This is expected because the diffuse β(1, 1) prior placed on p made sequence information essential to informing the evolutionary distance between samples and driving estimates of sampling proportion lower. Conversely, estimates with date-only data are much higher and more uncertain owing to the lack of sequence information. The evolutionary rate prior drives a similar trend as for R 0 , but its effect is diminished by the major effect of sequence data. Altogether, this highlights that different parameters can be drawn from date and sequence data to different extents. We do not present the marginal-prior data in figure 4 because each corresponding posterior is uninformatively broad (i.e., effectively uniform over the prior domain). However, this too supports our assertion of a trade-off between sequence information and posterior certainty in sampling proportion because of its contrast to the simulation studies and SARS-CoV-2 data sets where sampling proportion was fixed. The only information included in the marginal-prior treatment is the sample size, and this is most informative when the sampling proportion is fixed. However, in the H1N1 analyses, the diffuse prior on sampling proportion means sample size alone is insufficient to drive an informatively narrow posterior and sequence is the primary driver of inference.
Taken together, the results of the H1N1 data set demonstrate that uncertainty in the evolutionary rate prior partly determines the strength of effect for sequence data, since any information in sequence data is used in shifting from prior to posterior. In addition, sampling proportion demands sufficient sequence information for Counts are rescaled by a factor of 10 −5 . Vertical marks along the date axis represent sampling times in the data set. The vertical dashed line marks the change time between R e1 and R e2 . Asymmetrical violin plots give the posterior density for R e under each data treatment with the marginal prior omitted. Distributions on the left correspond to analyses with a fixed evolutionary rate and those on the right correspond analyses with a uniform prior on the evolutionary rate. The more diffuse clock prior decreases the influence of date data. (B) An asymmetrical violin plot of the posterior sampling proportion with the left and right corresponding to fixed and uniform clock priors identically to A. Sequence data drive inference for both clock prior configurations.

MBE
inference, and this cannot be substituted by sampling times as readily as has been shown for R e (Volz and Frost 2014;Featherstone et al. 2021).

Discussion
The results of our simulation study add clarity and practical understanding to previous results showing that sampling times contribute substantially to phylodynamic inference under the birth-death (Volz and Frost 2014;Featherstone et al. 2021). We furthermore demonstrate that sequence data are not always secondary in influence and can drive inference of R 0 in some instances. This affirms the sensitivity of the birth-death to the signal encoded in sequence data. We demonstrate that the number of site patterns is key to enabling sequence-driven inference but is not the sole factor. Over-sampling from a pathogen population can also dilute sequence information and pivot analyses to being date driven. The tendency for date data to drive inference as site patterns decrease and sampling proportion increases can be explained by the reduction in uncertainty that each data source offers. Dates impose hard bounds by restricting tree space to a subset of trees that agree with the chronology of sampling times. Conversely, sequence data inform topology through phylogenetic likelihood but do not definitively constrain tree space in the same way as date data. The net result is that sequence data must explore a larger tree space than date data. Thus, as the information in sequence decreases, such as via intensified sampling lowering the number of site patterns per sample, the disparity between the sequence and date data tree spaces increases and dates more often drive analyses. See supplementary figure S8, Supplementary Material online for visual explanation.
We note that having sufficient sequence information to drive inference is a related but ultimately different concept to the phylodynamic threshold, which is defined as the sampling timespan needed for a pathogen to display temporal signal (Duchene et al. 2020). Our analyses show that surpassing the phylodynamic threshold (i.e., having temporal signal present in the data) is generally not a sufficient condition for sequence-driven inference. Conversely though, we demonstrate with the H1N1 data set that increased prior certainty in the evolutionary rate increases the effect of sequence data in the analysis. This means that an assessment of the phylodynamic threshold is highly useful in practice to promote sequence-driven inference, but only retrospective analysis using the methods introduced here can conclusively appraise the contribution of the sequence data.

Practical Insights
Our central finding is that sampling time data increasingly drive phylodynamic analyses as data sets grow in the proportion of sampled cases. Although information from sample times can be a boon to phylodynamic analyses when the sampling process is properly modeled, dependence on this information can become a vulnerability as sampling processes deviate from the assumptions of phylodynamic models over time and space. Thus, the phylodynamic value of a data set must take into consideration the amount of sequence and date information available to direct inference over tree space, rather than the size of the data set alone.-Bigger is not automatically better. To this end, we recommend a cautious approach to sampling, with the recognition that each additional sample is increasingly costly to the analysis in requiring sequence data to navigate a disproportionately large tree space in comparison to sampling times. Therefore, samples that are both identical in sampling time and sequence, such as from a super spreading event, should be avoided when curating phylodynamic data sets. However, we make a distinction between sampling and data set curation. We do not at all advocate that sampling be reduced to avoid recovering clonal samples since these are undeniably hallmark of many outbreaks. Rather, when larger data sets are being curated specifically for phylodynamic analysis, care should be taken to avoid excessive duplication where possible. As a basic start point, having at least one site pattern per sample greatly increases the likelihood of a sequence-driven analysis (supplementary fig. S6, Supplementary Material online). Other practical considerations are summarized in table 3.
Repositories of pathogen genome sequences and curated data sets, such as GISAID and Nextstrain are indispensable assets to the phylodynamics research program, but careful curation of data sets with respect to both sampling times and sequence quality is essential to provide valuable phylodynamic results. Just as it is valuable to surpass the phylodynamic threshold in sampling, it is of equal value not to surpass the notional threshold of redundancy in sampling from phylodynamic studies.
Moreover, different parameters require date and sequence information for accurate inference to differing extents. The ability to infer sampling proportion is integral to the value of phylodynamics as it allows for estimation of how many unknown cases there may be in an outbreak. We show sampling proportion is highly sequence driven, which necessitates both sampling to maximize sequence information relative to date information, and measurement of the effect of sequence data in analyses to comment on the reliability of estimated sampling proportions.
Fixing or constraining parameters in analyses where it is justifiable is also critical to ensuring sequence information is reflected in posterior results. For example, fixing the evolutionary rate is ideal if possible. Providing topological constraints between samples would also help concentrate sequence information in the data by reducing the tree space analyses must traverse. These guidelines are complementary to recent work due to Louca et al. (2021), showing that at least one of the unidentifiable parameters of the birth-death model (R e , δ, and p) needs to be fixed for accurate inference from the state space of analyses (i.e., the tree and generative parameters).

MBE
A final practical message is that a phylodynamically ideal microbe would, for example, evolve at a fixed substitution (subs/genome/time) on an infinitely sized genome with consistent clock-like evolution to uniquely barcode every transmission event. Clearly, no such microbe exists. Phylodynamics instead operates over microbes with genomes between 10 4 and 10 8 bases in size and evolution rates from 10 −2 to 10 −8 (subs/site/year) (Biek et al. 2015). In this view, there will never be a perfect data set or pathogen to analyze. The challenge instead becomes sampling and constraining analyses parsimoniously enough to ensure the epidemiological signature encoded in sequence data is returned as well as measuring this effect, such as by the methods developed here.

Application
In general, we recommend that users follow the same workflow as presented here to tease apart date and sequence signals under the birth-death model. That is, separate data into each of the four treatments (marginal prior is optional), perform analyses for each, and compare posteriors using the Wasserstein distance. However, we note that as with all Bayesian phylodynamic analyses, convergence the posterior is never guaranteed, and our workflow may not be feasible such as in the case of R e 2 for the H1N1 data set. In these cases, the practical guidelines (summarized in table 3) still offer a valuable reference for the design and interpretation of phylodynamic analyses under the birth-death model.
The question remains as to what users should do if analyses appear insufficiently sequence driven. This depends heavily on the use case. For example, work involving ancient DNA or fossil taxa is unlikely to allow much scope for adding or removing data to enrich sequence signal. In these cases, our workflow can still be applied to offer useful insight into the information driving analysis, even if no new data are available to increase sequence information.
Our workflow is especially pertinent to public health applications involving large data sets that inform policy and decision-making. Here, it is especially important to test for the relative impacts of date and sequence as these directly inform or understanding of the data driving evidence-based decisions. For example, if results are heavily date driven, then less weighting may be given to the posterior sampling proportion and by extension, phylodynamic conclusions about the number of total number of cases. Consider the order-of-magnitude difference between date-only and sequence only/full data posterior sampling proportion for the H1N1 data set as an example. In addition, phylodynamic analyses may be employed to garner understanding about the origin of transmission through the posterior distribution of trees. Again the level of credence in results can be judged relative to the impact of sequence in analyses. At this point, researchers may opt to include extra sequence data to enrich sequence information, or exclude duplicate samples to hone in on samples that drive existing sequence signal. Fundamentally though, our workflow offers the chance to develop a more complete picture of the drivers of inference, so analyses and data can be better understood in the context from which they originated. From this new perspective, users can decide to add or remove data depending on their objective and data availability.

Future Directions
The work presents the beginning of deeper interrogations of the drivers of phylodynamic inference. Such an understanding is critical as phylodynamics becomes a mainstay of infectious disease surveillance.
Here, we have presented a method, based on the Wasserstein metric, for post hoc analysis of the drivers of phylodynamic inference. A critical future direction would be to develop predictive methods with explicit mathematical formulae for the information content of phylodynamic data. Such a toolkit would enable a theory of optimal sampling for phylodynamic inference to be explicitly described. However, we stress that a method to assay the drivers of

MBE
inference, such as the work presented here, needed to be developed first to help measure the effectiveness of any future predictive methods.

Data Archival
All scripts used to simulate and analyze data are available at https://github.com/LeoFeatherstone/phyloDataSignal. git. The Feast package, containing our date estimator, is available at https://github.com/tgvaughan/feast.git.

Simulation Study
We simulated 100 outbreaks of 500 cases under a birthdeath process using the Tree-Sim R package (Stadler 2019). The birth rate was set to 2.5, death rate 1, corresponding to R 0 = 2.5. Sampling probability was set to p = 1, resulting in trees with 500 tips. We then extended this to a set of 300 outbreaks by sampling again with probability p = 0.05 and p = 0.5, resulting in trees of 25 and 50 tips. We used a consistent seed such that each outbreak with p = 0.05 or p = 0.5 corresponds to a subsample of another with p = 1, allowing us to assess the effect of sampling proportion on inferring W † . For each outbreak, we used Seq-Gen to sequence alignments of length 20,000, which is roughly average for RNA viruses (Rambaut and Grass 1997;Sanjuán et al. 2010). We set an Jukes-Cantor model with evolutionary rate set to either 10 −3 or 10 − 5 subs/site/time using Seq-Gen (Rambaut and Grass 1997).
Our choice of evolutionary rates allows us to compare the effects of higher and lower sequence information with the former corresponding to about 20 substitutions per infection, and 0.2 for the latter. The above resulted in 600 alignments to test in the four treatments described above. We analyzed each under a birth-death model using BEAST v2.6 Bouckaert et al. (2019) with a Uniform[0, 5] prior for R 0 and all other parameters set to the true value for simplicity and to disentangle any impacts of parameter nonidentifiability (Louca et al. 2021).

Empirical Data
All empirical analyses were conducted using BEAST v2.6.6 (Bouckaert et al. 2019) and MCMC chains were run until all parameters of interest had ESS above 200 following burnin. The only exception is the marginal prior for the H1N1 data with a fixed clock. Achieving ESS > 200 for R e 2 required a prohibitively long chain. The lack of convergence in R e 2 could be attributed to an inconsistent signal for branching to inform R e 2 . However, this outcome is reflective of our experience that analyses under the no-data, sequence-only, and marginal-prior data treatments are sometimes nonconvergent. In these cases, we assert that the value of this work persists in the practical guidelines elucidated (table 3). Nevertheless, we found that convergent analyses required computation time similar to that for complete data, which naturally varies with the computational platform used. Analyses using date-only data are usually faster since they do not require calculation of the phylogenetic likelihood. Sequence-only and marginal-prior analyses are usually slower, which is consistent with prior knowledge that dates usually drive inference under the birth-death.

SARS-CoV-2
We analyzed two similar SARS-CoV-2 data sets taken from Lane et al. (2021). They consisted of 112 and 188 samples, respectively. We analyzed each data set under the four conditions above. In each, we placed a Lognormal(mean = 1, sd = 1.25) prior on R 0 and an Inv − Gamma(α = 5.807, β = 346.020) prior on the becoming-uninfectious rate (δ) following estimates of the duration of infection (=1/δ) in Lauer et al. (2020). We also fixed the sampling proportion to p = 0.8 since every known Victorian SARS-CoV-2 case was sequenced at this stage of the pandemic, with a roughly 20% sequencing failure rate. We also placed an Exp(mean = 0.019) prior on the origin, corresponding to a lag of up to 1 week between the index case and the first putative transmission event.

H1N1
We analyzed North American H1N1 swine flu samples taken from the 2009 swine flu pandemic studied by Hedge et al. (2013). We fit a birth-death model with two intervals for R e , before and after July 07, 2020. This date sits in a breakpoint in sampling intensity and was chosen to separate transmission dynamics between the early and late stages of the sampling timespan.
We fixed the becoming infectious rate to 91, corresponding to a duration of infection of 4 days. We also placed a β(1, 1) prior on the sampling proportion and used an HKY with four gamma categories. We compared the effect of two prior configurations for the evolutionary rate, using either a fixed rate a 4 × 10 −3 subs/site/year or a uniform prior (U(10 −4 , 10 −2 )). All other priors were left as defaults.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.