Model definitions and assumptions

All models were designed to be applied to a meta-population of a multi-strain opportunistic pathogen. Each population i corresponded to a matched sample of microbes from asymptomatic carriers, or the environment, and a sample from disease. It was assumed all cases of disease emerged independently from the carried or environmental population. Across the meta-population, the microbes were classified according to their type j, strain k, or the combination of both, j,k. The “type” may be defined by any phenotypic categorisation; analyses of S. pneumoniae typically use serotyping. The “strain” can represent any coherent genetic subdivision of the population; for instance, analyses of S. pneumoniae can use the Global Pneumococcal Sequence Clusters (GPSCs) [37], but alternatively-defined “clades” or “lineages” could be employed instead. If x denotes any one of these classifications, then the carriage rate in population i, ρ_i,x, was estimated based on the total number of samples (including negative results), η_i, and the number of samples positive for isolates of category x, c_i,x:

Modelling the number of positive samples using a Binomial distribution follows the precedent of Weinberger et al [39], and assumes all carriage samples are independently drawn from the circulating population.

A set of related statistical models were used to estimate the progression rate [2], ν, the hazard of developing a disease, given carriage of a microbe over a unit of time. As ν was assumed to be constant for each type x in each study i, the number of observed cases of disease caused by x in i, d_i,x, was modelled as the product of the number of potential hosts in i, N_i; the proportion carrying x in i, ρ_i,x; the per unit time probability of x causing disease, ν_x; and the duration of the study, t_i. Hence, the expected number of isolates from disease, , in all models was proportional to all these quantities:

Some models also adjusted the expected d_i,x between i through a scaling parameter, γ_i, that was relative to a reference population. This corrected for differences between studies both in the actual hazard of developing a disease, due to host population variation, and differences in the comprehensiveness of surveillance across the study populations. Hence, the expected number of isolates from disease, , in such models was proportional the quantities:

All these parameters were assumed to be constant within a study over t_i. As all disease cases were assumed to develop independently, they were modelled as a Poisson process. However, relevant unmodelled variation, or oversampling of local transmission chains in either carriage or disease surveillance [56], could increase heterogeneity in the observed isolate counts in either component of the paired samples. Hence d_i,x was additionally modelled as following a negative binomial distribution, in which overdispersion was quantified by the precision parameter, ϕ, relative to the mean, μ. Lower ϕ values indicate greater overdispersion, as the variance of the distribution was equivalent to:

Models of disease prevalence

Four different structures were used to model . For each, two different versions were fitted to the data. In the simpler version, the observed disease case counts were assumed to follow a Poisson distribution:

The more complex version assumed a negative binomial distribution of the observed disease case counts:

Null model

In this simplest model, the progression rate was independent of the pathogen’s type (ν_x = ν for all x) and host population (γ_i = 1 for all i). Therefore d_i,x is expected to be correlated with c_i,x:

Type-specific model

This model, similar to that Weinberger et al applied to a single population [39], allowed for variation in progression rate between types, but still assumed host populations to be homogeneous (γ_i = 1 for all i). Therefore d_i,x depends on both the carriage prevalence and progression rate of x:

Study-adjusted model

In this model, variation in disease prevalence between studies reflected differences in host population and surveillance, rather than differences between types. Hence the progression rate was independent of the pathogen’s type (ν_x = ν for all x), but varied between i due to γ_i:

Study-adjusted type-specific model

This model allowed for both variation in the progression rate between types x and study populations i:

Joint modelling of strain and type progression rates

For datasets in which information was available on both type (for S. pneumoniae, serotype) j and strain k, ν_x was modelled in five different ways:

Model priors and implementation

The characteristics of each study population (η_i, N_i, t_i) and observed counts of type x in carriage (c_i,x) and disease (d_i,x) were assumed to be accurate. As each recovery of a type x isolate from carriage contributes to the estimation of ρ_i,x, even if an individual carries multiple types, the prior distribution was:

This means each type’s carriage was independently estimated, and the sum of carriage prevalences (∑_x ρ_i,x) may be greater than one, if there is high level of multiple carriage in a population. The beta distribution was used, as it is the conjugate prior of the binomial distribution, although setting the scale and shape parameters to one make it equivalent to a uniform distribution bounded by zero and one.

Invasiveness has been estimated to vary over orders of magnitude in S. pneumoniae [29], and therefore a uniform prior was placed on the logarithm of v_x in all models. The lower bound was set at one IPD case per million carriers per unit time, as the largest studies used surveillance of host populations consisting of tens of millions of individuals (S1 Table). The upper bound was ten cases per carrier per unit time, which would effectively correspond to an obligate pathogen for a bacterium with the carriage duration of S. pneumoniae [57]:

As the scaling factors were relative measures of disease prevalence between populations, γ_i for the population with the largest sample size was fixed at one. The scaling factor for each other population i was allowed to vary higher or lower. The Cauchy distribution was assumed as the prior for the logarithm of γ_i, as this describes the ratio of two normally-distributed random variables. However, the Hamiltonian Monte Carlo algorithm used for model fitting does not efficiently sample heavy-tailed distributions, and therefore the prior was reparameterised to use a random variable Γ_i with a unit-sized uniform prior [58]:

The Γ_i values were then transformed to generate γ_i values with a Cauchy distribution using the location (μ_study = 0) and scale (τ_study = 2) parameters:

The τ_study value allowed for variation in observed incidence over multiple orders of magnitude between datasets, ensuring the model makes allowance for substantial heterogeneity in observed disease burdens [59].

A similar reparameterization of the Cauchy distribution was used for the models in which type and strain background independently contributed to a genotype’s progression rate (the type-determined strain-modified, and strain-determined type-modified, progression rate models). The categorisation modelled as determining the progression rate (v_x) had the same prior as when a single factor determined the progression rate. This was multiplied by a second value, v_y, determined by the classification modifying the progression rate. The prior distribution for the logarithm of this value was a truncated Cauchy, symmetrical around zero. This represented the expectation that the modification of the progression rate would be small (approximately one) in most cases, but may be substantial in rare instances. The truncation was found to be necessary for efficient sampling with the Hamiltonian Monte Carlo algorithm. This was parameterised with a random variable Υ_y, sampled from a uniform distribution with narrower boundaries than those for the analogous Γ_i variable:

The Υ_y values were then transformed to generate v_y values using the location (μ_mod = 0) and scale (τ_mod = 0.5) parameters:

This meant v_y values could vary over three orders of magnitude (between 0.0313 and 32.0), while enabling models to converge with reasonable MCMC lengths.

The precision parameter of the negative binomial distribution, ϕ, was sampled from an exponential distribution. Values of ϕ that are large relative to suggest there is little evidence of overdispersion, and therefore these parameter values should be associated with small prior probabilities. As must always be close to, or above, one in well-powered case-carrier studies, an exponential rate parameter of one was used to penalise the likelihood of the negative binomial model when ϕ was above this value. Yet this broad, always positive, distribution enabled estimation of smaller, albeit non-zero, ϕ values for representing highly overdispersed data:

Model parameters are summarised in Table 1. All models were implemented using R [60] and Rstan [61]. Each model was fitted using two Hamiltonian Markov chain Monte Carlo (MCMC) evaluations for 25,000 iterations, of which half were warmup iterations. Convergence was assessed through analysing the MCMC traces, using the criteria of requiring all values to be below 1.05, and no reports of divergent transitions in the sampling [58]. This was achieved using default MCMC sampling parameters for the serotype analyses, but the strain analyses required a smaller step size (0.01), higher acceptance probability (0.99) and larger tree depth (20). The data and code used for the analyses are available as an R package from https://github.com/nickjcroucher/progressionEstimation/.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of parameters used in the models.

https://doi.org/10.1371/journal.pcbi.1009389.t001

Model evaluation and comparison

The ability of models to recover known parameter values from synthetic datasets was tested by simulation, as described in S1 Text. Models were compared using likelihoods, independent of priors, through approximate leave-one-out cross-validation with the loo package [62]. The marginal likelihoods of different models, given the data, were compared with Bayes factors calculated using the bridgesampling package [63]. The same number of iterations were used to calculate the logarithmic marginal likelihoods as were used to fit the models. Evaluation of our results against odds ratio calculations used the metafor package [64], as described previously [29]. Model outputs were analysed and plotted with the tidyverse and ggpubr packages [65,66].

S. pneumoniae serotype data

The matched carriage and IPD serotype datasets were combined from two meta-analyses [29,36]. Twelve of the reported studies were omitted due to lack of publicly-available data [31,67–72]; difficultly in defining independent cross-sectional carriage samples [73,74]; small sample sizes once stratified by age [75], or study design biased towards particular serotypes [76]. If a PCV were introduced during a study, then samples were divided into pre-PCV and post-PCV datasets. This resulted in 21 systematically-sampled and comprehensively serotyped paired asymptomatic carriage and disease samples (S1 Table). Of these studies, 20 included data on child carriage and child IPD, and five included data on child carriage and adult IPD.

Data that were only resolved to serogroup level, based on the currently-known serotypes, were omitted from the analysis, although the overall number of samples taken (η_i) was not adjusted. However, data on serotypes in historical studies that are now known to correspond to multiple serotypes (e.g. serotype 6A now being resolved into 6A and 6C [77]) were included unaltered. Additionally, isolates of serotypes 15B and 15C were combined into the single 15B/C category.

S. pneumoniae strain data

Analysis of variation of invasiveness between S. pneumoniae strains in children used population genomic data from South Africa [37], a mixture of population genomic and genotyped data from the USA [37], and genotyped data from Finland [74], Oxford [31] and Stockholm [78]. The equivalent analysis of strain invasiveness using isolates mainly from adult disease used data from Portugal [54]. Some of the genotyped datasets were not as thoroughly documented as the serotyped datasets, and therefore it was not appropriate to fit models that lacked a study-specific scale factor that reduced their reliance upon accurate carriage sample information (S2 Text).

Bayesian meta-analysis of S. pneumoniae serotype invasiveness

Twenty matched IPD case and nasopharyngeal carriage datasets were used to quantify S. pneumoniae serotype invasiveness in children (S1 Table). From these, 7,340 carriage isolates and 2,851 disease isolates were extracted, across 72 serotypes (S1 Fig). Multiple models of invasiveness were proposed to analyse these data (see Methods). The null model assumed there was no systematic difference between serotypes’ invasiveness across datasets. The type-specific model assumed each serotype had a characteristic invasiveness that was consistent across locations. The study-adjusted model assumed invasiveness did not vary across serotypes, but could differ between studies. Finally, the study-adjusted type-specific model allowed invasiveness to vary between both types and studies. All four approaches to estimating the expected numbers of disease cases and carriage isolates were fitted assuming the number of disease cases would follow a Poisson distribution, and allowing for overdispersion by fitting a negative binomial distribution to these values, resulting in eight models overall. These models were able to accurately recover progression rates (S2 Fig) and infer overdispersion (S3 Fig) from simulated data (S1 Text), with the true values used to generate the synthetic datasets lying within the 95% credibility intervals of the estimates from the fit of the corresponding model.

There was limited power for distinguishing the different model structures when the counts of both carriage and disease isolates were low for a given type in a study. Therefore paired observations from an individual study were only used for the model comparison if either the number of carried, or disease, isolates was at least five. This reduced the dataset to 7,048 carriage isolates and 2,617 disease isolates across 45 serotypes (S4 Fig). The eight models were each fitted to the data using two Hamiltonian MCMC chains, both run for 25,000 iterations. All chains appeared to have converged, based on the traces of the posterior probabilities (S5 Fig) and all values being below 1.001 (S6 Fig). Comparisons of the observed and predicted isolate counts showed all models could accurately reproduce the observations from carriage studies, with the exception of the null and study-adjusted Poisson models (Fig 1). However, only the most complex models, adjusting for both S. pneumoniae serotype and study, were able to accurately reproduce the disease case counts.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Comparison of observed and predicted counts of each serotype within each study of S. pneumoniae invasiveness in children.

The points show the observed value on the horizontal axis, and the median predicted value on the vertical axis. The error bars show the 95% credibility intervals. The points and bars are coloured by the study to which they correspond; note that England and Wales is abbreviated to “E&W” (S1 Table). The red dashed line shows the line of identity, corresponding to a perfect match between prediction and observation. The left column shows the correspondence for carriage data (values of c_i,j), and the right column shows the correspondence for disease isolates (values of d_i,j). Each row corresponds to a different model: (A) null Poisson model; (B) null negative binomial model; (C) type-specific Poisson model; (D) type-specific negative binomial model; (E) study-adjusted Poisson model; (F) study-adjusted negative binomial model; (G) study-adjusted type-specific Poisson model; (H) study-adjusted type-specific negative binomial model.

https://doi.org/10.1371/journal.pcbi.1009389.g001

The deviations between the observed and predicted isolate counts can be quantified using the distributions of absolute residuals (S7 Fig). These show the negative binomial models replicate the carriage values closely, relative to the equivalent Poisson models, with correspondingly greater deviation from the disease counts enabled by the overdispersion permitted by the negative binomial distribution’s precision parameter, ɸ. That the ɸ values were below one for the three simplest negative binomial models suggested there was substantial unmodelled variation (S8 Fig). However, for the study-adjusted type-specific model, the ɸ value rose above one, indicating less overdispersion. Correspondingly, the similarly low residuals for the Poisson and negative binomial versions of this model suggested adjusting for study and serotype enabled progression rates to be estimated robustly with either distribution (S7 Fig).

The success of the study-adjusted type-specific models in reproducing the observations may represent overfitting by the most complex models. Therefore formal model comparisons were undertaken using leave-one-out cross-validation (LOO-CV) [62]. Although this suggested the study-adjusted type-specific models were the most appropriate, the large number of model parameters meant the Pareto k statistic diagnostic values were too high for these comparisons to be reliable (S2 Table) [79]. This was still true if the models were compared using only the likelihoods calculated from the disease counts, which were more constrained than those calculated from the carriage data (S3 Table). Hence models were instead compared using Bayes factors calculated from bridge sampling [63]. This approach was validated by demonstrating it was able to accurately assign simulated data to the model under which it was generated (S4 Table). The Bayes factors identified the null Poisson model as the most poorly fitting, whereas the study-adjusted type-specific negative binomial model was the most strongly favoured by the data (S5 Table).

Identification of highly invasive non-vaccine serotypes

This study-adjusted type-specific negative binomial model was therefore applied to the full dataset. Both the MCMC traces and statistics indicated the model fit converged after 25,000 iterations (S9 Fig). This enabled the estimation of invasiveness for 72 serotypes (Fig 2 and S1 Dataset), including all 24 included in vaccine formulations that are currently licensed or under development [20–22]. The analysis found the serotypes associated with high invasiveness and narrow credibility intervals were the three serotypes added to PCV7 to generate PCV10 (1, 5 and 7F), which are also present in higher-valency PCVs, and serotype 12F, included in PCV20. Other serotypes not included in current vaccine designs, but likely to be highly invasive, were 9L, 19C, 24B, 24F, 25A, 27, 28A and 46. However, the credibility intervals of some of these estimates were substantial, resulting from small sample sizes.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Invasiveness estimates for S. pneumoniae serotypes from the study-adjusted type-specific negative binomial model applied to the full dataset of serotype studies from child carriage and disease.

Points represent the median estimate, and are coloured according to the vaccine formulation in which they are found; all higher valency formulations encompass the serotypes found in lower valency formulations, with the exception of 6A not being present in PPV23. The point shape represents the sample size, summed across disease and carriage isolates from all studies, on which the estimates are based. The error bars show the 95% credibility intervals. (A) Estimates for all 72 serotypes in the full dataset. (B) Estimates when distinguishing between vaccine-type serotypes in unvaccinated and vaccinated populations. Vaccine serotypes in vaccinated populations are denoted with the suffix “_pv”.

https://doi.org/10.1371/journal.pcbi.1009389.g002

The serotypes associated with low invasiveness and narrow credible intervals were 11A (included in PCV20), 6C, 21, 23A, 23B, 34, 35B and 35F. Many other serotypes (e.g. 10F, 11F, 19B, 33A, 35A, 36, 39, 42 and 45) were associated with low invasiveness values, but these estimates had wide credibility intervals. The higher uncertainty generally corresponded with lower sample sizes, although some of these serotypes had a total sample size above 10 (e.g. 10F, 35A and 42), but were nevertheless rare in IPD. The absolute estimates of invasiveness relate to the reference population for the model fit, which was the post-PCV7 sample from the Navajo nation [39], as this study contributed the largest number of isolates (S1 Fig). The model fits found the maximal progression rates were above 10⁻² IPD cases per carrier per year, with the minimum point estimates below 10⁻⁴ IPD cases per carrier per year. However, the incidence of IPD in the Navajo population is high, with excellent disease surveillance [80], and therefore other studies were expected to report lower disease cases per carriage episode. Correspondingly, the study-adjustment scale factors were below one for most other datasets (S10 Fig), suggesting many locations would expect to detect around 30–50% of the number of IPD cases per carriage episode for the same carriage serotype composition.

The fitted study-adjustment scale factors were all associated with similarly narrow credibility intervals, demonstrating there were enough shared types between samples to robustly infer these parameters. For four locations (Atlanta, Bogota, the Netherlands and France), both pre-PCV and post-PCV samples were included in the meta-analysis, which might be expected to have similar scale factors, if the host population and disease surveillance were consistent across these eras. However, this was not the case for studies from Bogota. The model’s assumption that serotype invasiveness is constant across studies may be violated by a reduction in vaccine-type invasiveness following PCV introduction, if vaccine-induced immunity inhibited progression from carriage to IPD. Therefore the study-adjusted type-specific negative binomial model was refitted, treating vaccine-type serotypes as different types pre- and post-PCV (including 6A as a PCV7 type, and 6C as a PCV13 type [81]; S11 Fig). The point estimates of invasiveness for most vaccine serotypes dropped following PCV introduction, with the evidence strongest for the invasiveness of serotypes 14 and 7F falling after the introduction of PCV7 and PCV10, respectively (Fig 2). This fit further improved the consistency of scale factors across the Netherlands and Atlanta, but did not resolve the discrepancy between samples from Bogota (S12 Fig). These data suggest PCVs can reduce the invasiveness of vaccine types before they are eliminated by herd immunity, and therefore distinguishing between pre- and post-PCV invasiveness may help standardisation across meta-analyses.

Comparison of Bayesian modelling with odds ratios

The outputs of the study-adjusted type-specific negative binomial model were compared to those from estimating pneumococcal invasiveness from the same dataset using odds ratios (Fig 3). Odds ratios were combined across studies using either random (Fig 3A) or fixed effects (Fig 3B) models, and the results disaggregated based on the total number of isolates, summed across carriage and disease, for each serotype. Across all sample sizes, there was a significant positive correlation between the two measures of invasiveness. The range of invasiveness values extended over two orders of magnitude for both methods. Hence the study-adjusted type-specific negative binomial model produced similar relative estimates of invasiveness to odds ratio analyses, while having the advantage of providing both overall, and location-specific, absolute estimates of these progression rates that can be used in quantitative analyses.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Comparison of the serotype invasiveness estimates from the study-adjusted type-specific negative binomial model with those from meta-analyses of the same dataset calculated as logarithmic odds ratios combined using (A) fixed effects models and (B) random effects models.

Points are coloured based on the vaccine formulations in which serotypes are found. The horizontal error bars are the 95% credibility intervals from the Bayesian model, and the vertical error bars show the 95% confidence intervals from the logarithmic odds ratio analyses. Each plot is divided based on the sample size associated with each serotype, calculated as the sum of isolates from carriage and disease across all studies. The correlation between the two estimates is summarised as Kendall’s τ on each panel, with an associated p value.

https://doi.org/10.1371/journal.pcbi.1009389.g003

Fixed effects analyses assume there is a single invasiveness value across all locations, whereas random effects analyses allow for variation in the underlying invasiveness estimate between studies. The identification of between-location heterogeneity typically suggests random effects models are more appropriate for analyses of pneumococcal invasiveness [29,36,37]. This variation may either reflect genuine differences in serotype behaviour, or arise from standardisation relative to a variable mix of serotypes between locations. By accounting for this variation, random effects models often calculate wider confidence intervals than fixed effect models (Fig 3). Additionally, random effects models are typically not recommended where there are fewer than five studies in which a serotype features [82], and therefore meta-analyses of odds ratios may be a mixture of fixed- and random-effects models [37]. By contrast, this single, consistent Bayesian model can be applied across serotypes, regardless of their distribution across studies. Furthermore, this framework explicitly accounts for the uncertainty associated with overdispersed data through the negative binomial precision parameter, ϕ (Fig 1). This variation is quantified separately from the uncertainty in the progression rate estimate.

Deviations between the methods in the rank order of invasiveness point estimates were most evident at small sample sizes. This was most notable for serotypes 33A and 36, both of which were associated with relatively high logarithmic odds ratios (above zero) and relatively low Bayesian progression rate estimates (below 10⁻⁴ cases per carrier per year). Both were observed only in carriage: five isolates across four studies for serotype 33A, and one isolate in each of three studies for serotype 36. Hence the Bayesian model estimates provide a more informative representation of the limited available data.

At the upper end of the scale, the Bayesian model found the most highly invasive serotype to be 46, with credibility intervals indicating this serotype is unlikely to be of intermediate invasiveness [30]. By contrast, the confidence intervals for the odds ratios calculated for serotype 46 were larger relative to the variation between serotypes, with those calculated using random effects spanning the full range of invasiveness point estimates across the species. This serotype was observed in two studies, but never isolated from carriage, and there is further circumstantial evidence that serotype 46 is likely to be highly invasive (see Discussion). Hence this Bayesian framework can use case-carrier studies to provide informative estimates of invasiveness even from small sample sizes.

Odds ratios were previously found to correlate with another measure of invasiveness, the “attack rate”, which in turn negatively correlated with carriage duration [18]. This could result from IPD being most common shortly after the acquisition of a novel serotype in the nasopharynx, which occurs more frequently for serotypes carried for only a short period in each host. The invasiveness values from child IPD data were compared with recent estimates of carriage duration from multi-state modelling of longitudinally-sampled carriage studies in Maela, Thailand [57,83]. Although the highest invasiveness values were associated with short carriage durations, likely as longer carriage duration is expected to result in higher cross-sectional carriage prevalences, there was no strong overall relationship (S13 Fig).

Differences in invasiveness between host age demographics

Although S. pneumoniae bacteria primarily circulate between children, they frequently cause disease in adults. Five datasets were identified in which adult disease cases could be matched with child nasopharyngeal carriage datasets, to estimate serotype invasiveness in adults. Overall, 3,756 carriage isolates and 3,041 disease isolates were extracted across 53 serotypes (S14 Fig). As with the child IPD samples, this dataset was filtered for observations where either the number of carriage or disease isolates was at least five, to improve the power for identifying the best fitting model. This left 3,704 carriage isolates and 2,969 disease isolates across 40 serotypes for this model selection analysis (S15 Fig). As with the child dataset, two 25,000 iteration Hamiltonian MCMCs were used to fit each of the eight models, which converged based on the traces of the posterior probability (S16 Fig) and values being below 1.05 (S17 Fig). Similar to the analysis of invasiveness in children, comparisons of the observed and predicted isolate counts again showed that only the models adjusting for both S. pneumoniae serotype and study were able to accurately reproduce the adult disease case counts (S18 Fig). Bayes factors concurred that the study-adjusted type-specific negative binomial model was the most likely, given the data (S6 Table).

This model was applied to the complete adult disease dataset. MCMC traces and statistics again indicated the model fit converged after 25,000 iterations (S19 Fig). As with the child analysis, the reference population for the model fit was the post-PCV7 sample from the Navajo nation [39]. The high incidence of IPD in this population is again evident in the adult population, as the other studies had lower study-adjustment scale factors (S20 Fig). Most locations expected to detect 5–50% of the number of IPD cases per carriage episode for the same carriage serotype composition.

This analysis therefore enabled the estimation of invasiveness of 53 serotypes in adults (Fig 4A and S2 Dataset), including the 24 in vaccine formulations that are currently licensed or under development. The analysis showed the serotypes associated with high invasiveness and narrow credibility intervals not in currently-used PCVs included serotypes 8 and 12F, included in the PCV20 formulation; serotype 20, included in PPV23, and the non-vaccine type serotype 13. Other serotypes not in current vaccine designs, but likely to be highly invasive, included serotypes 9L and 12B. However, the credibility intervals of these estimates were substantial, resulting from small overall sample sizes. The non-PCV serotypes associated with low invasiveness (upper 95% credibility interval estimates <10⁻³ disease cases per carrier per year) were serotypes 6C, 18B, 21, 23A, 23B, 35B and 35F.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Invasiveness estimates for S. pneumoniae serotypes from the study-adjusted type-specific negative binomial model applied to the full dataset of serotype studies from child carriage and adult disease.

Points represent the median estimates, and are coloured according to the vaccine formulation in which they are found; all higher valency formulations encompass the serotypes found in lower valency formulations, with the exception of 6A not being present in PPV23. The error bars show the 95% credibility intervals. (A) Estimates for all 53 serotypes in the full dataset. The point shape represents the sample size, summed across disease and carriage isolates from all studies, on which the estimates are based. (B) Scatterplot comparing the invasiveness estimates for serotypes across adults and children.

https://doi.org/10.1371/journal.pcbi.1009389.g004

Comparisons with the corresponding child invasiveness analysis (Fig 4B) showed many PCV serotypes had similar invasiveness estimates in both age groups [39], with PCV7 and PCV10 types having mid-range and high invasiveness estimates for both children and adults, respectively. The invasiveness of serotype 6C, which is affected by PCV13-induced immunity against the 6A component, was low in both children and adults. The non-PCV serotype 9L appears to be highly invasive in both populations, although this value was associated with considerable uncertainty. However, other non-PCV types showed a discrepancy between the two age groups. For example, serotypes 12B and 29 were estimated to be highly invasive in adults but not children, though the credible interval ranges were large for both.

Serotypes included in the PCV15 and PCV20 formulations appeared to have mid-range invasiveness in both age groups, except for the rarely invasive serotype 11A, and the consistently highly invasive serotypes 8 and 12F. This suggests the removal of these latter two serotypes through herd immunity from higher-valency PCV infant immunisation programmes would likely be beneficial to the adult population [84,85]. However, this would concomitantly reduce the effectiveness of the current PPV23 adult vaccine, due to its overlap in serotype coverage with PCVs meaning it would protect against a lower proportion of the post-vaccine S. pneumoniae population [29,86]. PPV23’s residual effect would be determined by the prevalence and adult invasiveness of the final non-PCV PPV23 serotypes. Both serotypes 2 and 20 were estimated to be highly invasive for both age groups, albeit with broad credible intervals for adults for the former, due to it only being included in one adult dataset. The remaining non-PCV PPV23 serotypes (9N and 17F) also exhibited elevated invasiveness compared to many non-vaccine types, suggesting that PPV23 may offer protection against a limited number of higher-risk serotypes following the implementation of PCV15 or PCV20 infant vaccination programmes.

Invasiveness varies between pneumococcal strains and serotypes

S. pneumoniae populations are genetically diverse, and can be divided into discrete strains [52,87]. Isolates of a particular serotype are often found across multiple strain backgrounds, and a single strain may be associated with multiple serotypes through switching [53,88]. To test how each of these characteristics affected pneumococcal progression rates, five different models of isolates’ invasiveness were fitted using the study-adjusted type-specific model framework (see Methods). These corresponded to the hypotheses that an isolate’s invasiveness was determined by its serotype; by its strain background; primarily by its serotype, but modified by strain background; primarily by its strain background, but modified by its serotype; and determined by the combination of its serotype and strain background. All five models were fitted using Poisson and negative binomial distributions.

These ten models were used to conduct a meta-analysis of six studies in which both serotype and strain background could be determined, for at least some isolates (S21 Fig). Three of the studies (post-PCV7 and post-PCV13 South Africa; post-PCV7 USA) were modified versions of a comparison primarily conducted using genomic data [37]. These were combined with pre-PCV genotyped studies from Oxford, UK [31]; Stockholm, Sweden [78], and Finland [74] (see S2 Text and S7 Table). For each dataset i, isolates were grouped by both their serotype j and strain k. As with the child serotype-only analysis, models were fitted to a subset of these samples in which the count of each serotype-strain combination in either carriage or disease in a study (c_i,j,k or d_i,j,k) was at least five. To focus on testing for within-category variation, only strains associated with five different serotypes, and serotypes associated with five different strains, were included in the analysis. The final dataset comprised 11 serotypes, 35 strains, and 46 serotype-strain combinations (S22 Fig). All models converged on stable set of parameter estimates, as inferred from the traces of logarithmic posterior probability values (S23 Fig) and distribution of values (S24 Fig). Few observed frequencies of serotype-strain combinations in either carriage or disease were inconsistent with the 95% credible intervals calculated from any of the fitted models (S25 Fig). There was little evidence of negative binomial distributions improving the fit of these models, and the precision parameters were correspondingly high, albeit with evidence of greater overdispersion when serotype information was omitted from the model fit (S26 Fig). The lowest absolute residuals were associated with the more complex models that accounted for both strain background and serotype (S27 Fig).

Correspondingly, comparisons with LOO-CV found such models that combined genetic and serotype information to be the best-performing (S8 and S9 Tables), although the Pareto k statistic diagnostic values were again too high for this comparison to be reliable. Comparisons using Bayes factors concurred, identifying the most likely model as that in which invasiveness was primarily determined by serotype, but modified by strain, with the count of disease isolates following a Poisson distribution (S10 Table). Some carriage sample sizes had to be approximated in this analysis (for the Oxford and Finland studies in particular; see S2 Text), but the study adjustment factors were robustly estimated and significantly differed across studies, suggesting the model was able to compensate for this aspect of the data (S28 Fig). To test whether these problems could have affected relative model likelihoods, the model comparison was repeated with a 100-fold change in the carriage sample size values for those studies in which the precise number was not reported. The study-adjustment factor meant the comparisons were relatively insensitive to these changes, and the same model was identified as the most likely, given the data (S11 Table). Hence model comparisons demonstrated that neither serotype nor strain background alone determines a genotype’s invasiveness.

Identification of invasive pneumococcal genotypes

The best-fitting type-determined strain-modified Poisson model could be reliably fitted to the full dataset (S29 Fig). The invasiveness estimates were plotted to analyse variation by serotype within strains (Fig 5 and S3 Dataset). Some strains expressed serotypes of uniformly low invasiveness (e.g. serotypes 11A, 15A and 20A within GPSC22), whereas others expressed serotypes with consistently higher invasiveness estimates (e.g. serotypes 4 and 19A within GPSC27). By contrast, 19A had an elevated invasiveness relative to 19F within GPSC1, and relative to serotype 15B/C in GPSC4. Similarly, serotypes 8 and 33F were substantially more invasive than 11A within GPSC3. Hence there was considerable variation in estimated invasiveness both between, and within, strains.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Estimates of the invasiveness associated with serotypes inferred from genotyped isolates using the study-adjusted type-determined strain-modified model.

Serotypes are arranged by the strains in which they are found. Data are only shown for strains found in multiple serotype-strain combinations, each represented by at least ten isolates across the studies with genotype information. Points represent the median estimate, and are coloured by the vaccine formulations in which the corresponding serotype is present. The error bars represent the 95% credible intervals. The shape of the point represents the sample size on which the estimate is based.

https://doi.org/10.1371/journal.pcbi.1009389.g005

Similarly, the factors by which strain backgrounds modified serotypes’ invasiveness were plotted, grouped by the serotypes with which they were associated (Fig 6). The meta-analysis identified significant evidence to support many of the observations of within-serotype invasiveness variation noted in the individual analyses (S12 Table). These included elevated invasiveness being associated with GPSCs 1, 24 and 41. Heterogeneity was most strongly evident within some of the paediatric serotypes [19], including 6A, 6B, 14, 19F, 19A, 23F and 23A. This may represent these common serotypes being distributed across a wider diversity of strain backgrounds. However, the estimates for other low invasiveness serotypes disseminated across multiple backgrounds, such as 15B/C and 35B, were more consistent. Furthermore, some highly-invasive serotypes (such as 1 and 12F) were similarly invasive in multiple strains. This implies some serotypes may have a stronger effect on an isolate’s invasiveness than others.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Estimates of the coefficient by which strains modify the invasiveness of their expressed serotype, inferred from genotyped isolates using the study-adjusted type-determined strain-modified model.

Strains are arranged by the serotypes with which they are associated. Data are only shown for serotypes found in multiple serotype-strain combinations, each represented by at least ten isolates across the included studies. Points represent the median estimate, and are coloured by the vaccine formulations in which the corresponding serotype is present. The error bars represent the 95% credibililty intervals. The shape of the point represents the sample size on which the estimate is based.

https://doi.org/10.1371/journal.pcbi.1009389.g006

The overall invasiveness estimate was plotted for each serotype-strain combination in the full dataset, including low frequency types, to identify high-risk genotypes that might emerge as common causes of IPD post-PCV (S30 Fig). The genotypes identified as most concerning in the near-term were typically those expressing serotypes targeted by higher-valency PCVs and regarded as highly invasive, such as 8 and 12F. Yet there were also examples of genotypes with high invasiveness point estimates expressing serotypes not included in any planned PCVs (e.g. 15A, 18B, 25F, 28A, 33A and 33D), albeit these values were inferred from small sample sizes. Most of these genotypes are likely to remain rare, else their invasiveness estimates may regress to the population-wide mean as more data emerge. Ongoing surveillance will be required to identify whether any represent a potentially problematic genotype post-PCV at an early stage in their spread.

An additional genotyped dataset from Portugal (S31 Fig), in which the majority of disease cases were isolated from adults [54], was analysed with the type-determined strain-modified Poisson model (S32 Fig). There were few instances of multiple serotypes being isolated from the same strain, but several serotypes were distributed across multiple genetic backgrounds (S33 Fig). This found evidence of heterogeneity in invasiveness within serotypes 3, 6A and 23F. The elevated invasiveness of serotype 3 was associated with the GPSC6 strain, similarly found to have high invasiveness when expressing multiple serotypes in children (Fig 6). For serotype 6A, this heterogeneity involved strains that did not appear in the child dataset (Fig 6), suggesting further within-serotype variation in invasiveness may emerge as studies include a greater diversity of host age groups.