Evaluating the power of the causal impact method in observational studies of HCV treatment as prevention

Evaluation of HCV TasP using the CIM

Let t index the various waves of NESI, where t=1 corresponds to the 2008–2009 swap, t=2 to 2010, etc, and let i=0, …, n index the various units, where i=0 is the treated unit. In future NESI surveys (t > 5) it could be possible to evaluate the effect that HCV treatment scale-up had on virus prevalence by comparing p 0 t ( 1 ) , the prevalence at time t under the intervention in the treated site, to an estimate of the counterfactual p 0 t ( 0 ) , the prevalence that we would observe in the treated site if no intervention took place. That is, θ t = p 0 t ( 1 ) − p 0 t ( 0 ) , where θ _t is the causal effect of HCV TasP on prevalence at time t (t>5).

The CIM makes use of the data in the pre-intervention period (t≤5) and post-intervention data in the control sites to obtain estimates p ̂ 0 t ( 0 ) of the counterfactuals for t>5. It fits a Bayesian structural time-series model to the outcome in the pre-intervention period. Following standard modelling practice with prevalence data we choose to model y 0 t ( 0 ) = log p 0 t ( 0 ) 1 − p 0 t ( 0 ) instead of p 0 t ( 0 ) directly, and further assume that

where μ _t is the temporal local level component, y t = y 1 t , y 2 t ⊤ is the outcome (logit-prevalence) on control sites at time t,^[2] β = β 1 , β 2 ⊤ are the regression coefficients, ϵ t ∼ N ( 0 , σ ϵ 2 ) and δ t ∼ N ( 0 , σ δ 2 ) . In the model of Eqs. (1) and (2), the local level component μ _t induces temporal correlations in y 0 t ( 0 ) , and the regression component exploits correlations of y 0 t ( 0 ) with outcomes on the control sites.

The parameters of model (1) and (2) are estimated using Markov chain Monte Carlo (MCMC) techniques (Brodersen et al. 2015). Posterior simulations are simplified using the following conditionally conjugate prior distributions. We let σ ϵ − 2 , σ δ − 2 ∼ G a m m a ( ν 2 , s 2 ) . Brodersen et al. (2015) explain that ν can be thought of as prior sample size and s can be chosen such that s/ν is a guess for the variance. In practice, we can set s = ν ( 1 − R 2 ) σ ̂ y 2 , where σ ̂ y 2 is the sample variance of y 0 t ( 0 ) and R ² is the proportion of the variability in y 0 t ( 0 ) we expect to be explained by the regression component.

For β , a spike-and-slab prior (Chipman, George, and McCulloch 2001; George and McCulloch 1993, among others) is used. This prior assumes that for each β _i , there is a binary γ _i such that β _i ≠0 when γ _i =1 and β _i =0 otherwise. We let each γ _i ∼ Bernoulli(q _i ), where q _i is the prior probability that the coefficient of unit i is non-zero. The expected number of units with γ _i ≠ 0 using this prior is ∑ i = 1 n q i . Hence, we can set q _i =k/n for all i to encourage only k control units with γ _i =1 (β _i ≠0). Conditionally on γ = γ 1 , … , γ n ⊤ , let β _γ include the elements of β for which γ _i =1. We assume that β γ ∼ N 0 , σ ϵ 2 Σ β , where Σ _β is some prior covariance matrix (e.g. the identity matrix). The use and careful tuning of the spike-and-slab prior is important in epidemiological applications for two reasons. Firstly, by setting some γ _i =0 at each MCMC iteration, the method excludes controls whose data are not predictive of y 0 t ( 0 ) , and thus reduces the total number of parameters that need to be estimated. This is useful since the total number of pre-intervention time points is typically similar (or even smaller) than the total number of control units. Secondly, by calculating the posterior inclusion probability (the posterior mean of ρ _i ) for each control unit, it allows us to identify the ones that mostly contribute to the estimation of the counterfactual.

As additional data (t>5) become available, counterfactuals can be obtained by extrapolating the model of Eqs. (1) and (2). Let μ ̂ 5 , β ̂ , σ ̂ ϵ 2 , σ ̂ δ 2 be a sample from the posterior distribution of parameters μ 5 , β , σ ϵ 2 , σ δ 2 , respectively. A sample from the posterior predictive distribution of p 06 ( 0 ) , the counterfactual prevalence in the first post-intervention survey, will be p ̂ 06 ( 0 ) = exp y ̂ 06 ( 0 ) 1 + exp y ̂ 06 ( 0 ) , where

with δ ̂ 6 ∼ N ( 0 , σ ̂ δ 2 ) and ϵ ̂ 6 ∼ N ( 0 , σ ̂ ϵ 2 ) . Assume that we draw L such samples, p ̂ 06 , ℓ ( 0 ) (ℓ=1, …, L). Then, L samples from the posterior distribution of the causal effect at t=6 will follow as θ ̂ 6 , ℓ = p 06 ( 1 ) − p ̂ 06 , ℓ ( 0 ) , from which we obtain a point estimate (the mean of θ ̂ 6 , ℓ ) and a credible interval (the 2.5 and 97.5% percentiles of θ ̂ 6 , ℓ ).

Setting

Our objective is to assess the potential of the CIM for estimating the effect of HCV TasP using the existing UK HCV data (Section 2) combined with post-intervention data that will be collected. More specifically, we investigate the performance of the estimator of the causal intervention effect provided by the CIM method and identify the characteristics of the data that most affect the quality of the estimates of θ ̂ t . Our evaluation will also inform the potential of CIM for similar datasets. To achieve these goals, we performed a series of simulations.

First, we note that the performance of the estimator of θ _t depends solely on the performance of the estimator of p ̂ 0 t ( 0 ) , since θ t = p 0 t ( 1 ) − p 0 t ( 0 ) and p 0 t ( 1 ) is observed. More specifically, if p ̂ 0 t ( 0 ) is an unbiased estimate of p 0 t ( 0 ) , then θ ̂ t is an unbiased estimate of θ _t . Furthermore, then the 95% CI of θ _t will also include the true intervention effect if and only if the 95% CI of p 0 t ( 0 ) includes the untreated counterfactual. Therefore, it suffices to evaluate p ̂ 0 t ( 0 ) . We did this by considering the following performance measures at each post-intervention time point: (i) the mean (over simulated datasets) of the prediction error (MPE),^[3] where in each simulated dataset the prediction error is defined as the difference p 0 t ( 0 ) − p ̂ 0 t ( 0 ) ; (ii) the standard deviation (over simulated datasets) of the prediction error (sd-PE)); (iii) the mean (over simulated datasets) width of credible intervals of p 0 t ( 0 ) (CIW); (iv) the % false discovery rate (over simulated datasets) (FDR), where in each dataset a false detection occurred when p 0 t ( 1 ) = p 0 t ( 0 ) (i.e. θ _t =0) and the 95% CI of p 0 t ( 0 ) did not include p 0 t ( 1 ) ; and (v) the % detection rate (over simulated datasets) (power), where in each dataset a detection occurred when p 0 t ( 1 ) < p 0 t ( 0 ) (i.e. the intervention reduced prevalence) and the lower bound of the 95% CI of p 0 t ( 0 ) was higher than p 0 t ( 1 ) .

We simulated 10,000 datasets, each consisting of HCV logit-prevalence measurements y _it for n + 1 units and T times points. At each time point t, we drew y t = y 0 t ( 0 ) , y 1 t , … , y n t ⊤ ∼ M V N m , S that is the mean, variance and correlation of the outcomes remained constant over time. We drew the elements of m from a U n i f o r m m min , m max ; m _min/m _max is the minimum/maximum logit-prevalence found in the NESI dataset presented in Section 2. We set S 11 = σ y 2 . The remaining diagonal elements of S are drawn from a U n i f o r m s min 2 , s max 2 ; s min 2 / s max 2 is the minimum/maximum over all i of s i 2 , where s i 2 is the sample variance of the time-series of unit i in the NESI dataset. To obtain the off-diagonal elements of S , it suffices to pick the values of ρ _ij , the degree of correlation between the data on units i and j, where i , j ∈ 0 , … , n and i≠j. We set ρ _0j =ρ for all 1≤j≤k ₂ and ρ _0j =0 when j>k ₂. That is, the treated unit is only correlated to the first k ₂ controls units. Further, for all i, j such that 1≤i, j≤k ₂ and i≠j, we set ρ _ij =0.8ρ. Finally, for i>k ₂ (i.e. the k ₁=n − k ₂ control units that are not correlated to the treated unit), we have that ρ _ij =0 for all j≠i i.e. these units are not correlated to any other unit in the dataset. For each simulated dataset, we introduced intervention effects (i.e. obtained p 0 t ( 1 ) ) by reducing p 0 t ( 0 ) in the post-intervention period by a certain %. More specifically, we introduced 21 different effects from 0 to 50% with increments of 2.5%.

We attempted to generate data that mimic the HCV TasP application of Section 2. Therefore, we used the following simulation parameters in our baseline setup. The variance σ y 2 was set equal to the variance of the logit-prevalence measurements of the treated unit in the motivating dataset. Let T=t ₁ + t ₂, where t ₁ is the total number of pre-intervention data points per unit and t ₂ is the total number of post-intervention observations. We set t ₁ to be 6 and 12, t ₂=3 and n=8. In practical applications, we expect that only a small proportion of the control units to be correlated with the treated unit. Hence, in the baseline simulation we set k ₁=6 and k ₂=2. For the k ₂ ‘useful’ controls, we assumed that ρ=0.8.

In order to identify the features of the data that mostly affect the quality of causal estimates provided by the CIM, we performed several sensitivity analyses. In each sensitivity analysis we repeated the baseline simulation altering a single characteristic of the dataset and re-evaluated the five performance measures. The characteristics that we considered are (I) the variability of the outcome, of the treated unit σ y 2 ; (II) the total number of observations in the pre-intervention period, t ₁; (III) the total number of control units whose outcomes are not correlated with the outcome of the treated unit, k ₁; (IV) the total number of useful controls, k ₂; (V) the level of correlation between y _0t and the outcomes of the useful controls, ρ; and (VI) the hyperparameters of the spike-and-slab prior on the regression coefficients β . The values that we use for characteristics I–V are shown in Table 2.

Results

The MPE, sd-PE, CIW and FDR for the baseline simulations are shown in Table 4 of Appendix A. As can be seen in Table 4, the MPE at each post-intervention time point in the baseline setup was negligible (compared to the sd(MPE)), for both t ₁=6 and t ₁=12. This fact implies that, over the 10,000 simulated datasets, the estimates p ̂ 0 t ( 0 ) coincided on average with the corresponding ‘true’ values p 0 t ( 0 ) . It is confirmed in Figure 4 of Appendix A, where we plot the simulated values of p 0 t ( 0 ) against the estimated causal effect θ ̂ t . However, we see that there is positive correlation between p 0 t ( 0 ) and θ ̂ t i.e. the effect of the intervention is overestimated when the prevalence in the treated unit is high and underestimated when it is low. This correlation is expected due to the use of the logit transformation and drops with higher t ₁. We also see that the FDR is very close to the nominal 5% for t ₁=12, and slightly inflated for t ₁=6.

The power that we obtained at each t>t ₁ in the baseline simulations can be seen in Figure 1. As expected, the power increased with the % decrease in prevalence due to the intervention, and reached 100% when the intervention reduced prevalence by half. Lower drops in prevalence were associated with lower power. For example, a 10% decrease in HCV prevalence is only detected with probability 25 and 30% for t ₁=6 and t ₁=12, respectively. The power achieved was comparable at all three post-intervention time points. Nonetheless, it decreased with t. This decrease is due to fact that the variance of the random walk component is ( t − t 1 ) σ δ 2 (t>t ₁) which leads to wider credible intervals as t increases. Generally, in practical applications we expect that the uncertainty in the estimates of p 0 t ( 0 ) provided by the CIM will increase with t unless the time-series component has no contribution (this could be the case, for example, when all of the variability is attributed to the regression component).

Figure 1:

Baseline simulations results. The plots shows the power of detecting an intervention effect obtained by the CIM, as a function of the intervention effect magnitude. The left panel shows results for t ₁=6 and the right panel for t ₁=12. All results are based on 10,000 simulated datasets.

One of the advantages of the Bayesian approach is that several quantities of interest can be calculated directly from the posterior distribution of the model parameters. For example, rather than testing if θ _t is zero at each t>t ₁, one can use a summary to test for an overall effect. We examined the average causal effect in the post-intervention period defined as

A credible interval for ϑ that excluded zero was considered as evidence of an overall intervention effect in the entire post-intervention period. Figure 1 presents the power that we obtained when we tested an overall intervention effect, when this effect was constant (i.e. θ _t =θ for all t>t ₁). As expected, there were big gains in power when we summarised the information across all t ₂=3 post-intervention times. For example, for t ₁=6, a prevalence decrease of 20% was detected with probability 80% when we used ϑ to test for it but only with probability 60% when we examined each post-intervention time point individually. Hence, in practical applications, it is worth monitoring the outcome of the treated units on multiple time points after the intervention is introduced, as this can increase the chances to detect an intervention effect. Moreover, it might be worth considering the average effect only in the last s<t ₂ post-intervention time points since some interventions might not be effective immediately after introduction.

The results of our sensitivity analyses are summarised in Table 4 of Appendix A, Figure 2 and Figures 5 and 6 of Appendix A. Figures 2, 5, and 6 plot the power achieved by the CIM against the magnitude of the intervention effect at post-intervention times t ₁ + 1, t ₁ + 2 and t ₁ + 3, respectively. The MPE was negligible (compared to its standard error) across all sensitivity analyses and therefore is not further discussed. For the remaining performance measures, the results that we obtained for t=t ₁ + 1 were similar to the results obtained for t=t ₁ + 2 and t=t ₁ + 3. Hence, for the remainder of this section we focus attention to the first post-intervention time point.

Figure 2:

Results of the simulation study for the second post-intervention time point t=t ₁ + 1. The figure presents the power for detecting an intervention effect achieved by the CIM as a function of the intervention effect magnitude, in all baseline settings and sensitivity analyses. All results are based on 10,000 simulated datasets.

The value of t ₁ largely affected the performance of the CIM. As expected, increasing t ₁ caused all sd-PE and CIW to decrease (Figure 2(a)), since the parameters were estimated with higher accuracy. Further, the FDR was inflated for low values of t ₁. One possible explanation is that when t ₁ was low, it was more likely to observe strong correlations between the treated and a control unit by chance, thus assigning non-zero regression coefficients to control units whose outcomes were not truly correlated to the outcome of the treated unit. As expected, the variance of the outcome σ y 2 was crucial for the performance of the CIM. Larger outcome variance led to larger sd-PE and CIW (Table 4), and to a substantial drop in power (Figure 2(b)). Nonetheless, for fixed t ₁, the values of the FDR were similar across all values of σ y 2 considered.

The sd-MPE, CIW and power were not very sensitive the total number of ‘unrelated’ controls k ₁ (Figure 2(c)). With increasing k ₁, sd-PE and CIW increased, whereas the power dropped. The reason could be that there was a need to estimate more regression parameters as k ₁ increased. This effect was more prominent when t ₁=6. However, this drop in performance was negligible. We believe that this robustness to the addition of controls whose outcomes are not informative of the outcome of the treated unit is due to the spike-and-slab prior which successfully identifies these controls and, on average, set their coefficients to zero. This finding suggests that in real problems, since the expected drop in power is negligible, it is preferable to include all the available control units and allow the CIM to identify the ones that are important.

Increasing k ₂ slightly improved power but the gains were small, since each additional control could only explain a small proportion of the variability in y 0 t ( 0 ) that was not already explained by the existing ‘useful’ controls. Another factor that affected the quality of the causal estimates was the level of correlation between the outcome of the treated unit and the outcomes of ‘useful’ controls. This is expected since the method uses the regression component to exploit linear relationships in the data. Therefore, the stronger these relationships were, the higher was the proportion of the variability of y 0 t ( 0 ) explained by the regression component. As a result, sd-PE and CIW decreased, leading to an increased power. For small intervention effects, satisfactory power was only achieved for large values of ρ (Figure 2(e)).

Our final sensitivity analysis aimed to demonstrate the effects of the specification of the prior. Therefore, we repeated the baseline simulations, using a different spike-and-slab prior for the regression coefficients β . To this end, the informative spike-and-slab prior presented in Section 2 was replaced by the software default. Figure 2(f) presents the power under the two prior distributions. The power to detect an effect substantially dropped under the software default prior. This was the case because this prior was less informative compared to the prior that we initially used, and thus led to greater posterior uncertainty and therefore wider credible intervals for the untreated counterfactuals.

Effect on the performance of the CIM

The CIM assumes that the outcome of interest is observed without error in both the treated and the control units. This assumption is plausible in many real life problems, e.g. the one considered by Brodersen et al. (2015) where the outcome of interest is the total number of daily visits in various web-pages (the units) which can be precisely enumerated. Other examples of outcomes that can be measured without error (or estimated very precisely) include the daily sales of a product in a geographical region, the total number of deaths due to a disease in a hospital and the annual GDP of a country. However, in many epidemiological studies, it might not possible to observe the outcome without some error. For example, in the motivating application of Section 2, the true HCV prevalence in each unit is unknown and it is estimated through surveillance data as p ̃ i t = k i t N i t , where k _it and N _it represent the total number of infected individuals and the total sample size from the surveillance study in unit i at time t, respectively. Note that we refer to p ̃ i t as imprecise prevalence in order to distinguish it from the estimated prevalence p ̂ i t ( 0 ) obtained from the CIM.

To assess the impact that the use of imprecise outcomes (instead of the true, unknown outcomes) had on the performance of the CIM, we re-analysed the same 10,000 simulated datasets that we analysed for the baseline simulation of Section 3, when t ₁=12 and t ₁=24. Instead of implementing the CIM to y i t = log p i t 1 − p i t , we implemented it to y ̃ i t = log p ̃ i t 1 − p ̃ i t , where p ̃ i t = k i t N i t and k _it were simulated from a B i n N i t , p i t distribution. We evaluated the performance considering the same performance measures as in Section 3. We only present the power at the first post-intervention time point because we found that the results were very similar in the remaining post-intervention time points. We artificially introduced the intervention effects that were non-zero, by drawing the k _0t from a B i n N i t , p 0 t * , where p 0 t * were obtained by reducing the original prevalence p _0t . The sample size across units and time points was constant, i.e. N _it =n for all units i and times t. We simulated n=50, n=50 and n=150.

Table 3 (sd-PE, CIW and FDR) and Figure 3 (power) summarise the results for t ₁=12; the results for t ₁=24 are similar and therefore not shown. For comparison, we also show the results that we obtained in the baseline simulation when we assumed outcomes were measured without error. As expected, the use of the imprecise outcomes p ̃ i t instead of the perfectly measured outcomes p _it degrades the performance of the CIM substantially. Table 3 shows that both the sd-PE and CIW increase when the CIM is implemented using p ̃ i t . For example, the CIW obtained when the sample size n=50, was approximately twice the width that we obtained using the original data. As a result, there was also reduced power to detect an intervention effect (Figure 3). For instance, the power to detect a 25% decrease using the approximated outcomes and n=50 was roughly 37%, as opposed to 75% for the exact prevalence outcomes.

Figure 3:

Effect of measurement error on the performance of the CIM. The figure presents the power achieved for t=t ₁ + 1 when the CIM and CIM-EIV are applied to the imprecise outcomes p ̃ 0 t . For reference, we show results for the CIM when implemented to the true outcomes (CIM, n=∞). The left and right panels correspond to t ₁=12 and t ₁=24, respectively. Results are based on 10,000 simulated datasets.

The increase in uncertainty (and therefore loss of power) occurred because y ̃ i t are noisy observations of y _it and therefore the correlations between y ̃ 0 t and y ̃ i t (i>0) were weaker than the correlations between y _0t and y _it in the original simulation study. As a result, the estimates of regression coefficients of the k ₁ predictive control units were biased downwards and the estimates of σ ε 2 were biased upwards. In the classic linear regression setting this phenomenon is known as regression dilution, see e.g. Frost and Thompson (2000).

In addition to the increased uncertainty in the estimates of the causal effect, the use of imprecise measurements also led to an increased FDR (Table 3). More specifically, for n=50, the FDR at t=t ₁ + 1 was roughly 16%, more than triple the desired nominal level of 5%. A potential explanation is that some control units appeared to be highly correlated with the treated unit in the pre-intervention period by chance, because of the error in p ̃ i t . As a result, the coefficients of these units were over-estimated, leading to inaccurate prediction of the untreated counterfactual in the post-intervention period.

Both of these problems, i.e. increased uncertainty in the causal estimate and increased false positive rate, became more profound when the sample size n was reduced. The reason is that as n decreased, the y ̃ i t become more variable (i.e. the measurement error increased).

An errors-in-variables causal impact method

The simulation study of Section 4.1 shows that measurement error has an adverse impact on the performance of the CIM, reducing power and increasing the false positive rate. The former is expected and inevitable when it is not possible to measure the outcome precisely. However, an increased false positive rate is an undesirable property which can reduce the reliability of a significant finding obtained using the CIM, especially when the estimated intervention effect is small. In this section, we extend the CIM in an attempt to deal with this problem.

We propose a two-level Bayesian hierarchical. At the first level, we have the data. Let k 0 t ( 0 ) be the total number of infected individuals in the treated units when there is no intervention. We have that k 0 t ( 0 ) = k 0 t when t≤t ₁ and missing for t>t ₁. Further, let k i t ( 1 ) = k i t (t>t ₁) be the total number of infected individuals in the treated unit when the intervention is in effect. We assume that

Equation (5) relates the logit-prevalence to the observed data thus acknowledging that there is uncertainty regarding its true value. The smaller N _it is, the larger the uncertainty regarding the true value of y _it . Depending on the application, one might need to adopt observation Eq. (5) in order to account for more complex relationships between the observed data and the logit-prevalence (e.g. when there are data from multiple sub-populations of individuals).

At the second level, we have the unknown prevalence parameters. Similar to the CIM, we assume that the untreated logit-prevalence y 0 t ( 0 ) in the treated unit can be written as

where ε t ∼ N o r m a l ( 0 , σ ε 2 ) for all t. For the treated prevalence p 0 t ( 1 ) in the treated unit we a priori assume that p 0 t ( 1 ) ∼ B e t a ( 1,1 ) for all t>t ₁. For σ ϵ 2 and β we use the same prior specifications as in Section 2.1. The intercept α _t arises from an AR(1) process i.e.

where ϕ ∈ (−1, 1) is the persistent parameter and η t ∼ N ( 0 , σ η 2 ) . For the AR hyperparameters μ, σ η 2 and ϕ we use similar priors as Kastner and Frühwirth-Schnatter (2014). More specifically, we let μ ∼ n(0, 10³), σ η 2 ∼ G a m m a ( 0.5 , 0.5 ( 1 − R 2 ) σ ̂ y 0 2 ) and ϕ + 1 2 ∼ B e t a ( 1,1 ) , where R ² and σ ̂ y 2 are defined as in Section 2.1.

Samples θ _t,ℓ (ℓ=1, …, L and t>t ₁) from the posterior distribution of the causal effects are obtained as p 0 t , ℓ ( 1 ) − p 0 t , ℓ ( 0 ) . The p 0 t , ℓ ( 1 ) are drawn from their Beta(1 + k _0t , 1 + n_0t − k _0t ) posterior distributions. The p 0 t , ℓ ( 0 ) are drawn from their posterior predictive distributions via MCMC. The proposed algorithm is a block Gibbs sampler that is on each iteration, one parameter (or block of parameters) is drawn from its full conditional distribution given the remaining parameters and data. The indicator variables γ _i are drawn one at a time, see e.g. Sutton (2020). The AR hyperparameters μ, σ η 2 and ϕ are jointly updated using a Metropolis-Hastings step (Kastner and Frühwirth-Schnatter 2014). The unknown logit-prevalence y i = y i 1 , … , y i t 1 ⊤ are drawn one at a time from their normal full conditionals; for this to be possible we make use of the Pólya–Gamma representation of the Binomial likelihood as proposed by Polson, Scott, and Windle (2013). The remaining model parameters α = α 1 , … , α t 1 ⊤ , β and σ ε 2 have conjugate prior distributions and are therefore easy to update. The code that we used has been made publicly available.^[4]

In the model of Eqs. (5)–(7), the covariates y _t are random variables. Hence, the model is similar in spirit to errors-in-variables (EIV) models often used to deal with the problem of regression dilution in practice, see for example Dellaportas and Stephens (1995). We therefore refer to this model as the CIM-EIV approach. However, it is more general than an EIV model as it allows for the response variable y _0t to be measured with error as well.

Application to the simulated data

We applied the proposed CIM-EIV method to the data that we simulated for the experiment of Section 4.1. We used the same R ² and σ ̂ y 2 (for the priors on the variance parameters) as we did for the CIM. The prior distributions for the spike-and-slab parameters were the same as for the CIM with the exception that Σ _β was set to 10³ I .

The sd-PE, CIW and FDR are presented in Table 3. We see that for fixed n, the sd-PE obtained by the CIM-EIV was lower compared to the one obtained by the CIM. However, the proposed method successfully adjusted for the uncertainty regarding the true values of the prevalence thus leading to wider credible intervals. As a result, we see that the proposed EIV approach reduced the FDR compared to the CIM, and that the benefits were more apparent when n was small. When we increased n, the magnitude of the difference y ̃ i t − y i t decreased and therefore the CIM-EIV did not improve much compared to the CIM (whose performance in terms of the FDR was already satisfactory). Therefore, we recommend that the CIM-EIV method is used especially in cases where the problem of dilution is expected to be high.

Note that the power of the CIM-EIV method was lower compared to the power of the CIM (see Figure 3). This is expected, since the CIM-EIV relaxes the assumption of the CIM that the outcomes are known precisely. In order to increase the power, one can combine post-intervention time-points as explained in Section 3.

Main findings

Using an HCV treatment as prevention intervention as a case study, our paper presents a series of simulations studies to investigate the potential of the CIM for use in observational epidemiological/public health studies aiming to estimate the causal effect of an intervention on an outcome of interest using aggregate time-series observational data. Overall, our experiments show that if the untreated outcome of the treated unit is linearly related to the (untreated) outcomes of some of the controls units and the effect of the intervention is effective, then the method will provide satisfactory power. We have found that the main characteristics of the data that affect the ability of the CIM to detect a non-zero intervention effect are the length of the time-series in the pre-intervention period, the variability of the outcome and the strength of the linear relationships between the pre-intervention data on the treated unit and the control units.

This work has demonstrated some of the potential merits of adopting a Bayesian approach for this problem. In particular, we have shown that it is possible to improve power by summarising information from all post-intervention time points rather than considering each one separately. Moreover, our simulation experiments suggest that if the prior distributions for the CIM model parameters are not chosen carefully then the method may provide misleading results. Finally, we have studied the implications of prevalence being measured with error on the performance of the CIM. Specifically, our simulations show that when the prevalence is estimated based on a small sample of individuals, the power of the method drops substantially and the false positives rates are inflated. In such cases, it might be preferable to use the proposed CIM-EIV approach.

Our work has important implications for HCV elimination initiatives and HCV TasP researchers. Theoretical modelling studies have shown the substantial potential benefits of scaled-up HCV treatment for PWID on reducing HCV chronic prevalence and incidence (Cousien et al. 2014; Durier, Nguyen, and White 2012; De Vos and Kretzschmar 2014; Hellard et al. 2014; Martin et al. 2011; Martin, Miners, and Vickerman 2012a; Martin et al. 2013a, 2013b; Martin et al. 2016a, 2016b, 2016c; Rolls et al. 2013; Vickerman, Martin, and Hickman 2011; Zeiler et al. 2010). However, empirical studies are needed to confirm that HCV treatment as prevention expansion can yield population declines in prevalence and incidence. As randomized controlled trials testing HCV TasP may be logistically difficult, prohibitively expensive, or ethically questionable, observational studies may provide alternative evidence for a TasP effect. Our findings, that the CIM method is a robust method for detecting a TasP intervention effect using surveillance data from the UK, provide an important methodological tool for use in empirical evaluations of HCV TasP using observational data. Indeed, ongoing observational studies of HCV treatment expansion among PWID such as occurring in Dundee and across the UK as part of the EPiTOPE (Hickman et al. 2019) study will, when combined with CIM methods, shed important new information on the effectiveness of HCV TasP in the real-world.

Since our simulated data have been generated based on an existing UK HCV dataset regarding the effectiveness of TasP against the HCV, we expect that our conclusions will be relevant to other public health applications in clinical practice.

Limitations and future research

This work has limitations. First, in our simulation experiments we have assumed that the mean untreated logit-prevalence remains constant over time in both control and treated units. Further, we have assumed that the correlation between the logit-prevalence of the treated unit and the logit-prevalence of control units lalso remains constant. Hence, in the future, it is worth studying the performance of the CIM (in terms of both bias and power) under data generating mechanisms where these assumptions do not hold. This could be done, for example, by introducing a declining trend in a subset of the units. Second, in future research, it is worth comparing the performance of the CIM with other existing methodologies, such as difference-in-differences and generalised linear mixed models, since the results from existing comparative studies (Gobillon and Magnac 2016; Kinn 2018; O’Neill et al. 2016; O’Neill et al. 2020) may not generalise to the type of data that we consider.

There are many ways in which the proposed CIM-EIV approach can be improved. One idea is to account for the fact the unknown prevalence in control units is likely to show serial correlations. For example, one could assume that the logit-prevalence in control units is an AR(1) process. Another option is to account for correlations between controls units. Both of these extensions are likely to improve the precision of the causal estimates provided by the method.

Finally, we note that we use UK surveillance data to construct our case study, which incorporates regular, routine surveillance among PWID. In many settings, surveillance among PWID occurs more sporadically, or among fewer sites, or does not occur at all. In these settings, CIM methods may not generate sufficient power to detect an intervention effect, or the observational period may need to be lengthened. Further studies in different settings with alternative surveillance systems are warranted.