Local prevalence of transmissible SARS-CoV-2 infection : an integrative causal model for debiasing fine-scale targeted testing data

Targeted surveillance testing schemes for SARS-CoV-2 focus on certain subsets of the population, such as individuals experiencing one or more of a prescribed list of symptoms. These schemes have routinely been used to monitor the spread of SARS-CoV-2 in countries across the world. The number of positive tests in a given region can provide local insights into important epidemiological parameters, such as prevalence and effective reproduction number. Moreover, targeted testing data has been used inform the deployment of localised non-pharmaceutical interventions. However, surveillance schemes typically suffer from ascertainment bias; the individuals who are tested are not necessarily representative of the wider population of interest. Here, we show that data from randomised testing schemes, such as the REACT study in the UK, can be used to debias fine-scale targeted testing data in order to provide accurate localised estimates of the number of infectious individuals. We develop a novel, integrative causal framework that explicitly models the process underlying the selection of individuals for targeted testing. The output from our model can readily be incorporated into longitudinal analyses to provide local estimates of the reproduction number. We apply our model to characterise the size of the infectious population in England between June 2020 and January 2021. Our local estimates of the effective reproduction number are predictive of future changes in positive case numbers. We also capture local increases in both prevalence and effective reproductive number in the South East from November 2020 to December 2020, reflecting the spread of the Kent variant. Preparations for future epidemics should ensure the rapid deployment of both types of schemes to accurately monitor the spread of emerging and ongoing infectious diseases.


Introduction
The spread of the new severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and the ensuing outbreaks of coronavirus disease 2019 (COVID- 19) have placed a significant burden on To understand the ascertainment bias problem and a statistical approach to correction, it is helpful to consider a simplified causal model for Pillar 1+2 data. This is represented by a directed acyclic graph (DAG), shown in Figure 1(a), that charts the dependencies of an individual from infection status to test result. The circles indicate the binary (Yes/No) states of an individual. 100 The DAG characterises the joint distribution of the major factors leading to the observed data.
Throughout the paper we use the term targeted testing data to refer to data gathered under some ascertainment process distinct from (stratified) random sampling, with an exemplar being selection for testing of the subpopulation with COVID-19 symptoms, which comprises a sizeable proportion of Pillar 1+2 tests. The Pillar 1+2 DAG can be compared to that of a randomised surveillance 105 study (shown in Figure 1(b)). The randomised nature of the test allocations in REACT renders Tested conditionally independent of Symptoms given Infected, yielding unbiased estimates of Infection rates. The DAG explicitly characterises statistically why we cannot use Pillar 1+2 data directly. The DAG also points to a potential solution if the statistical dependencies as indicated by the arrows in Figure 1(a) can be modelled, then Pillar 1+2 data can be used. In this paper, we 110 describe an approach allowing characterisation and adjustment for the ascertainment bias inherent in Pillar 1+2 data.
In addition to prevalence, there are a number of epidemiological parameters that may be useful for informing localised NPIs. For example, one particular variable of interest is the (time-varying) effective reproductive number R t , defined roughly as the average number of infections caused 115 by an infectious individual; when R t > 1, the epidemic will continue to spread. Estimation of these parameters relies on careful mathematical modelling supported by relevant data such as surveillance testing results or number of hospitalisations [15]. Incorporating multiple sources of data can produce more reliable parameter estimates [16], though doing so in a computationally efficient manner may be nontrivial [17]. For example, a stochastic epidemic model of the 2009 120 influenza outbreak in Finland that used data on hospitalisations, lab tests and vaccination data provided insights into the time-varying nature of R t but took over a month of compute time to run [18]. Given the time-sensitive nature of the current epidemic, one important modelling consideration is the timely inference of parameters; the work we present here has been developed with both accuracy and computational efficiency in mind. 125 The current pandemic has spurred the development of a number of models that also aim to incorporate multiple sources of data in order to estimate important epidemiological parameters, in particular R t [19]. Abbott et al. [20] generate daily estimates of R t at a national and PHE region level by incorporating case counts and death notifications, building on a model to estimate R t from incidence time series [21]. Birrell et al. [22] estimate daily, PHE regional R t using ONS CIS data, Targeted test-and-trace data (Pillar 1+2); and (b) Randomised surveillance data (e.g. REACT).
In (b), randomisation breaks the causal link between COVID-19 symptoms and swab testing. The nodes represent binary (yes/no) states for an individual in the relevant population.
5 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The research referenced so far infers epidemiological parameters at spatially coarse scales, such as PHE region. To extend this body of work, we focus here on accurate estimation of prevalence (along with other epidemiological parameters) at a more local level, such as the LTLA.
There are two very useful websites providing LTLA-level up-to-date estimates and predictions of some epidemiological parameters (but not prevalence). 1 A team at Oxford has produced a 140 local Covid map 2 as part of the Royal Society Data Evaluation and Learning for Viral Epidemics (DELVE) initiative. 3 Their methods take as input Pillar 1+2 daily counts and commuter flow data, and output local estimates and predictions of R t and positive case numbers. An Imperial College team has produced a COVID-19 UK map and table. 4 Their methodology takes as input: daily cases data, weekly deaths data, as well as daily infections from the ONS CIS and REACT 145 data sets [26]. They output estimates and predictions of R t , positive case numbers, and change in new infections. Their results are based on the epidemia software [27], which is an extension of the Bayesian semi-mechanistic model introduced in [28], though detailed methods are not yet available. We have downloaded their R t estimates and find a high level of consistency in local R t estimates between our model and theirs. A group from Lancaster University has estimated 150 daily case prevalence (proportion of infected population), incidence, and R t at the Local Authority District level in England (315 areas in total) by building an epidemic model incorporating measures of human mobility and Pillar 1 and 2 tests across England [29]. An important aspect of this approach is that it assumes each infection is eventually reflected in the Pillar 1 and 2 case reports and so does not account for possible ascertainment bias in targeted testing. Furthermore, the 155 model requires substantial computational resources to obtain timely estimates.
Within this urgent and fast developing area of research, it is clearly important to define the aspects in which our method contributes novelty. Firstly, we have developed methods to infer local prevalence, I t , accurately from targeted testing data. Here we work with weekly period prevalence, and explicitly target the number of infectious individuals via a correction to the estimated PCR-160 positive numbers. This is all novel and important in its own right -being able to estimate local prevalence accurately from targeted testing data adds an important facet to existing COVID-19 monitoring capabilities. Second, our method outputs bias-adjusted cross-sectional prevalence likelihoods p(n t of N t | I t ), where n t and N t are positive and total targeted test counts. This allows prevalence information from targeted data to be coherently embedded in a modular way into 165 complex spatiotemporal epidemiological models, including those synthesising multiple data types.
We exemplify this by implementing an Susceptible-Infectious-Recovered (SIR) model around our ascertainment model likelihood. Third, our local ascertainment model is based on targeted testing data alone with, uniquely to our knowledge, both the number of positive and total tests being modelled (n t and N t ). This has two important benefits: spatiotemporal variation in testing uptake 170 1 Technical details of the methodology driving these websites is not yet available at the time of writing, but in both cases the peer-review process is underway.

Results
Correcting for ascertainment bias in targeted testing data With reference to the causal DAG in figure 1, we define the essential bias parameter, δ, as 180 δ := log Odds(Tested | Infected) Odds(Tested | Not Infected) (1) i.e. the log odds ratio of being tested in the infected versus non-infected populations. Larger values of δ generally correspond to higher levels of ascertainment bias, i.e. a higher chance of an infected individual being selected for testing, relative to a non-infected.
Our approach combines randomised surveillance data (REACT) and targeted surveillance data (Pillars 1 and 2) to infer δ at the coarse geographical level (PHE region). We then integrate this 185 information by specifying a temporally smooth empirical Bayes (EB) prior on δ 1:T , applied to each constituent local region (LTLA) in the local prevalence analyses. Figure 3 shows the resulting EB priors on δ; there is potentially more variation in δ across regions early on in the sampling period (pre-September 2020), though the prior credible intervals (CIs) are quite broad and often overlapping. The data provide more information on δ from October 2020 onwards, and there is a 190 consistent upward trend for all nine PHE regions.

Cross-sectional local prevalence from targeted testing data
De-biased likelihood for modular sharing of prevalence information Equipped with a coarse-scale (PHE-region level) EB prior on bias δ we evaluate a fine-scale (LTLA-level) δ-marginalized likelihood of the form p(n t of N t | I t ,ν t ), as described at (17)

Longitudinal local prevalence and transmission
The cross-sectional de-biased likelihood can be introduced modularly into a wide variety of downstream epidemiological models. We illustrate this by using the likelihood as an input to a simple SIR epidemic model (see Methods-Full Bayesian inference under a stochastic SIR epidemic model and Figure 12). Figure 5(a) plots estimated prevalence against R t number at the most recent time . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) Relating local prevalence and transmission to spread of the UK variant One striking feature of the maps in Figure 6 is the increasing prevalence in the London area throughout November to December 2020. This is consistent with the known arrival of the UK variant of concern (VoC) 202012/01 (lineage B.1.1.7), that emerged in the South East of England 230 in November 2020 and which has been estimated to have a 43-90% higher reproduction number than preexisting variants [30].
We investigate this hypothesis similarly to [30], by characterising the relationship between estimated local R t and the frequency of VoC 202012/01, as approximated by the frequency of S gene target failure (SGTF) in the Taqpath sequencing assay used over this time period [31].
235 Figure 8 illustrates the spatial distributions of VoC 202012/01 against estimated prevalence and estimated R t from mid-November 2020 to mid-December 2020. The increase in frequency of the VoC was initially isolated to the South-East but then spread outwards, accompanied by a corresponding increase in both local estimated prevalence and R t . We observe a strong positive association between the local VoC frequency and estimated local R t , consistent with the increased 240 transmissibility identified in [30].

10
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

11
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Validation 1 (accuracy)
We qualitatively assess the performance of de-biased fine-scale (LTLA-level) prevalence estimates by measuring how well they predict LTLA-level REACT data. The validation is best described in terms of coarse-scale REACT training data and contemporaneous fine-scale REACT test data.

Validation 2 (prediction)
The effective reproduction number, R t , measures whether the number of infectious individuals 255 is increasing, R t > 1, or decreasing, R t < 1, in the population at time point t. Figure 9 compares LTLA R t estimates with the future change in local case numbers. For validation purposes, here we are doing one-step-ahead at a time prediction and comparing predictions with out-oftraining-sample observed statistics (fold change in raw case numbers from baseline). The results 12 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

13
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Validation 3 (external)
We extracted estimates of effective reproduction number R t based on our de-biasing model likelihood implemented within a standard SIR model, illustrated in Figure 12. We compare the results 270 to the local R t estimates outputted by at the Imperial College COVID-19 website. 5 A crossmethod comparison of longitudinal traces of R t for a subset of LTLAs is shown in Figure 10.
Encouragingly for both approaches, the estimates generally display good concordance, with credi- Effective R value Figure 10: Comparison of R t estimates between de-biasing model and Imperial model [32].
For each of the nine PHE regions, we present the constituent LTLA whose name is ranked top alphabetically.
ble intervals overlapping appropriately, despite being based on different data and models. 6 .

275
Observational models for surveillance data  6 The imperial model uses daily cases data, weekly deaths data, as well as daily infections from the ONS CIS and REACT data sets 15 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint and this allows direct, accurate statistical inference onĨ, the proportion of the population that would return a positive PCR test.
Focusing prevalence on the infectious subpopulation PCR tests are sensitive, and can detect the presence of SARS-CoV-2 both days before and weeks after an individual is infectious. It is usually desirable for prevalence to represent the proportion of a population that is infectious. We can obtain a likelihood for the number of infectious individuals I as follows: where I andĨ are the number of infectious and PCR-positive individuals respectively.
where P(Tested | Infected) and P(Tested | Not Infected) are the probabilities of an infected (respectively non-infected) individual being tested on date t.
We introduce the following parameters: δ := log Odds(Tested | Infected) Odds(Tested | Not Infected) (5) ν := log Odds(Tested | Not Infected) , 16 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint leading to the targeted swab testing likelihood being represented as P(n of N | I, δ, ν) = Binomial n | I, logit −1 (δ + ν) × Binomial(N − n | M − I, logit −1 ν) . (7) The unknown parameter requiring special care to infer is δ, i.e. the log odds ratio of being tested in the infected versus the non-infected subpopulations. The other parameter, ν, is directly estimable 310 from the targeted data:ν := logit[(N −n)/M ] is a precise estimator with little bias when prevalence is low.
Test sensitivity and specificity.
The likelihood at (7) assumes a perfect antigen test. If the test procedure has false-positive rate α, and false-negative rate β, the targeted likelihood is instead where z denotes the unknown number of truly infected individuals that were tested. The first term in the sum at (8) is obtained by substituting z in (7), while the second term is with n β denoting the number of false-negative test results. An analogous adjustment can be made to the randomised surveillance likelihood at (2).
Cross-sectional inference on local prevalence 320 We leverage spatially coarse-scale randomised surveillance data to specify an EB prior on bias parameters p(δ) at coarse-scale (PHE region), and thereby infer prevalence accurately from targeted data at fine scale (LTLA j within PHE region J j ). We explicitly use the superscripts LTLA (j) in PHE region (J j ) in step 4 below where notation from both coarse and fine scale appear together.
All quantities in steps 1-3 are implicitly superscripted (J j ) but these are suppressed for notational 325 clarity. For computational efficiency we handle prevalence in a reduced-dimension space of bins as described in Supplementary Information (SI) section Interval-based prevalence inference -set-up and assumptions. The method in detail is as follows: 1. Infer prevalence from unbiased testing data. At a coarse geographic level (PHE region J j ), estimate prevalence from randomised surveillance data u t of U t . Represent the posterior 17 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint 2. Learn δ t from accurate prevalence. At a coarse geographic level, for each t ∈ T , we estimate bias parameter δ t by coupling biased data n t of N t with accurate prevalence infor- where a moment-matched Gaussian approximation is performed at (12). The posterior density in the sum at (11), p(δ t | n t of N t , I t ,ν t ) is conjugate under a Beta(a,b) prior on logit −1 (ν t + δ t ) ≡ P(Tested | Infected), and so can be evaluated as 3. Specify smooth EB prior on δ 1:T . A smooth prior on δ 1:T is specified as follows where N(δ | 0, Σ δ ) imparts a user-specified degree of longitudinal smoothness, thereby sharing information on δ across time points. Ignorance of δ t , in the absence of random surveillance data, is encapsulated in a Gaussian with large variance σ 2 flat . A standard choice for N(δ | 0, Σ δ ) corresponds to a stationary autoregressive, AR(1), process of the form with a diffuse Gaussian prior c ∼ N(0, σ 2 flat ) and with smoothing tuned by 0 < ψ < 1 and 345 white noise variance σ 2 ε . The normalised form of the prior at (14) is with (μ, diagonal matrix D T ×T ) having elements (μ t ,σ 2 t ) for t ∈ T and (0, σ 2 flat ) for t ∈ T .

4.
Infer cross-sectional local prevalence from biased testing data. At a fine-scale geographic level (LTLA j in PHE region J j ), having observed n t ], the likelihood in the integral at (18) is available at (7), and the prior p(δ

Debiasing lateral flow device (LFD) tests with PCR surveillance (or vice versa)
The methods can be adapted straightforwardly to the situation in which the randomised surveil-355 lance study uses a different assay to the targeted testing. For a concrete example we could use 18 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint REACT PCR prevalence posteriorp t (Ĩ t ) from (10) to debias Pillar 1+2 LFD test data n t of N t . Equation (11) can be adjusted to estimate ascertainment bias δ pertaining to LFD data as follows: whereĪ t andĨ t are the unobserved LFD-and PCR-positive prevalence respectively, and the con-

Full Bayesian inference under a stochastic SIR epidemic model 365
The cross-sectional analysis described in Cross-sectional inference on local prevalence generates the (17), at each time point for which targeted data are available. These likelihoods can be used as input for longitudinal models to obtain better prevalence estimates and to infer epidemiological parameters such as R t . Priors on R, I, R + We place priors on I, R + measured as a proportion of the population; this proportion then gets mapped to prevalence intervals on subpopulation counts as described in Interval-based prevalence

19
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021 Reproduction value, R t  Figure 11. Inference is for a region, e.g. an LTLA, based only on targeted test data collected in this region, n t of N t . A prior on δ t parameterized (μ t ,σ 2 t ) brings information on the Pillar 2 ascertainment bias learned from randomized surveillance testing data available for the PHE region in which the LTLA lies. The T × T covariance matrix Σ δ imparts temporal smoothness on δ 1:T . Effective reproduction numbers are denoted R 1:T , number of infectious individuals by I 1:T , and the number of immune individuals by R + 1:T .

20
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint inference -set-up and assumptions. Specifically, we use truncated, discretized Gaussian distributions on the proportion of the population immune and infectious. For example, on number of infectious individuals I t at each timepoint t, we specify the prior (suitably normalized over its support) with an example weakly informative hyperparameter setting being µ I = 0.5%, σ I = 1%, p min = 0%, p max = 4%. To ensure meaningful inference on R + 1:T , we place an informative prior that reflects the state of knowledge of the immune population size; we do this using an informative truncated Gaussian prior on R + 1 , and non-informative priors on R + 2:T . We place a noninformative 380 uniform prior on each R t , e.g. a Uniform(0.5, 2.5).

MCMC sampling implementation
We perform inference under the model represented in the DAG at Figure sampling I new from p(I | R, n, N ), and then R + new from p(R + | I new ).
Sampling from p(I | R, n, N ) The sampling distribution on prevalence can be expressed: which is an HMM with emission probabilities taken from the δ-marginalised likelihood at (18), and transition probabilities taken from (37).

Sampling from p(R | I)
The prior joint distribution of R 1:T is modelled using a random walk: where σ 2 R is a user-specified smoothness parameter. The update involves sampling from We discretize the space of R t into an evenly spaced grid and sample from the HMM defined at (24) [36]. The transition probabilities are given by (23) (suitably normalised over the discrete R t space) and the emission probabilities given by (37).

Data
With the exception of the variant of concern 202012/01 analysis, all data underlying the results 400 presented here are publicly available. Randomised surveillance data comes from the REACT study

Analysis scripts
The R scripts used to generate the results in this manuscript are available at https://github. com/alan-turing-institute/jbc-turing-rss-testdebiasing.

Discussion
We have introduced and applied an integrative causal model allowing accurate inference from 410 community testing data. The flexible probabilistic framework allows simultaneous and coherent incorporation of a number of important features, including: • Adjustment for ascertainment bias caused by preferential testing based on symptom status, or on other confounders.
• Allowing for heterogeneous testing capacity by modelling the total number of tests conducted 415 locally.

22
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021 • The model outputs week-specific debiased prevalence with uncertainty (via a marginal likelihood), that can be incorporated modularly into more complex models.
• An SIR epidemic model implementation allowing estimation of R t and adjustment for vaccination in the immune population Because of the extensive Pillar 1+2 testing effort, there is a large amount of information con-425 tained in these targeted data at LTLA level, even for a single weekly timepoint, as shown by the narrow width of CIs in Figure 4. However, due to the strength of this information, the targeted data also have the potential to introduce bias into more complex models if they are incorporated as observed nodes without an appropriate ascertainment correction. Equipped with a well specified prior on δ, however, precise and accurate estimation at even finer scales may be feasible. Given 430 the high volume of data, despite focusing on high resolution geographical units such as LTLAs, estimates do not seem to be affected by the classical issue of small area estimation, that is to say sample sizes too small to carry on inference without borrowing strength from neighbouring units.
If interest is focused on ultra fine-scale geography and/or when prevalence is much lower, then spatial borrowing would be beneficial, with such additional smoothing subject to a variance-bias 435 trade-off.
In addition to estimating at finer spatial scales, it may be desirable to estimate local prevalence from targeted data stratified by factors such as age and ethnicity. This would involve modelling the relevant confounders in Figure 1, rather than marginalising them out. One approach to doing this within the existing framework would be to perform a stratified analysis (e.g. performing the 440 whole analysis using weekly PHE region and Pillar 1+2 Data stratified by age bands). A more sophisticated approach could model δ t semi-parametrically with some spatiotemporal smoothness assumptions on the effect of age and other confounders. Of course, any approach would require appropriate metadata on any confounders that affect randomised surveillance and/or targeted sampling.

445
For cross-sectional prevalence estimation, a key dependency is the availability of a regular, upto-date stream of randomised surveillance data at some level of spatial resolution. Here we deployed REACT data at the coarse PHE-region scale. The UK has led the way internationally in having regular national surveillance randomised surveys like REACT and ONS. This modelling work shows the importance of having both targeted testing and also a rolling randomised surveillance 450 survey to be able to better track the epidemic. The methods are transferable beyond the UK wherever randomized testing data are being gathered. can be applied in other This could be built in an integrated way from the start as preparedness for pandemics, in particular for diseases where asymptomatic transmission plays an important role.

23
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

25
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
Coronavirus  28 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint A full list of model parameters, along with either their prior distribution or the value at which they were fixed, can be found in Table S1.  PCR false-negative rate, β, Fixed, β = 0.05, taken from [39] Expected time to recovery, 1/γ T recovery ∼ Exponential(γ), with γ = 1 week Effective reproduction number, R t Random walk:

Discussion of methodological assumptions and caveats
Interval-based prevalence inference -set-up and assumptions The full prevalence state space comprises all potential numbers of infectious individuals in the

29
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint and make three assumptions to allow computationally efficient inference on the B-dimensional space of bins, denoting these assumptions Interval-1:3 as follows: Interval-1 The testing data likelihood, conditional on prevalence bin, is evaluated at the bin mid-615 point: Interval-2 Prevalence I is uniformly distributed within each bin: otherwise.
Interval-3 The distribution of new infections, conditional on prevalence bin, is evaluated at the bin midpoint (with the same assumption applying to new recoveries): Ascertainment bias model -assumptions and caveats 620 Debias-1 Spatial homogeneity of δ across LTLAs within a PHE region. The fact that we see relatively low variation in δ at each time point across PHE regions in Figure 3, particularly after October 2020, is consistent with a finer-scale spatial homogeneity assumption being reasonable.
Debias-2 We handle prevalence in a reduced-dimension space of bins as described in SI section 625

Interval-based prevalence inference -set-up and assumptions
Debias-3 (In)stability of ascertainment mechanism. It is clear from Figure 3 that the ascertainment effects captured by δ can change rapidly and without obvious cause over time.
Contemporaneous randomised surveillance data, such as REACT or ONS CIS, allow estimation of δ. However, when predicting prevalence forward in time beyond availabil-630 ity of randomised surveillance data, we are making the implicit assumption that the ascertainment bias remains stable forwards in time, and such results should therefore be interpreted with caution.

PCR+ to infectious mapping -assumptions and caveats
For full details please see Supplementary Information-PCR positive to infectious mapping -635 method details.
Infectious-1 Pillar 1+2 positive test counts, across a four-week period, are used as an approximation to the true relative incidence over that time interval at coarse-scale level (e.g. PHE region).

30
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
Infectious-3 The infectious interval for an average individual is defined to span days 1 to 11 post infection, based on Figure 1A of Ferretti et al. [33].

SIR model -discussion, assumptions and caveats
The illustrative epidemic model we implement here has one of the simplest SIR compartmental

Gaussian approximation for δ
We approximate the cross-sectional component of the EB prior for δ using a moment-matched Gaussian approximation (see (12)). Figure 13 illustrates the suitability of this approximation for 680 PHE regions London and the North West across nine weeks.

SIR model details
We implement a DTMC SIR epidemic model based on the standard model as described in ([35], Chapter 3). As we choose ∆t to be a day/week, we allow multiple infections and recoveries in a time interval width ∆t; this requires derivation of Markov transition probabilities between all 685 states (rather than just neighbouring ones), which we do below having established some notation.

33
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021 The probability in (31) can be parameterised by the effective reproduction number, R t : We approximate (33) with a Poisson distribution as follows [45]: 11 Distribution of the number of new recoveries ∆R t The number of new recoveries, denoted ∆R t , occurring in the time interval ∆t up to time t is 11 According to Rule 2 in [45], the Poisson approximation is reasonable when both of these inequalities hold: Of the two, the first is the least likely to obtain, but is still reasonable under most circumstances. For a simple example, if we set γt = 1 and R t−1 = 1, the number of infectious individuals I t−1 > 5 is sufficient for the approximation to be reasonable.

34
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint this and so the conditional distribution for ∆I t follows from (34) and (35): Interval-to-interval transition probabilities are evaluated as where the first term in the sum at (37) follows from Assumption 2 at (28), and the second term where V t is the current number of vaccinated individuals in the population (with ∆V t ≡ V t − V t−1 ).
The total number of immune, i.e. vaccinated and/or recovered, individuals at time t (denoted R + t ) can then be represented by the recurrence This leads to the Markov conditional distribution for R + t via convolution of (35) with (38) The above treatment of immunity assumes individuals are made permanently immune immediately 715 through either vaccination or infection. It would be straightforward to relax the above formulation to allow for more sophisticated treatment of immunity, for example specifying (a) a delay in vaccine effects, (b) incomplete vaccine efficacy (e.g. in the case of novel variants), or (c) decaying immunity over time.

35
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The basic reproduction number at time t, R 0 t is related to the effective reproduction number R t by the following equation, where M is the total number of individuals and S t is the number of susceptible individuals at time t. Recall that S t ≡ M − R + t − I t where R + t is the number of immune individuals and I t is the number of infectious individuals, both of which are estimated by our DTMC SIR model. We can 725 plug in these estimates into (40) to estimate R 0 t for a given LTLA. Figure 15 plots R 0 t and R t for a selection of LTLAs.

PCR positive to infectious mapping -method details
To target the P(Infectious | PCR positive) success probability in (42), we introduce the following notation: Infectious t ≡ Individual is infectious in week t PCR+ t ≡ Individual is PCR positive from swab taken in week t (45) and proceed as follows: 12 where, at (49), we assumed conditional independence between Infectious t and PCR+ t conditional on Infected t−k . Also, at (48), we assumed that testing PCR positive implies that an individual 730 was infected at most four weeks prior to being swabbed, which is consistent with Figure 1A of [34] (data input 2 below). We import three distinct data inputs to estimate the various terms in (49). 12 We use ∧ to denote logical AND.

36
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

(which was not certified by peer review)
The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint Data input 1 -Infectious interval Figure 1A of Ferretti et al. [33] shows the estimated probability density function of the serial interval for SARS-CoV-2 transmissionwe denote this density function f Fer (d). Noting the support of this density to be approximately [1,11], we specify that an average individual is infectious between days 1 to 11. Formally we define, independently for each individual in the population, P(Infectious on dth day post-infection) := where X denotes an individual selected uniformly at random from the population. We can use this to estimate the P(Infectious t | Infected t−k ) term appearing in the numerator of (49) as follows Data input 2 -PCR positive interval and use it to estimate the term P(PCR+ t | Infected t−k ) appearing twice in (49), evaluating the following estimator for each k = 0, . . . , 3: P Hel (PCR+ | swabbed day d after becoming infected) (52) Hellewell et al. [34] helpfully provide reproducible scripts 13 and we use these to extract the pos-735 terior distribution on P Hel (PCR+ | swabbed day d after becoming infected) from their Figure 1A, whose uncertainty we propagate to estimator (52) and onwards to (49), yielding a distribution on P(Infectious t | PCR+ t ) which we take forward approximated by a moment-matched Beta distribution (at each week t) to be used as an EB conjugate prior on the success probability in (42).

Data input 3 -Pillar 1+2 incidence 740
For the purposes of adjusting the PCR positive map to changing incidence, we use the raw regional weekly positive test counts n 0:T , where we denote weeks by t = 0, . . . , T . We use this data input 13 https://github.com/cmmid/pcr-profile 37 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) which can be directly substituted for P(Infected t−k ) in top and bottom of (49) with the second term on the right of (55) cancelling between numerator and denominator, and therefore not requiring evaluation. We note that we are using raw counts to model relative incidence over a relatively short period (four weeks), which is making the assumption that the bias is relatively stable over this timeframe (see Assumption Infectious-1 in SI-PCR+ to infectious mapping -assumptions and 745 caveats).

Sensitivity analyses
Prior hyperparameters for δ The EB prior for δ depends on two hyperparameters: σ controls the variance of the white noise associated with each individual time point, while ψ controls the degree of autocorrelation from one 750 time point to the next . Figures 16 and 17 show the estimates for prevalence and R t respectively of infectious individuals using different values of these two hyperparameters. Note that in the main text, we present results using σ = 1 and ψ = 0.99.

Sensitivity and specificity of PCR tests
PCR tests are not perfect and are subject to both false positives and false negatives. In our 755 analysis, we account for imperfect testing via the false positive rate, α, and the false negative rate, β (see (8)). Figures 18 and 19 show the estimates for prevalence and R t respectively of infectious individuals using different values of these two hyperparameters. Note that in the main text, we present results using α = 0.001 and β = 0.05. 14 We use ∨ to denote logical OR.

38
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint Figure 14: EB prior on P(Infectious | PCR positive) by week and PHE region. The top panel shows raw weekly Pillar 1+2 incidence for the nine PHE regions; this is to provide intuition for the Varying incidence model in the panels below. The bottom nine panels display the prior we place on P(Infectious t | PCR+ t ), which is specific to week and region for the Varying incidence model, but is constant across weeks/regions for the Uniform incidence model (see legend in panel 2). Error bars (at left of panel for Uniform incidence; around curve for Varying incidence) represent 95% credible intervals. 39 . CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

40
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

41
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

42
. CC-BY 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted May 18, 2021. ; https://doi.org/10.1101/2021.05.17.21256818 doi: medRxiv preprint