A Bayesian approach to estimate changes in condom use from limited HIV prevalence data

In this section we present the part of our modelling framework that aims to capture the dynamics of HIV transmission in an environment with features as in the Mysore district in Southern India. The model is based on that of Pickles et al. (2010) and can be described as follows: a stable but open population is considered, consisting of FSWs and their clients. In line with Vickerman et al. (2010), we assume that the remaining members of the population (eg. companions of clients or sex workers), that are not involved directly in sex work, have little influence on the dynamics of the epidemic and can therefore be ignored. The model structure allows for two different groups of FSWs, i.e. high-risk ($H$) and low-risk ($L$), that may have different numbers of clients ($C$). Each individual in these three groups can either be susceptible to HIV infection, infectious, or retired either due to death or ceasing commercial sexual activity. Moreover, individuals that either decease or stop being involved in commercial sex are assumed to be replaced by susceptible ones, such that the size of the risk population remains constant. A graphical representation of the model is given in Figure 1 via a flow-diagram for each of these three population groups.

A key component of the epidemic model is the force of infection $\beta(t)$ which is assumed to be the product of the number of clients multiplied by the probability of HIV transmission in at least one out of $n$ sexual acts per client encounter. A sexual act can lead to HIV transmission when either no condom is being used, occurring with probability $1-CU(t)$, or when a condom is not being used effectively, occurring with probability $CU(t)(1-e)$. Hence, the combined probability of a sexual act that can lead to HIV transmission, also referred to as risky sexual act, simplifies to $1-eCU(t)$. The probability of HIV transmission from the sex worker to the client (or from the client to the sex worker) in a single risky sexual act is given by $p_{S}$ (or $p_{C}$), thus providing the expression of the force of infection; e.g. for high risk FSWs we get $C_{H}[1-\{1-p_{S}\times(1-eCU(t))\}^{n}]$, where $C_{H}$ is the number of clients for high risk FSWs:

where the following notations are used:

• •

  $C_{H}$ and $C_{L}$: number of clients of FSWs per month, which differs for high risk and low risk FSW

• •

  $n$: number of acts per client encounter

• •

  $p_{S}$ or $p_{C}$: probability of HIV transmission from client to sex worker or sex worker to client respectively during an unprotected sex act

• •

  $e$: efficacy of condoms in protecting against transmission of HIV per sex act

• •

  $CU(t)$: proportion of FSWs’ sex acts that are protected by condoms, that we allow to vary in time (it is the main parameter that we want to estimate)

Additionally, the transmission dynamics of HIV will depend on the following lengths of time:

• •

  $\mu_{S}^{-1}$ or $\mu_{C}^{-1}$: average length of sexual activity as a sex worker / client

• •

  $\alpha^{-1}$: average life expectancy with HIV

Information on the above quantities is available from various serial cross-sectional bio-behavioural surveys (IBBAs) in Mysore district in southern India (Ramesh, et al., 2008). However, there is still uncertainty around these biological and behavioural parameters, which is reflected on the estimates of CU by using a Bayesian approach (De Angelis et al., 1998).

The model can also be defined via the following set of differential equations:

In the equations above, the state of the epidemic over time is described by the number of susceptibles among clients ($S_{C}$), among the low-risk female sex workers ($S_{L}$) and among the high-risk female sex workers ($S_{H}$), as well as the corresponding numbers of infected individuals ($I_{C}$, $I_{H}$ and $I_{L}$). Three types of constant parameters are involved: the initial prevalence among the different groups of interest in 1985 ($\theta_{i.c.}=\{S_{H}(1985),I_{H}(1985),S_{L}(1985),I_{L}(1985),S_{C}(1985),I_{C}(1985)$), time-invariant parameters describing the biological and behavioural determinants of HIV transmission ($\theta_{tr.}=\{\mu_{S},\mu_{C},\alpha,C_{H},C_{L},p_{S},p_{C},n,e$}) and the parameters that play a role in the trajectory priors for CU between 1985 and 2010 ($\theta_{CU}$). All three components $\theta_{i.c.},\theta_{tr.}$, and $\theta_{CU}$ are integrated into a global vector of constant parameters, denoted by $\theta$. Under this notation, the trajectory $X_{1:n}$ of the space vector $X(t)=\{S_{H}(t),I_{H}(t),S_{L}(t),$ $I_{L}(t),S_{C}(t),I_{C}(t)\}$ is defined as a deterministic function of $\theta$ and CU ($X_{1:n}=f(\theta,CU)$) through an HIV transmission model, and is compared with the available observations, denoted by $y_{1:n}$. Note that the function $f(.)$ is not available in closed form but can be obtained given the trajectory of $CU$ by solving the above ordinary differential equations (ODE). More specifically, we introduce a time discretisation with equidistant points of time step $\delta$ resulting in a discretised skeleton of CU denoted $CU^{discr}=\{CU(t_{0}),CU(t_{0}+\delta),CU(t_{0}+2\delta),..,CU(t_{n})\}$. The partition of the CU trajectory can be made arbitrarily fine by the user-specified parameter $\delta$ to limit the approximation error induced by the time discretisation.

Assigning a model for the observation error provides the likelihood of the observation $y_{1:n}$ conditional on the CU trajectory $p(y_{1:n}|\theta,CU)$. In this paper we use a binomial distribution, considering that prevalence estimates are derived from a random sample of 425 FSWs or clients in the Mysore district. We denote $h_{i}(X)$ the value of the prevalence estimated through the HIV transmission model at time $i$: $h_{i}\{X(t_{i})\}=\{I_{H}(t_{i})+I_{L}(t_{i})\}/2$ when prevalence among FSWs is observed, $h_{i}\{X(t_{i})\}=I_{C}(t_{i})$ if prevalence among clients is observed. With these notations, the likelihood function is defined from the following:

Overall, the model appearing in Figure 1 and Equations 1 and 2 is a simplified version of the one in (Pickles, et al., 2010). This was done mainly for parsimony reasons; models of increased complexity can be used provided that there is adequate information on their parameters. More details on informative priors are provided in Section 2.3.

From a modelling perspective, the main contribution of this paper lies in the introduction of various formulations for the evolution of the CU, aiming to explore its space of trajectories rather than evaluate a limited set of scenarios. The proposed trajectory priors are motivated by two different objectives. First we assign a trajectory prior aiming to impose little information to its shape. We proceed with a Brownian motion for CU, transformed to take values in the real line, although integrals of Brownian motion can also be used. Loosely speaking these models can be linked with a smoothing splines approach (Wahba, 1990). The second objective aims for a slightly more informative formulation regarding the evolution of CU via sigmoid-shaped growth curves. Growth curves provide natural models for smoothly growing quantities in different contexts such as biology (Zwietering et al., 1990), marketing (Lessne and Hanumara, 1988) and epidemiology (Omran, 1971), while they were also mentioned in Pickles et al. (2010). Moreover, the posterior of the CU trajectory, based on the model with Brownian motion as a prior, revealed a sigmoid shaped growth. Hence, the second trajectory prior, denoted by dBR, is based on the generalised Bertalanffy-Richards model; see for example (Garcia, 1983; Yuancai et al., 1997). In order to enrich this context and address estimation issues that can be encountered with the dBR (Lei and Zhang, 2004), we also consider an alternative empirical sigmoid curve (dSigm). In what follows, we denote with $x$ the latent process that drives the CU trajectory, which in turn provides the link with the prevalence observations through the model in Section 2.1.

The parameters contained in $\theta_{i.c.}$ and $\theta_{tr.}$ cannot be identified from the prevalence observations only, so we assign informative priors on them. These are summarised in Table 1 and are similar with the priors used in Pickles et al. (2010). They were either based on previous literature regarding general quantities as transmission probability for unprotected acts, or life expectancy with HIV, whereas the ones concerning quantities that are more sociologically and geographically specific were estimated from cross-sectional individual-based surveys (IBBAs) in Mysore.

The parameter vector $\theta$ includes an additional component, $\theta_{CU}$, that contains the parameters for different models describing the CU trajectories, $\theta=\{\theta_{i.c.},\theta_{tr.},\theta_{CU}\}$. Although there is some information in the data for $\theta_{CU}$, the posterior will depend on the prior to a large extent. As mentioned earlier, there is very little information on the volatility parameter of the BM formulation. Throughout this paper we used a Uniform prior between $0$ and $2$ $months^{-1}$. As explained in more detail in Section 3.1, the parameter of main interest in this study is the quantity $\Delta CU=CU(2009)-CU(2003)$. Simulations suggest that, if we combine the BM approach with a $Unif(0,1)$ prior for $CU_{0}$ (CU in 1985), this results in a symmetric prior on $\Delta CU$ that is centered around 0 with 2.5% and 97.5% points at $\pm 0.6$ respectively. We considered it as a reasonably vague prior for $\Delta CU$ and evaluated the performance of the resulting model via the simulation experiments of Section 3. More diffuse priors can also be used by setting a larger value for the upper limit of the Uniform prior for $\sigma$. Regarding the parameters of the sigmoid curves, we used vague priors that are also shown in Table 1.

The joint posterior distribution can be obtained up to proportionality through the HIV transmission model of Section 2.1, which links the prevalence observations with the CU trajectories, the trajectory priors of Section 2.2 and the remaining priors of Section 2.3. For the dBR and dSigm trajectory priors, it can be put in a non-linear regression framework, with the non-linear function being the solution of the ODE, and can therefore be implemented with standard software such as WinBUGS through WBDiff (Lunn, 2004). However this is not possible for the BM case where more involved techniques are required. Since the posterior probability density function is intractable, a data augmentation scheme can be utilised. This inference problem poses some challenges due to the high dimension of the discretised representation of $CU_{1:n}$ and its strong correlation with the vector of constant parameters, $\theta$. This correlation imposes problems to Gibbs schemes on $\theta$ and $CU_{1:n}$, leading to extremely poor mixing and convergence properties. The Particle MCMC algorithm (PMCMC, see Andrieu et al. (2010)) algorithm offers a solution by updating the two components jointly, thus reducing the problem to a small-dimensional MCMC on $\theta$. Implementation is based on the estimates of the likelihood $\hat{p}(y|\theta)$ that are provided by a particle filter.

The particle filter and the PMCMC are described in Algorithms 1 and 2, in terms of the quantities introduced in the previous sections. More details about this algorithm and its practical implementation can be found in Dureau et al. (2012) (through an application in a similar context) and Andrieu et al. (2010).

Regarding the choice of $Q(.)$ in Algorithm 2, we use a random walk Metropolis-Hastings algorithm in a transformed parameter space (log or logit) to ensure positivity. Each iteration of the MCMC algorithm requires an execution of the particle filter, which induces substantial computational cost if the importance sampling covariance matrix $\Sigma$ is ill-adapted. Adaptive approaches (Roberts and Rosenthal, 2009) can be used to tune $\Sigma$ but they require lengthy explorations of the target space. We propose to speed up this process by pre-exploration of a proxy posterior density $p^{EKF}(\theta|y)$ relying on a Gaussian approximation of the dynamic system and the Extended Kalman filter methodology (Dureau et al., 2012). A simple bootstrap version of the particle filter is used as it is not straightforward to consider data-driven transition proposals given complex observation regime of our model. Note that given the short length of the observed time series, simpler alternatives to PMCMC may perform reasonably well. For example, the resampling step can be omitted and the particle filter output can be used to approximate $p(CU_{1:n}|y_{1:n},\theta)$. In our application however, it turns out that additional particles are needed for this approach, thus not offering a great reduction to the computational cost when compared to PMCMC. We therefore suggest the use of PMCMC as a robust computational tool that still does not requires a large amount of time in applications of this type. In order to meet with the computational requirements of the ensemble simulations at stake in the present article, and to facilitate future applications of the methodology we are proposing, the calculations were made using the SSM inference package. This package, targeted to time series analysis via state space models, generates executables for likelihood optimisation and Bayesian inference algorithms in parallelisable C for any stochastic compartmental model expressed in a high-level modelling grammar (see Dureau et al. (2014) for more information).

By fitting each of the previously introduced models we obtain samples from the marginal posterior density $p^{meth}(CU(t)|y)$ ($meth\in\{dBR,dSigm,BM\}$). However, our interest mainly lies in the amplitude of the shift in CU between 2003 and April 2009 measuring the estimated increase in CU during the intervention, henceforth denoted by $\Delta CU$. The posterior draws of CU trajectories can be transformed to provide samples from the posterior of this parameter of interest. The samples can then be used to form an estimator $\hat{\Delta}CU^{meth}$ of $\Delta CU$ such as the posterior median of $\hat{p}^{meth}(CU(t)|y)$. In what follows we explore the frequentist properties of this estimator derived from each of the trajectory priors.

It may also be of interest to assess the estimating capabilities, given the limited amount of data, for the hyperparameters of the various CU priors ($CU_{0}$, $\eta$, $r$, $m$, $t_{in}$ and $\sigma$). As it turns out there is information for some of them ($CU_{0}$, $\eta$, $t_{in}$), whereas some others are hard to estimate and are determined mostly by their prior ($r$, $m$ and $\sigma$). Nevertheless, from a subject matter point of view, interest lies mainly on $\Delta CU$, whereas the remaining quantities (in CU priors) can be regarded as nuisance parameters. Another appealing feature of $\Delta CU$ is that it appears in all models and therefore provides an omnibus quantity for comparison. Hence, inference properties of these parameters ($CU_{0}$, $\eta$, $r$, $m$, $t_{in}$ and $\sigma$) are only studied indirectly through inference properties of $\Delta CU$.

The performance of each estimator $\hat{\Delta}CU^{meth}$ in estimating $\Delta CU$ is evaluated from the following criteria (where $L=1000$, the number of simulations):

In addition to the quantities above we are also interested in assessing the discriminative ability of each model in detecting increases in CU. Focus is given to increases in CU that are at least as high as a pre-specified threshold T, which reflects the minimum practical increase. When analysing the data a researcher may decide that CU did increase more than T if the value of the estimator $\hat{\Delta}CU^{meth}$ is higher than a user-specified threshold $t$. Each decision mechanism may lead to different types of error and is therefore associated with a particular sensitivity and specificity. More specifically we can define the true and false positives in the following way

• •

  Sensitivity (true positives rate) for $t$: $\frac{\sharp(\hat{\Delta}CU^{meth}>t,\Delta CU>T)}{\sharp(\Delta CU>T)}$

• •

  Specificity (1 - false positives rate) for $t$: $\frac{\sharp(\hat{\Delta}CU^{meth}<t,\Delta CU<T)}{\sharp(\Delta CU<T)}$

We proceed by first reporting sensitivities and specificities corresponding to the case of $t=T$. This corresponds to saying that $\Delta CU$ is higher than $T$ if its estimator is higher than $T$. We then use a range of different $t$’s and obtain the sensitivity-specificity pair that corresponds to each of them. A lower detection threshold $t$ will increase the sensitivity of the method, but it also increases the risk for false positives, and vice versa. These pairs are combined to form the Receiver Operating Characteristics ROC curve by plotting sensitivity versus 1-specificity. The area under the ROC curve (AUC) provides an overall measure of discriminatory power as it reflects the probability of correctly classifying a randomly chosen positive instance as higher than a randomly chosen negative one (Fawcett, 2006). For example, an AUC value of 50% indicates no power (i.e random choice) This detailed procedure is repeated to assess the ability to detect two different levels of increase in CU, with $T$ set to $20\%$, $30\%$ and $40\%$ respectively.

The performance of the estimators derived from the different models is measured using a set of simulated experiments where CU trajectories are sampled from a given growth curve model ($dSigm$), and parameters from $\theta_{i.c.}$ and $\theta_{tr.}$ are sampled following their prior distributions. To maximise the utility of this test procedure for future application of this methods to help evaluate Avahan in different districts (manuscript in preparation), only plausible and realistic CU trajectories are considered: cases with prevalence in 2010 between 2% and 40% and with CU shifts that occurred after 1995. Furthermore, the test trajectories have been sampled so that $\Delta CU$ regularly spans the $[0;0.9]$ interval. Alternatively, a second set of simulations has been generated based on step-wise CU trajectories (constant CU up to $t_{in}$, and then constant CU after $t_{in}$) to explore the impact of model misspecification and to assess the sensitivity of results to the underlying model used for simulations ($dSigm$).

For each of these experiments, an epidemic is simulated to provide observations $(y_{i}^{sim})$ replicating the observation framework applied in Mysore: three prevalence estimates among female sex workers and one among clients, concentrated during the period of the intervention. From these observations, the MCMC algorithm is applied to each method $meth$ to sample from $p(CU^{meth}(t_{1}:t_{n})|y^{sim}_{1:n})$. Then, given the posterior CU samples the estimators $\hat{\Delta}CU^{meth}$ can be computed, and compared to their true counterparts $\Delta CU$ by calculating the measures of performance of the previous subsection.

This section contains the results of the simulation experiments that were conducted as described in the section 3.3. In all cases each model is fit to each of the $1000$ simulated datasets and the estimator is formed via the posterior median at the points of interest. Various frequentist properties of these estimators are then assessed such as bias, the standard deviation and the MSE. Moreover, in Table 3, we examine the ability of each model to classify the simulated instances of the $\Delta CU$ parameter, and assess the risk of overstating versus understating. In other words, we assess the ability of model-driven estimators to address questions such as ’was the shift in CU during the intervention over 0.2, 0.3, or 0.4?’, by looking at their sensitivity and specificity as well as the resulting AUC.

We begin with results corresponding to cases where the data were simulated from the dSigm model. Table 2 reports bias, standard deviation and MSE for each model-driven estimator. It suggests that all models underestimate $\Delta CU$. More precisely, the dBR model tends to strongly understate the shift in amplitude, by 0.33 in average. The biases of the dSigm model is only slightly smaller (-0.31), and optimal results are obtained with the Brownian Motion model (-0.26). Similarly, in terms of MSE, the performance of the BM model is better. Figure 2 plots the bias against the true value of $\Delta CU$ for each model and reveals an increasing association. Another interesting note is that the ranking of the different models is consistent across the entire range of $\Delta CU$ values. If, for example, the true shift is between 40% and 50%, it is on average underestimated by 0.21 with the best method (BM) and more than 0.30 points with the dBR and dSigm methods.

The bias in estimating $\Delta CU$ under each of the CU priors can be attributed to a large extent to the prior implied by each formulation on $\Delta CU$. As mentioned in Section 2.3 the BM approach results in a symmetric prior on $\Delta CU$ that is centered around 0 with 2.5% and 97.5% points at $\pm 0.6$ respectively. The posterior median is therefore pulled towards zero, thus resulting in conservative estimates. The amount of shrinkage depends on the upper limit of the Uniform prior on $\sigma$. The corresponding priors under the dBR and dSigm formulations result in priors for $\Delta CU$ that put more mass around zero and this heavily depends on the priors placed on $k$ and $m$. It would perhaps be interesting in the future to consider alternative priors for $k$ and $m$ to reduce the bias in $\Delta CU$. Regarding the dBM model, the Uniform prior assigned on $\sigma$ appears to work reasonably well.

Table 3 suggests that all estimators based on the median of the posterior density of $p(\Delta CU|y_{1:n})$ have good distinguishing power: the AUC is between 0.87 and 0.93 in all cases. In line with the results of Table 2, the estimates provided by the BM model achieve better sensitivity (53%, 29% and 14%) than the other models (between 0% and 29%). The strong negative bias observed in Table 2 is also reflected by an optimal 100% specificity of the estimators. The performance, in terms of sensitivity, decreases as the level of $\Delta CU$ increases. Figure 3 also depicts the ROC curve for the dBM model, obtained as described in section 3.3.

In an additional simulation experiment we examine the frequentist properties of the models in the case where underlying CU trajectories were generated from a step function. As shown in Table 4 a very similar picture is obtained with the BM model possessing the smallest bias and MSE, while dSigm being slightly better than dBR.

Overall, the results presented in these tables provide an informative qualitative assessment for the ability of the different models to capture $\Delta CU$ from limited prevalence data on an important and diverse set of likely scenarios (1000 experiments with true amplitude in CU shift spanning between 0 and 0.9). First of all MSE and AUC figures suggest that although the number of prevalence observations is low and some elements of the transmission process are uncertain, it is still possible to extract information on our time varying parameter and provide estimates of the amplitude of the shift in CU during the intervention. Furthermore, there seems to be a possibility to control the risk of overstating these quantities by analysing the outputs of the three models that offer different levels of compromise between sensitivity and specificity. Thus, although the procedure may fail to identify some shifts in CU, there seems to be some reassurance that once a shift is detected it will most likely be true. In other words the obtained model-driven estimators seem to be conservative, yet robust on detecting CU trends.

The two models with the higher overall performance, BM and dSigm, are quite different in nature: the dSigm trajectories are smooth, whereas under the Brownian motion prior they are non-differentiable. Hence, the choice between the two models can also be based on prior beliefs of the researcher regarding the smoothness of the CU trajectories.

Mysore is one of the districts targeted by the Avahan intervention, and Avahan was the first HIV prevention intervention in this region. Four HIV prevalence estimates have been obtained between 2003 and 2009, three among female sex workers, and one among clients. Results from the inference procedure using a Brownian motion model are shown in Figure 4, suggesting a strong impact of the intervention. The purpose of this paper was to determine what level of increase of CU between 2003 and 2009 can be inferred while controlling the risk of overstating it. As it was shown in section 4.1, dSigm models could provide a good alternative to the BM formulation. Hence, we also present here results obtained with this model for the Mysore dataset (see Figure 4). An interesting point to note here is the lower posterior uncertainty around the CU trajectory in the case of the dSigm model. This is in line with our expectations given the fact that this prior imposes a particular growth structure and therefore introduces additional information. Loosely speaking this may be viewed as a bias-variance tradeoff. The choice of each CU formulation is likely to affect the inference procedure; therefore it is useful to explore its implications, for example via the simulation experiments of the previous section. Table 5 shows the estimates of $\Delta CU$ for each of the three presented models. The results indicate a positive increase in all cases. In particular, for the BM and dSigm models the corresponding posterior medians are 0.54 and 0.49 while the 95% credible intervals are [0.06;0.95] and [0.09;0.97] respectively.

A stronger conclusion regarding a lower bound for the CU shift between 2003 and 2009 can be made by comparing the posteriors medians to the results of Table 3. If the underlying set of simulations is to be considered realistic, an argument in favour of a CU increase being at least 0.4 can be made. Since the posterior medians are more than 0.4 under both BM and dSigm models (.54 and .49 respectively), Table 3 suggests that a statement for $\Delta CU>0.4$ will be correct with probability given by the specificity of each model (100% in both cases). While being more informative than the credible intervals obtained directly from the posterior densities (over 0.06), these numbers are heavily dependent on the assumption that the simulations of Section 3 provided an adequate approximation of the reality. Finally, Figure 4 and Table 5 show that the results obtained from the deterministic Sigmoid and Brownian motion models strongly coincide: they suggest that CU was stable over the 1985-2003 period, remaining below 0.5, sharply increased between 2003 and 2007, and kept growing until 2009.

Overall, the results of this section support and strengthen the findings of various papers(Boily et al., 2007; Deering et al., 2008; Boily et al., 2008; Lowndes et al., 2010; Pickles et al., 2010; Boily et al., 2013; Pickles et al., 2013) which relied on very informative prior information from self reported estimates of condom use to assess the impact of the Avahan initiative on HIV infection averted. The implications of our analysis via the proposed methodology of the paper is the ability to draw inference on the space of CU trajectories without informing the model about the timing of the Avahan initiative or relying on other sources of data such as FSWs surveys of condom distribution data.