Data sources

For many generalised epidemics, such as those in sub-Saharan Africa, the primary source of time series data on HIV prevalence is a system of sentinel surveillance sites that undertake anonymous HIV-prevalence testing of blood samples from pregnant women attending antenatal clinics (ANC).13 Use of ANC data has been justified based on the argument that prevalence among pregnant women serves as a reasonable proxy for prevalence in the general adult population, although some concerns have been noted about the representativeness of the sample and the possible effects of changing rates of contraceptive use over time.14 15 More recently, results from national population-based surveys have indicated that ANC data tend to overestimate prevalence in the general population.16 Since these surveys are only available for 1 or 2 years in most countries, the Reference Group model uses ANC data to inform the shape of prevalence over time and calibrates this curve to match the level of the national population-based survey estimates when available.5 As we are primarily concerned in the present study with characterising the shape of prevalence curves generated under alternative models, we omit this calibration step here in order to facilitate visual inspection of model output in relation to empirical ANC data. Our approach, however, generalises easily to a more detailed overall modelling approach that includes calibration to estimates from national surveys.

We fit the model separately to a total of 18 populations, distinguishing urban from rural regions in nine countries: Angola, Burundi, Ethiopia, Gabon, Kenya, Lesotho, Rwanda, Uganda and Zimbabwe. Rwanda and Uganda offer two important examples in which the Reference Group model fails to capture recent trends in prevalence.

Proposed model

Our proposed alternative to the current EPP model explicitly represents changes in infection risks over time to better approximate patterns of HIV spread in generalised epidemics. We simplify the population structure of the Reference Group model by collapsing the ‘susceptible, at-risk’ and ‘susceptible, not at-risk’ compartments into one uninfected (and at-risk) sub-population. In line with the modified structure, we remove two corresponding parameters from the Reference Group model — f₀, which governs the initial split between at-risk and not-at-risk entrants, and φ, which modulates this split in mature epidemics. In assessing the performance of the spline-based approach, we limit our analysis to scenarios that do not include ART, although we note again that our approach generalises easily to a more detailed specification including ART.

The model can be described by three differential equations, following the general notation for the Reference Group model used in EPP 200917:dY/dt=V(t)−∫0tV(x)·g1(t−x)dx−μY(t)dU/dt=∫0tV(x)·g1(t−x)dx−∫0tU(x)·g2(t−x)dx−μU(t)dZ/dt=E(t)−V(t)−μZ(t)whereV(t)=r(t)·(Y(t)+U(t))·Z(t)/N(t)

The key difference in the parameterisation of this model compared with that in the current Reference Group model is the use of a time-varying function for infection risks, r(t), modelled using splines, rather than assuming constant r. A number of other key elements of the model are consistent with the current Reference Group formulation. The model distinguishes early-stage from late-stage infection (eg, those having CD4 counts above or below 200 cells/μl) in order to accommodate the modelling of treatment eligibility. Y represents the number of HIV-positive individuals at the earlier stage, U is the number of HIV-positive individuals at the later stage, and Z is the number of uninfected individuals (Y+U+Z=N). E(t) is the number of new adults entering the population at time t and is a function of HIV prevalence 15 years earlier, which accounts for differential fertility and child survival by maternal HIV status. The functions g₁(x) and g₂(x) determine HIV disease progression, from the early to the more advanced stage, and from the advanced stage to AIDS-related death, respectively. Both functions are parameterised as Weibull densities, with shape and location parameters α=1.69 and β=10.40 for g₁(x) and α=1.00 and β=2.52 for g₂(x), yielding median residence times of 8.4 and 1.8 years, respectively.

We make another important modification to the Reference Group model by changing the way in which the model captures the initial phase of an epidemic. Fitting the Reference Group model currently requires estimation of the year (t₀) in which a fixed seed value of Y_t0=0.0025 initiates the epidemic. In our proposed specification, an epidemic is initiated in 1975 with a small pulse of infected individuals Y₀, which is a parameter to be estimated with a broad prior uniform distribution spanning the range 10⁻¹³ to 0.0025. Given the large uncertainty about when HIV first appeared in human populations,18 and the inevitable artifice of simulating the introduction and initial spread of HIV as a discrete event, this new approach avoids having to estimate the ‘start year' of the epidemic, while acknowledging substantial uncertainty about the early dissemination of HIV through the population. It also ensures a constant number and position of splines for the r function and improves efficiency, as the fitting algorithm no longer needs to be re-run for each potential value of t₀.

There are several possible choices of spline functions for modelling infection risks over time, r(t).19 20 We propose to use third degree B-splines with a difference penalty to enforce smoothness (also referred to as ‘P-splines’).21 B-splines are easy to generate (the basis functions consist of polynomial pieces) and have stable boundaries. With this formulation, a relatively large number of basis functions can be set at equal intervals, and a difference penalty on the coefficients of adjacent basis functions determines the amount of smoothing.21 Here we use a second-degree difference penalty, which enforces smoothness in the second derivative of the spline function (ie, smooth changes in the slope of r over time in the epidemic model). Figure 1 shows an example of how one trajectory for r(t) (thick solid line) is generated from seven B-spline basis functions (dotted lines), as used in our analysis. The function r(t) is constructed as a linear combination of the basis functions, with the spline coefficients (denoted β₁ through β₇) serving as weights.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1

Example of a spline-based curve for infection risks over time (solid line) using third degree B-splines comprised of seven evenly-spaced basis functions (dashed curves labelled 1 through 7). The spline curve is a linear combination of the basis functions and their coefficients (labelled β₁ through β₇).

In line with the current strategy for estimating model parameters and characterising uncertainty, we adopt a Bayesian approach to fitting the model to surveillance data,22 specifically using a Bayesian analogue to B-splines with difference penalties.23 For the second-degree difference penalty used here, the coefficient for basis function i is expected to follow the trend of the previous two coefficients, and its prior can be expressed as: β_i=β_{i -1}+(β_{i -1}–β_{i -2})+u_i. The priors for β₁ and β₂ are assumed to be proportional to a constant. Note that if all of the β coefficients are equal, the resulting ‘curve’ will be a horizontal line. On the other hand, if a spline coefficient differs greatly from its neighbouring coefficients, we expect the curve to change in value rather quickly over the range of that basis function. By letting u_i ∼ normal(0,τ²), the amount of smoothness is determined by the variance parameter τ², as smaller values for τ² will produce smaller values for u_i, which in turn imply a smoother function. The parameter τ² is estimated during model fitting, and we use the common uninformative prior τ² ∼ inverse-gamma(0.001,0.001).24 In preliminary analyses we tested between four and 10 evenly-spaced basis functions and found that seven basis functions were sufficient to obtain adequate flexibility for r(t) to match a wide variety of different empirical prevalence trends. When fitting the model, we add a further constraint that r must always be positive.

Short-term projections

While the primary objective of the model described here is to estimate past epidemic trends for reporting on the current status of HIV/AIDS epidemics, one of the required functions of the model is to make short-term projections beyond the last observed data point in settings that lack data for the most recent reporting year(s). There are a range of methodological challenges in extrapolating time series beyond the last observed data point, including both general challenges and certain challenges specific to spline models.25 26 In order to address the challenge of extrapolation, we incorporate a prior probability distribution for future values of r derived from observed prevalence levels in the most recent years of data. The specification of the prior draws on the mathematical theory of infectious disease dynamics.27 Briefly, if we let I be the proportion infected and S be the proportion uninfected, with infection rate r and mortality rate μ, the change in the infected population is given by dI/dt=rIS–μI. At later stages in an epidemic, we expect prevalence to approach equilibrium, all else being equal. Under this formulation, dI/dt=0 when r=μ/S, which we approximate as r≈1/(1– prevalence)×1/duration. (This is an approximation given a non-constant mortality rate for HIV/AIDS.) With the mortality assumptions in EPP 2009, mean duration with HIV/AIDS is approximately 11.5 years. We implement this prior for r(t) by assuming that future values of r are normally distributed, with mean r(t)=1/(1–Y(t^*)/N(t^*))×1/11.5, where Y(t^*)/N(t^*) is modelled prevalence in the last year of data. Variance of the distribution is set equal to the mean of the squared deviations between r(t) and 1/(1–Y(t)/N(t))×1/11.5 for the last 10 data-supported years of the projection. The prior is included in the model fitting process and serves to guide the spline curve (which is generated as a continuous function for the entire estimation and projection period) beyond the last observed data point.

We note that this prior is intended only to ensure that most curves for r reside within a reasonable range for later-stage, fully generalised epidemics, and not necessarily to capture long-term dynamics in future predictions. By ‘reasonable range’, we suggest, for example, that an annual value of r=1, which implies an average of more than 11 secondary infections over the lifetime of an infected person, would be unreasonably high for a late-stage epidemic as illustrated by Uganda in 2010, given that the average lifetime number of sexual partners has been estimated at around 2–3 for women and 10 for men across several African studies,28 and that most partnerships do not result in infection (the average probability of transmission in a single sexual act has been estimated to be as low as 1 per 1000).29 A recent modelling study based on long-term follow-up of a cohort in Rakai, Uganda estimated that approximately two new infections would be generated by each infected individual based on random mixing in sexual partnerships,30 which underscores that values for r should be much less than 1 and within the range implied by the prior used here for future values of r.

Model fitting

Following the work of Bao and Raftery, we use incremental mixture importance sampling (IMIS) for model fitting, which is a Bayesian algorithm for simulating posterior distributions.31 We employ the same likelihood as that used currently in EPP, which assumes that prevalence is normally distributed on the probit scale, with random effects for ANC sites.32 The spline-based model was programmed in the software package R, with some functions written in C to increase computational speed, and the model differential equations were approximated as difference equations in quarter-year increments. For comparison with the existing Reference Group model, results were obtained by running the uncertainty analysis option in EPP 2009. For the spline-based model, we added the additional constraints that r <20 and prevalence <95% to increase the efficiency of the IMIS algorithm. These constraints allowed us to avoid calculating the likelihood for proposed parameter sets that led to epidemiologically implausible epidemic trajectories. The choice of r <20 was a very weak constraint; we did not observe values of 10 or greater in the posterior distributions of r for any country. We started the IMIS algorithm with 900 000 initial draws, with 900 draws for each importance sampling iteration.

Model evaluation

We compared the alternative models in terms of their observed fit to surveillance data from all sites. The performance of the models was also evaluated based on more formal criteria. To compare the relative fit of the spline-based and Reference Group models to ANC data, we computed Bayes factors33 along with density plots of the marginal likelihoods for both models. We further evaluated the spline-based model with posterior predictive checks of in-sample fit by comparing posterior predicted prevalence to observed prevalence for each site-year observation. We performed these posterior predictive checks after constructing case-deleted posteriors for each ANC site's random effect via sample reweighting,34 and we obtained predicted random effects for ANC prevalence via rejection sampling.32 To assess the model's performance in out-of-sample projection, we truncated the last 3 years of data, re-fit the model and determined if the 95% credible intervals for prevalence covered the median projection from the model fit to the full ANC time series. We tested the robustness of the IMIS algorithm by re-running our standard projections with a different random seed and comparing the posteriors for r, incidence and prevalence. We also investigated the sensitivity of the results to the prior on τ², and we visually compared posterior distributions of each model parameter to the distribution of initial draws for the IMIS algorithm.