key: cord-0477462-gun2wvov
authors: Storvik, Geir; Palomares, Alfonso Diz-Louis; Engebretsen, Solveig; Ro, Gunnar Oyvind Isaksson; Engo-Monsen, Kenth; Kristoffersen, Aja Braathen; Blasio, Birgitte Freiesleben de; Frigessi, Arnoldo
title: A sequential Monte Carlo approach to estimate a time varying reproduction number in infectious disease models: the Covid-19 case
date: 2022-01-19
journal: nan
DOI: nan
sha: 7fc56b614f9fe74bd39eb9c785d0ca23e4b377fa
doc_id: 477462
cord_uid: gun2wvov

During the first months, the Covid-19 pandemic has required most countries to implement complex sequences of non-pharmaceutical interventions, with the aim of controlling the transmission of the virus in the population. To be able to take rapid decisions, a detailed understanding of the current situation is necessary. Estimates of time-varying, instantaneous reproduction numbers represent a way to quantify the viral transmission in real time. They are often defined through a mathematical compartmental model of the epidemic, like a stochastic SEIR model, whose parameters must be estimated from multiple time series of epidemiological data. Because of very high dimensional parameter spaces (partly due to the stochasticity in the spread models) and incomplete and delayed data, inference is very challenging. We propose a state space formalisation of the model and a sequential Monte Carlo approach which allow to estimate a daily-varying reproduction number for the Covid-19 epidemic in Norway with sufficient precision, on the basis of daily hospitalisation and positive test incidences. The method is in regular use in Norway and is a powerful instrument for epidemic monitoring and management.

We propose a dynamic approach for the estimation of a time-varying or instantaneous reproduction number for a mathematical infectious disease spread model. We apply our method to the Covid-19 pandemic in Norway. Like in most other countries, the pandemic has been tackled with a combination of non-pharmaceutical interventions, from social distancing to partial lock-down, imposed or advised at various time points. Various viral variants with different characteristics have been competing in the population. Vaccination has also been gradually introduced. As a consequence of these changes, which can emerge both abruptly and smoothly, the reproduction number varies. Instantaneous estimates of reproduction numbers are useful for situational awareness. Being able to estimate such changes rapidly is important in guiding decision makers in future policy planning.

A reproduction number is precisely defined within a mathematical model of transmission. A large class of models, which has been shown to be very appropriate, is the so called SEIR model (S=susceptible, E=exposed, I=infected, R=recovered) . In particular, stochastic compartmental models are preferable when there are relatively low number of infections (Keeling and Rohani, 2011, chapter 6) . SEIR models are parametrised so that a meaningful estimate of the reproduction number (basic and effective) can be derived. Such compartmental models are based on many latent variables, in particular the number of individuals in each compartment within each region at every time point, and depend on many epidemiological parameters, including the transmission strengths which are a key component of reproduction numbers. All unknowns must be estimated from data, which in a realistic situation are scarce and incomplete. For this reason, inference is in general very difficult, because of a high dimensional parameter space, rather flat likelihoods or posterior distributions and, often, weak identifiability (De Angelis et al., 2015) .

Moreover, data which carry information about the transmissibility of the virus appear with an inevitable delay. In this paper we use two data sources, the daily number of hospitalised Covid-19 patients and the daily number of positive laboratory-confirmed RNA test cases. Both these data sources carry information about the transmission of the virus in the society, but with a random time delay from transmission. Therefore, the uncertainty of the estimates of a daily reproduction number increases in the last period of data. All this is particularly challenging during an emerging epidemic. For these reasons, instantaneous reproduction numbers are rarely assumed in SEIR models. In this paper we propose and test a sequential Monte Carlo (SMC) approach to inference which allows the efficient estimation of instantaneous reproduction numbers for the Covid-19 pandemic from actual data.

While our SMC approach is generic, we implement it on top of a stochastic SEIR model which we developed for the Covid-19 epidemic (Engebretsen et al., 2021) and which is in regular use by the Norwegian health authorities. This model assumes a spatial scale resolved on county level, and uses mobile phone mobility data for the geographical spread of the virus in temporal steps of six hours. In Engebretsen et al. (2021) , the transmissibility parameter, which represents the probability of transmission upon a contact times the contact rate in the population, is assumed to be constant in time, and is only changed at designed change-points. Inference in Engebretsen et al. (2021) is performed by a version of Approximate Bayesian Computation (ABC). MCMC convergence was essentially impossible to reach, because of the difficulty to design parameter perturbations which would lead, through the stochastic SEIR, to minor and controlled changes of the posterior distribution. It is very difficult to use ABC when the number of parameters is large, as is the case when including daily varying reproduction numbers.

In this paper we propose to perform Bayesian inference in combination with SMC. Static parameters related to the dynamic process for the reproduction number are estimated sequentially through SMC using sufficient statistics (Fearnhead, 2002; Storvik, 2002) . Application of SMC methods is challenging because the latent processes are of high dimensions, the SEIR model is only available as a computer algorithm, and data are very limited. A further very important practical aspect when working in real time during the pandemic, comes from the continuous need to improve the model, to change epidemiological assumptions, to use different data registries, and to improve the computational efficiency of algorithms. For these reasons, it has been very important to develop an SMC which is very flexible, so that changes in model specifications are easy to implement. This paper shows how a careful design of combination of model and algorithm reaches these aims and leads to a good fit to both data sources. We produce an analysis of the Norwegian pandemic, quantifying how interventions impacted the transmissibility and showing how our estimated instantaneous reproduction numbers are capturing changes rapidly enough. Based on the dynamic model, we perform future predictions of the situation for the next three weeks. We discuss the quantification of uncertainty in forecasts, which is of paramount importance for decision making.

There is an important literature on time varying reproduction numbers applied in various epidemics, see for example Cauchemez et al. (2006) ; Viboud et al. (2018) . A Bayesian framework for estimating time-varying reproduction numbers was proposed by Cauchemez et al. (2006) , and applied to the SARS epidemic. In Cori et al. (2013) a time varying reproduction number is defined as the ratio of the number C(t) of infected in day t over k w(k)C(T − k), where w(·) is the distribution of the generation time of Covid-19, which is often set to be a gamma distribution with mean and standard deviation estimated from specific studies (Ferretti et al., 2020) . There is a very successful R-package implementing this Bayesian method, called EpiEstim (Cori et al., 2021) . We include a comparison between EpiEstim and our approach based on only one of the data sources.

Several papers have applied MCMC algorithms for estimation of parameters in compartmental models (e.g Gibson and Renshaw, 1998; O'Neill and Roberts, 1999; O'Neill et al., 2000) . The recent paper Birrell et al. (2020) studies various SMC approaches in a SEIR model for influenza. It demonstrates the superiority of SMC methodology compared to MCMC for such dynamical models. The paper is also useful as a reference to SMC in epidemic modelling. Our work shares many similarities to this approach, including the use of several sources of data. A difference is our use of a dynamic model for the reproduction numbers. Also, a stochastic delay between infection and observation time is included in our setting. Another general inference framework is implemented in the R-package pomp (King et al., 2016) . The package contains multiple different implementations of estimation procedures, including Sequential Monte Carlo, for inference for partially observed Markov process models. There are several examples of applications of the package to epidemic inference, see for example (King et al., 2016; Stocks et al., 2020; Blackwood et al., 2013) .

The outline of the rest of the paper is as follows. In Section 2 the context and the data are described. Section 3 describes the full model, formulated as a state space model. In Section 4 we discuss how SMC algorithms can be applied for the inferential problem, including estimation of several static parameters. A simulation study and experimental results are reported in Section 5 specifically for the Norwegian Covid-19 pandemic. Additional results, including sensitivity analysis, are collected in the supplementary material. We conclude the paper by a summary and discussion. In the supplementary material details about experimental settings, algorithmic specifications. Data and code availability are available on a GitHub repository.

We start by setting the scene of the inferential task. The core is an existing model of the epidemic which has as input a set of parameters and variables, including daily reproduction numbers, and as output a series of time series of infection incidence. In our case, the model is a stochastic compartmental SEIR-type model that produces numbers of susceptible, exposed, pre-symptomatic, symptomatic and asymptomatic infectious and recovered at every time point. We also keep track of the disease incidence. We use two data time series to inform the SEIR model: the daily number of new hospital admissions of Covid-19 patients, and the daily number of laboratory-confirmed positive PCR tested cases. In order to exploit these data, we furthermore model the process of hospitalisation and testing of the SEIR output, in particular of the daily incidence of infected. The inferential task is to make inference on the input parameters.

The hospitalisation data contain admission to all hospitals in Norway and of all patients who are diagnosed with Covid-19 as the main cause. Admission on a certain day informs us of a transmission event that has occurred some days before. This time gap can differ between individuals. We make several assumptions on various time lags, as specified in supplementary material. For example, the number of days between symptom onset and admission to hospital is estimated to be negative binomial distributed with parameters estimated in a separate study of the Norwegian Covid-19 registry (Whittaker et al., 2021) . On average, for a patient being hospitalised, the time gap between infection and hospitalisation is estimated to be approximately 14 days. For concreteness, in this paper we assumed the distributions of the various time lags to be given.

The second data set is the time series of daily number of positive PCR tests. Again, there is a time gap between onset of symptoms and testing, which we estimate through a fixed distribution of delay with mean about 4 days. The reason is that it is important for inference that the two data sets are as aligned as possible.

We use two additional data sets, which enter the SEIR model as input variables: the daily number of positive cases who have been tested in Norway but infected outside of Norway, so called imported cases; and the total number of PCR tests made in Norway, as a surrogate of the effort made to detect positive cases. We start with the population of Norway, distributed according to the national census in the eleven counties (see www.ssb.no/statbank/table/12871/). Like in Engebretsen et al. (2021) , we seed our model continuously with positive cases imported from abroad on the day of recorded symptom onset or, if not available, when detected by testing. Imported cases that are hospitalised are not counted in the time series of hospitalisations, because they do not inform the model about the transmissibility of the virus in Norway. Because not all imported cases are likely to be discovered, we assume that each imported case stands for an unknown number of further undetected imported cases. We model this latent import with an additional Poisson distributed number of cases per observed imported case, with Poisson mean estimated from the data during calibration. We call this mean the amplification factor.

A final aspect, which is not central in this paper but that we mention for completeness, is that we use a geographical SEIR model on county (regional) level, so that the various compartments are geographically defined. Individuals are moved at random between the eleven counties of Norway using a mobility matrix, which is obtained every six hours from the movements of mobile phones, as explained in Engebretsen et al. (2021) . All parameters in the model are however shared between counties. Even if hospitalisation and test data are available at county level, in this paper we use only nationally aggregated data, because of the heterogeneity in the population size among the regions. A very important further aspect is the need to obtain inferential results as rapidly as possible, in at most a few hours, so to be able to publish results quickly just after release of the data update.

Let y t be the vector of hospitalisation and test data on day t. Let s t be the output vector of compartmental variables in the SEIR model at time t, for example the number of individuals in each county who are infected and symptomatic. Here we consider the model generically as an algorithm which outputs the compartmental variables s t at each time point t. R t is the unknown reproduction number at time t. We consider the following state space model:

for the reproduction number, (1a) s t ∼P s (R t , s t−1 ; ds t );

representing the SEIR process, (1b) y t ∼H y|s (s 1:t ; dy t );

for the hospital and test data. (1c)

To simplify notation, we do not include the dependence of the models on the set of static parameters θ. The distribution P R needs to be available analytically and easy to sample from. The distribution P s is assumed to be only available through a computer algorithm and we are only able to simulate from this distribution. In certain situations, this distribution can be available as a huge and complex Markov process. However this is not often the case, for example because of the complexity of the code or because of the lack of availability of sensitive data, like the mobility matrices in our case. The dimension of s t is large while R t is low-dimensional. In this work we consider a common scalar R t for all counties. Note that the data y t depend on the whole history s 1:t making y t only weakly informative about (R t , s t ). This is due to the fact that there is a random delay from transmission to being tested and possibly hospitalised. A graphical representation of the model is given in the left panel of Figure 1 . Our aim is to construct an efficient SMC method for the computation of p(R t , s t |y 1:T ) and we are interested in estimating the current status (t = T ), in smoothing (t < T ) and in forecasting t > T . The stochastic process {R t } is assumed to be Markov. We suggest three alternative prior models for P R . Let ε t ∼ N (0, σ 2 R ) be independent from all other variables. The first model is a random walk on log-scale, the second model extends this to an AR(1) structure. In the third

y t−1 y t y t+1 Figure 1 : Graphical representation of model (1) (left) and its reformulation (right). The hyperparameters θ are not included and can influence all conditional distributions. The dashed arrows illustrate the dependence due to the random delays between infections and testing and possibly hospitalisation.

model R t is assumed to be piece-wise constant, with jumps occurring at random.

For all models, we assume a constant start: Until a given time point D, R t is constant and equal to R 0 . The date D is set to the day when the social-distancing implementation started. Before D the reproduction number was not likely to change significantly. In this paper we set D equal to March 8, 2020, when teleworking started in many companies and universities in Norway. The Norwegian government announced the first package of interventions on March 12 2020. The models above describe the dynamics after D. All the static parameters in the three models must be estimated. The proposed method is not dependent on a specific model (1b), and would work with any epidemic model. In this work we use the particular SEIR model described in Engebretsen et al. (2019 , and whose algorithm is available though the spread package . In particular, this SEIR model has six compartments in each region (county) i: susceptible (S), exposed and not infectious (E 1 ), presymptomatic and infectious (E 2 ), infectious symptomatic (I), infectious asymptomatic (I a ), and recovered (R). The dynamics is described by the following equations

where δ t is 6 hours. We assume random mixing in each county in each 6 hour period, and move individuals between counties at the end of each such period according to a mobility matrix. In Engebretsen et al. (2021) mobile phone data are used to estimate such mobility matrices. These matrices report the number of individuals moving from county A to county B during each period, which we pick at random among the currently present in A, but favour the residents of B to return to B, to capture commuting. This rules makes the computational complexity of the SEIR model quadratic in the number of counties, due to the need for storage of both current visited location and residence of individuals. The reproduction number R t is related to β t through the equation

Finally we describe the likelihood model for the data (1c). A main difficulty is the link between the latent process {s t } and the observation process {y t } because of the unknown stochastic delays between infection and observation time, making the computation of H y|s (s 1:t ; dy t ) hard. We introduce an auxiliary process z t = (z H t,0:L H , z T t,0:L T ) with two components, one dedicated to the hospitalisation and the other to the test data. The auxiliary variable z H t,v is defined as the number of individuals who are infected at time t and hospitalised v days later. The time lag v is assumed to vary in {0, ..., L H }, for some appropriate L H . Similarly z T t,v is the number of individuals who are infected at time t and tested positive v days later. We rewrite (1c) as

where G(z t−1 , s t ; dz t ) is a Markov transition distribution assumed to be easy to simulate from. Now H y|z (z 1:t ; dy t ) is easy to compute. In more details, define y H t to be the number of daily Covid-19 admissions to hospital. An individual who is infected at time t is hospitalised with probability p H u at time t + u. The time lag u from infection until hospitalisation is assumed to follow a discrete distribution on the integers {0, 1, ..., L H }. Let ρ H u be the probability of delay u. To structure the model further, we attach to each infected individual its potential time lag until hospitalisation. So z H t,u is the number of individuals infected at time t and who could possibly be hospitalised u time-units later. Let I t be the total number of infected individuals at time t, which is available through the SEIR model as a component of s t . We can formulate the model as

where we here consider a Beta-Binomial distribution for a patient being hospitalised to take into account variability in hospitalisation between regions. The β H t parameter is specified indirectly through a time-varying probability p H t such that β H t = α(1 − p H t )/p H t where p H t is predefined using the age-structure of the individuals having tested positive. By storing and sequentially updating the quantities L u=0 z H t−u,u as well, we obtain a first order Markovian state space structure as illustrated in the right panel of Figure 1 .

Regarding the test data, as in Engebretsen et al. (2021) we consider the probability that an infected case is tested by means of a PCR test. We ignore that tests can lead to false positive responses. The logit of this probability is assumed to be linear in the total number of daily tests N T t in day t in addition to a time independent intercept. We write the detection probability ρ T t at time t as

The time lag between infection and testing is assumed to follow a discrete distribution on {0, 1, ..., L T } for an appropriate L T . The approach for handling this delay is exactly as for the hospital incidence, with p T t now playing the role of p H t . We introduce a new set of auxiliary variables for the test data {z T t,u }, similarly to the ones introduced for the hospitalisation data. Defining z t = (z H t,0:L H , z T t,0:L T ), we are within the model formulation (3).

Our aim is to perform inference on the whole set of latent variables x 1:t = (x 1 , ..., x t ) as well as on static hyper-parameters θ at each time-point t, by means of the posterior Algorithm 1 Auxiliary SMC with resampling at each time step. Operations involving index b must be performed for b = 1, ..., B. Here P x t denotes the transition distribution for x t , while Q t is the proposal distribution for x t . The indices A 1:B t defines the ancestral particles at time t after resampling.

Estimate of p(y 0 ) 4: for t = 1 to T do 5:

Estimate of p(y t |y t−1 ) 9: end for distribution p(x 1:t , θ|y 1:t ). A description of the SMC algorithm with resampling at each step is given in Algorithm 1. See also Chopin and Papaspiliopoulos (2020) , which contains more general algorithms. Although in this paper we focus on inference for R t , also s t will be of interest. In our setting, the main computational burden is the sampling from P s (R t , s t−1 ; ds t ) which has been parallelised in our implementation. For resampling, residual resampling (Liu and Chen, 1998) has been applied. However, the resampling step is both hard to parallelise and requires message passing, resulting in that a too high number of cores can decrease performance.

The SMC algorithm is by design sequential so that by storing values of x t obtained at the previous day, updates can easily be performed as new data arrive. A main challenge here is that the state x t at time t heavily depends on future observations y t+h because of the delay in hospitalisation. Although the reformulated model reduces immediate dependence, there are still strong correlations backwards, as illustrated in Figure 1 . There are clearly possibilities to develop more efficient proposal distributions, despite the availability of the SEIR model only as a computer algorithm. Because the main purpose of this paper is to investigate the utility of SMC methods for the estimation of daily varying reproduction numbers, we use only simple bootstrapping proposals. We do however take into account the delay aspect by using fixed lag smoothing, using data ahead of current time, that is

Fixed lag smoothing with a lag of l t = 24 days has been used in the Covid-19 runs in this paper. At the end of the time-series, l t = min{T − t, l} is used. The estimates on the last days will be more uncertain.

Algorithm 1 assumes that the parameters θ are known. Now we describe Bayesian inference for some of the static parameters. We denote by θ R the set of parameters in the model for {R t }, θ s the ones in the {s t } process, and θ y the parameters appearing in the data model. As mentioned, some of these parameters are fixed based on other data sources, and here for simplicity we do not propagate their uncertainty. For other parameters, sequential updates of their estimates are desirable. In principle, all parameters θ could be included as part of the state vector x t where the propagation of these static components just keeps them fixed. However, repeated resampling will quickly give degenerate samples for these parameters.

A review of parameter estimation in SMC is given in Kantas applications, but require repeated runs of the SMC routine. Although much smaller number of particles can be applied in such settings, some experiments with our models indicate that at least 250 particles are necessary, in which case one run uses about 10 minutes using 4 cores on Linux server and more cores did not help much in this case. Some experiments with an implementation of the Particle Metropolis-Hastings algorithm, based on the pseudo-marginal method by Andrieu et al. (2009) , is reported in supplementary section D. We will however focus on online methods for parameter estimation here. In cases where sufficient statistics v t (x 1:t ) for the parameters are available, the SMC algorithm can be easily updated to target p(x 1:t , v t (x 1:t )|y 1:t ) instead (Fearnhead, 2002; Storvik, 2002) . This is the case for the θ R parameters. Simulations of θ R at each time point can then be obtained from p(θ R |v t (x 1:t )) which then again can be used to obtain new samples of x t+1 (step 6 in Algorithm 1). A crucial step is that v t (x 1:t ) can be recursively updated. Section A in the supplementary material gives details on how θ R can be updated by this approach. Note that these methods can suffer from the degeneracy problem (Andrieu et al., 2005) . In supplementary section D we validate the parameter estimates obtained by this procedure both through comparisons between different runs and by using the samples obtained by the Particle Metropolis-Hastings algorithm.

For the analysis of the Norwegian Covid-19 case, we have used the 11 counties as our spatial scale. Mobility data and imported cases are used at regional scale while hospital incidence data and test data are used at national scale. Also the {R t } prior process is assumed to be common for all regions. The hospital incidence data are from the Norwegian national Beredt-C19 registry and the Norwegian Intensive Care and Pandemic Registry, and the test data are from the Norwegian Surveillance System for Communicable Diseases (MSIS registry). The reproduction number is assumed constant until March 7 2020. The results are based on data up to July 1, 2021 with test data included from August 1, 2020, when testing capacity in Norway was scaled as needed, and after which testing criteria had become rather stable. The number of parameters is large: the dimension of s t is 3157. In addition we have the reproduction numbers R t and the auxiliary variables z t 's. The prior distributions assumed for the parameters involved in the dynamics of R t are given in the first column of Table 1 . When running the model, we did not estimate the parameters in the SEIR model (2) nor the parameters π 0 and π 1 in (5), because they were estimated separately as described in Engebretsen et al. (2021) . Also the parameters related to observations in equation (4) were pre-estimated through other data sources. All further details on the model are given in supplementary material section B. Each run is based on B = 20 000 particles. One run of the 500 days considered here, using a linux server with 128 cores, took approximately 5 hours, which is appropriate for practical real-time purposes. Figure 2 shows our estimates of R t , with a 7-precedent-days moving-average smoothing (daily estimates are given in the supplementary material Figure S3 ), using the three considered models for the reproduction number. Posterior medians and symmetric 95% credibility intervals of the estimated parameters are reported in the second column of Table 1. The autoregressive model for R t was simplified by fixing µ = 0: when we estimated µ, we obtained an estimate very close to zero, but also some difficulties in the estimation of σ 2 . In addition, the marginal log-likelihood (MLIK) was slightly higher. We therefore opted for the simpler autoregressive model.

In Figure 2 we see that the different models lead to rather similar estimates for the time dynamics of {R t }. The autoregressive model (central panel) seems to have more uncertainty. On the other hand, this model is to be the preferred in terms of MLIK (Table 1 , last rows). Note that all models capture quickly the dramatic reduction in transmission in mid March 2020 when Norway had the first lockdown, which was implemented between March 9 and March 14 2020. In the second half of April 2020, society was re-opened, including schools and kindergarten, the 2-meters distance rule was reduced to 1-meter: {R t } appears to increase to around 1. A new peak appears in the end of July and beginning of August 2020, in correspondence with the end of the Norwegian vacations, the returning of Norwegians from abroad and the arrival in Norway of tourists. In particular there have been several isolated clusters, for example in two cruise ships and the return of students to university campuses. These clusters were well controlled by contract tracing, and the reproduction number dropped rapidly below 1. The autumn 2020 was also characterised by a series of larger local outbreaks, which were rapidly controlled. The growth of R t starting in October 2020 was also due to local outbreaks, and the incoming winterseason which affects the viral transmission. The number of cases was so large that contact tracing became less effective. The Norwegian government imposed a second national intervention, first in the end of October 2020 and then again in early November 2020. This reduced the reproduction number again below 1, in two pushes, where we can see that the second strengthening of the intervention was in fact needed for this purpose. We can see that the autumn 2020 interventions allowed the reproduction number to fall from approximately 1.5 to below 1 in about three weeks. Interestingly, we see a peak of R t just around Christmas 2020, in connection with vacation travel and intensive shopping. January 2021 marks the arrival of the alpha variant of the virus, which was more transmissible and which increased the risk of hospitalisation. The alpha variant was predominant in Norway at the end of March 2021, and we see R t just below 1.5 again. This increase happened despite at the end of January Norway introduced the strictest lockdown rules during the whole pandemic so far, including essential closure of all borders, and vaccination of the elderly started. During March and April 2021, R t started to decrease again, to remain below or around 1 until the middle of May 2021. Governmental interventions were reduced from May 2021, and the reproduction number stabilised around 1. At the end June 2021, approximately 50% of the adult Norwegian population had been vaccinated at least once, and approximately 30% twice. The effective reproduction number R t reflects the effect of vaccination which is included in the SEIR model. It is remarkable how the estimated reproduction number quantify the history of the epidemic so precisely. Another feature is the cyclic behaviour of R t , with a drop following an increase. The AR prior model on log scale also attracts towards 1. In Norway this is expected because of the rapid intervention strategy of the government (named "control") whenever R t was growing rapidly above 1, and the rapid reopening when the epidemic appeared under control. Local outbreaks were frequent, also visible in the raw data, but they were rapidly controlled by appropriate contact tracing and other successful local interventions.

The red lines in Figure 2 give 50% centred credibility intervals based only on hospitalisation data. We investigate in this way the value of the test data as a second source of information. Because we used test data only from August 1 2020, the estimates are essentially identical until then. After August 1 2020, the estimates based only on hospital data are smoother, indicating that the test data contain information about transmission (and its change) that is not transferred to the hospitalisation. One reason for this is that the younger generations have been infected in the autumn more than the elderly ones, who are most at risk for hospitalisation. The test data also contribute to a more precise estimate of the daily prevalence of infected in Norway. We also observe some misalignment in time between the two estimates, probably because the time lags are not stationary, while we assume them constant during the whole epidemic.

The last green horizontal line in Figure 2 corresponds to the last day with data, July 1, 2021, u=0 R t−u based on the random walk (RW, upper), autoregressive (AR, middle) and piece-wise constant model (PC, lower). We used hospital data from March 8 2020 to July 1 2021, test data from August 1 2020 to July 1 2021. We used B = 20 000 particles. Corresponding not-averaged daily estimates are given in Figure C. 2. The first vertical green line corresponds to the date at which test data is included in the analysis, the second green vertical line corresponds to the last date with observations, after that Bayesian prediction is performed. Blue shadow bands indicate posterior uncertainty quantiles. The red curves correspond to 50% median-centred credibility intervals only based on hospitalisation data. after which we perform a three weeks prediction, by simply running the model forward in time.

The predictions for the three models have similar means but differ in uncertainty. Both the random walk and the piece-wise constant model are non-stationary, which explains the high increase in uncertainty after the last observation point. On the other hand, the autoregressive model shows a more stable prediction performance, as expected. Figure 3 gives predicted number of new daily infected cases. Before the last green dashed line these credibility intervals give predictions based on the observed data, after this line, the predictions are obtained by running the SEIR model forward three weeks in time using the predicted reproduction numbers. When predicting forward in time, mobility matrices, imported cases and total number of tests are needed as input. We here re-use the values in the previous 21 days in the forecasts. We see the three waves which hit Norway in March/April 2020, from November 2020 to January 2021 and in March 2021. The number of cases is estimated around 1000 per day in the peaks (1300 in the third wave). The number of cases would in fact grow during July 2021, as here correctly predicted. Table 1 provides the estimates of the parameters in the three dynamic models of {R t }. How these estimates are learned over time is shown in Figures C.3-C.5 in the supplementary material. The plots also allow a comparison of the prior (which is the first time point to the left) and the final posterior estimate (last time point on the right). The estimates stabilise nicely. We report about some limited validation of the parameter estimates in section D. Note that the variance of the random walk dynamics is estimated to be smaller than the one of the autoregressive model, as is also clear in Figure 2 . On the other hand, the variance related to the piece-wise constant model is considerably smaller even if the variability seems to be smaller in Figure 2 . This is due the fact that for most time points there are no discontinuities while when changes occur there may be large discontinuities.

In Figure 4 we use the estimated parameters, including the instantaneous reproduction numbers, to simulate the daily hospitalisation incidence and the daily number of positive tests. We propagate uncertainty and produce probabilistic estimates, which we compare with the actual data. These plots show that we are able to fit both data sources well. We note the weekly structure in the test data. These plots also show three weeks ahead forecasts. The superimposition of the actual data of these three weeks, which were not used in the analysis, show that predictions were good.

There are several parameters that are fixed in the current implementation:

• Hyper-parameters in the prior model for {R t }, see Table 1 .

• Parameters related to the observations: the age-dependent probabilities of being hospitalised and detected positive by testing, and the parameters describing the distribution for delay from infection to testing or hospitalisation.

• Several parameters inside the SEIR model.

• Number of particles used in the SMC runs.

We did not perform a systematic sensitivity analysis, but focus on the SMC design parameters. The difference in the estimated marginal likelihoods is small compared to the differences between models (about 5.0). This figure (and other similar ones, not included) shows that the results are quite stable and that B = 2 000 might have been sufficient for estimation of {R t }. However, if interest is also in the latent structure s t or marginal likelihood values, more particles have shown to be necessary.

EpiEstim (Cori et al., 2013) is a popular method for estimation of the reproduction number R t based on case incidence as well as imported cases. It does however not allow for the incorporation of multiple data, so that hospital data are not taken into account, nor mobility. In Figure 5 we compare estimates of R t obtained by our SMC approach, though only using test data and this time from the start of the epidemic, with the results obtained by using the EpiEstim package in R (Cori et al., 2021, version 2.2.4) . For EpiEstim we assumed a serial interval (the time between the onset of symptoms of a primary case and the onset of symptoms of secondary cases) with mean 7.5 days and a standard deviation of 3 days. We see a very good agreement between the two estimated curves in general. However, the confidence bands are much narrower for EpiEstim. There are some differences between the estimates in March 2020 and during the summer 2020 and in November 2020.

Based on an estimate {R * t } obtained by the AR model from the real data, we simulated (s t , z t , y t ) from the model for all t. Twenty independent data sets were generated and for each case the simulated set {y t } was used for estimating {R t }. Figure 6 shows the results from one of these simulated data sets, demonstrating that the method is able to capture the main structure quite u=0 R t−u using B = 20 000 (blue) based on one simulated data set. Same setup as Figure 2 , hospitalisation data from March 12 2020 to July 1 2021, test data from August 1 2020 to July 1 2021. The red curve corresponds to the assumed true {R * t } process.

well, although in periods where the true R * t is considerably unstable (July/August 2020) estimates are worse. Figure C. 6 in the supplementary material shows results based on all 20 data sets, demonstrating the strong information content in such type of data about the {R t } process.

A time-varying transmissibility allows to quantify the effects of interventions and changes in the behaviour of people, in real-time. This is of key importance to policy makers, as the interventions often have immense societal and health costs. Understanding, while an epidemic develops, whether the implemented interventions are sufficient or not, or if interventions could be lifted, is essential. The Norwegian government's strategy was to control the epidemic, and this was achieved by multiple national and local non-pharmaceutical interventions, which are reflected in the temporal variations of the reproduction number. Our approach is the first which allows to monitor a dailyvarying reproduction number when using a complex compartmental model informed by multiple streams of data. The fact that our estimates of R t react rapidly to changes in the test data, means that the situation is captured only with a delay given by the generation time of the disease under study and the time gap between transmission and testing. For Covid-19, this amounts to about a week, because of a generation time of about 5 days and a delay between transmission and test of about 2 days. Picking up an exponential growth (R t > 1) before the epidemic grows out of control is essential for surveillance. The possibility of our method to validate the efficacy of contact tracing, to lead back R t to below 1, or not, is also very important.

We have shown how daily reproduction numbers and the latent compartment-wise populations in an SEIR model can be put into a state space model, so that an SMC technique for inference can be used. Obtaining unbiased estimates of the marginal likelihoods also makes it possible to do parameter estimation within a particle MCMC framework, although more work is needed here to make this computationally efficient. Compared to a parallel effort using Approximate Bayesian Computation (ABC) (Engebretsen et al., 2021) , the SMC approach is much faster and also easier to modify with respect to model changes, confirming the findings in Birrell et al. (2020) . Our implementation is modular, does not depend on the specific epidemic model (here SEIR), so that alternatives can easily be tested.

So far only simple bootstrap filters have been applied. This can be improved by utilising more efficient proposals, alternative algorithms such as the resample-move algorithm (Gilks and Berzuini, 2001) which was used in Birrell et al. (2020) or the recent promising ideas of iterated auxiliary filters and twisted models (Guarniero et al., 2017; Heng et al., 2020) . We expect these approaches to be very useful when we expand the SEIR model to have different reproductions numbers in each of the eleven Norwegian counties.

We compared three simple dynamics for R t and found that the autoregressive model was slightly better than the others. More work needs to be done here to compare the models in terms of prediction. In our approach it is easy to predict the future hospitalisation incidence and the number of positive cases tested. Note however that such simple dynamic models for {R t } are mostly suitable for now-casting and for short term forecasting, because of the lack of stationary that interventions and feedbacks imply. There are then interesting questions on how to use our probabilistic prediction of the time varying transmission strength to propagate uncertainty, in the context of variable planned interventions.

Estimation of (static) parameters is a challenging task in SMC. Several parameters related to the SEIR model, as well as parameters related to the observation processes, were pre-estimated based on external data sources. In (Engebretsen et al., 2021) , we used a version of ABC. Parameters in the {R t } dynamic process were estimated online based on the procedure of sufficient statistics, a method that can lead to degeneracy. We have also tested out the particle Metropolis-Hastings algorithm by Andrieu et al. (2010) . However convergence was slow and challenging, because of the computational cost of running the SMC algorithm even for a small number of particles. We validated our estimates in the supplementary material, and find that our online estimation procedure worked reasonably well for the given models and data. However, some experiments with the AR model when also including estimation of the parameter µ, did cause degeneracy problems as have been reported in in Andrieu et al. (2005) . A more efficient SMC algorithm might be better in utilising the potential of such algorithms to estimate static parameters. However, our estimates of {R t } appear to be relatively robust with respect to changes in these parameters. Other parameters related to the SEIR model and the observation processes can be more important.

Communicating uncertainty of estimates and the effect of stochastic and uncertain time lags from data back to infections, is a major challenge. Our current strategy has been to report estimates of the reproduction numbers one week back in time as the most reliable estimates, and this needs to be studied further.

As pointed out by one of the reviewers, it is quite easy to extend the modelling approach (and the algorithm) to include dummy variables describing interventions made by the government. In this case, one would need to estimate a time delay between the time point in which the interventions are decided and when they effect viral transmission. This delay might change in time and be intervention specific. In our model, interventions appear in the data after a delay, and are then reflected in a change of the reproduction number. It is possible to interpret changes in the estimated reproduction numbers in the light of the interventions set in place.

Finally, it is interesting to compare the estimates of the instantaneous reproduction number produced by this SMC model, with the estimates obtained by models which keep the reproduction number constant over longer time intervals, for example over four weeks. The SMC-based R t is able to capture changes at shorter time scales significantly better, but possibly with larger uncertainty than estimates of reproduction numbers assumed constant over longer time periods, if the transmission has been stable during such periods. Comparing prediction power is a further aspects that can be compared. Results from our SMC model are currently used in the weekly reports of the Norwegian Institute of Public Health, see www.fhi.no/en/publ/2020/weekly-reportsfor-coronavirus-og-covid-19/.

The data sets analysed in this paper come from the national emergency preparedness registry for Covid-19, owned by the Norwegian Institute of Public Health. The preparedness registry is temporary and comprises data from a variety of central health registries, national clinical registries and other national administrative registries. Further information on the registry, including access to data, is available at www.fhi.no/en/id/infectious-diseases/coronavirus/ emergency-preparedness-register-for-covid-19/.

An R package called smc.covid is available at github.com/geirstorvik/smc.covid. It contain scripts for running the SMC algorithm on some data examples. Because data are sensitive, not all data sets are public. The number of individuals in hospital and the number of positive cases per day are available, but the number of imported cases and the mobility matrices are confidential. For the two latter data sets we have therefore provided some simulated data sets which resemble the true ones.

Note first that

This gives

Further, for any value of a, we have p(τ |R 1:t ) ∝ p(τ )p(R 1:t |τ, a) p(a|τ, R 1:t ) and using a =â t (a), we obtain p(τ |R 1:t ) ∝ τ ατ −1 τ 0.5t exp(−β τ τ ) exp(−0.5τ g t (â)) τ 0.5 ∝ Gamma(α τ + 0.5t, β τ + 0.5g t (â t )).

The full model includes several parameter specifications that are fixed throughout all runs. Some sensitivity analysis is provided in Section C.

• R t was set to a constant R 0 from February 17 2020 to March 8 2020. We estimate R 0 and used a N (1.192, 0.1) prior for log R 0 (giving an expectation of R 0 equal to 3.3, which is compatible with the literature).

• The amplification factor for imported cases is set to 2.8 from (Engebretsen et al., 2021) .

• In order to take into account regional variability, the number of hospitalised individuals is assumed to follow a beta-binomial distribution with expectations equal to the probabilities given in Table 2 and with a shape parameter equal to 8. The time-dynamic changes in the probabilities are based on (Engebretsen et al., 2021) .

• The distribution for the delay from infected to hospitalised is assumed to follow a negative binomial distribution with parameters (8.87, 3.40) until August 1 2020 and thereafter with parameters (7.60, 3.34), as estimated from registry data. Table 3 : On the left, estimates of marginal log-likelihoods based on B = 5 000 particles using five repetitions of the SMC algorithm. On the right, p-values based on a t-test for testing the hypothesis on equality for the marginal log-likelihoods between the models.

• The distribution of the time between infection and testing positive is assumed to follow a discrete distribution with probabilities 0.05, 0.080, 0.16, 0.4, 0.3, 0.01 for delays of 1-6 days.

• In the model (5) for the detection probability, we used π 0 = −1.017, π 1 = 0.00013. The number of detected individuals was assumed to follow a beta-binomial distribution with expectations given by ( 5) and a shape parameter equal to 8. Figure 2 ). For the red curve, we changed the prior for τ = 1/σ 2 to a Gamma(1.2, 0.14) distribution. For the blue curve, we changed the values of (π 0 , π 1 ) to (-1.175,0.000074) (estimates obtained from another model on the same data). The green curve corresponds to a change in the shape parameters in the negative binomial distribution for hospitalisation and test data from 8 to 3. In all cases the results are stable, which was also the case when we compared 50% credibility intervals (not shown). Figure C .2 shows the fixed-24-days-lag smoothed estimates of daily R t values for the three models. This plot is not averaged over each preceding week, as is done in Figure 2 in the main text.

Figures C.3-C.5 show (filtered) parameter estimates of the dynamics of {R t } for the three models considered. All parameter estimates are based on the use of sufficient statistics (Fearnhead, 2002; Storvik, 2002) in the SMC algorithm. Table 3 shows that although there is some Monte Carlo variability in the estimation of the marginal likelihoods, the differences between the models are statistically significant. Table 1 , the red curve is obtained by changing (π 0 , π 1 ) to (-1.175,0.000074), the green curve is obtained by changing the prior for τ to Gamma(1.2, 0.14) and the blue curve corresponds to changing the α parameters in the Negative-binomial distributions for hospital and test data to 4.

As described in Kantas et al. (2015) , there are many different approaches to estimate static hyperparameters, both online and offline. In the paper we focused on online methods for estimation of parameters in the {R t } process, based on the use of sufficient statistics (Fearnhead, 2002; Storvik, 2002) . As pointed out by one of the referees, this method can suffer from degeneracy and converge towards different limit distributions in different runs. Such degeneracy problems can in particular occur for long time series. We validate the estimates obtained in two ways. We concentrate on the AR process, both due to the fact that this is our best model and that it turned out to be the most problematic one, with respect to parameter estimation. Figure D .1 shows the parameter estimates obtained from two independent runs. The agreement of these runs indicates that the estimation procedure performs quite well in our case. We do however emphasise that even though the problems with degeneracy did not seem to be present for the model considered here, if we considered an AR process where also µ is estimated, indeed different runs of the SMC algorithm resulted in somewhat different parameter estimates. The problems in that case occurred on the day where test data were introduced.

An important feature of SMC methods is that unbiased estimates of marginal likelihoods are available, making the pseudo-marginal approach by Andrieu et al. (2009) possible. Given a current parameter θ with corresponding estimate of the marginal likelihoodp(y 1:T ; θ), a proposal θ ∼ q(θ|θ) is generated and an estimatep(y 1:T ;θ) is obtained by a new run of the SMC algorithm. The new proposal will then be accepted with probability min{1, p(θ)p(y 1:T ;θ)q(θ|θ) p(θ)p(y 1:T ; θ)q(θ|θ)

},

corresponding to the Particle Metropolis-Hastings Algorithm in Andrieu et al. (2010) . Practical use of this algorithm requires tuning of the proposal distribution q(θ|θ) and of the number of particles used in each SMC run. Concerning the latter, the number of particles in each SMC run needs to be reduced dramatically compared to a single run of the SMC algorithm since these runs have to be repeated many times. Although Andrieu et al. (2010) showed that the number of particles should increase with time, we have only considered a fixed number of particles in our testing. Even with many trials on the tuning parameters, we struggled to obtain convergence within reasonable time (running the algorithm for several days). However, we stored all the proposed parameter settings together with their likelihood values. We then fitted a two-dimensional GAM model to the loglikelihood values using a and σ as covariates. More than 500 samples of (a, σ) were obtained through multiple runs of the PMCMC algorithm with different starting points (running in total for several days). Figure D. 2 shows the fitted GAM curve with the red circles corresponding to the samples obtained from the SMC algorithm. The SMC samples seem to be within the higher likelihood region of the fitted curve. Note that due to lack of converge, the PMCMC algorithm has probably not been able to explore the full posterior properly. Further, the "observed" values displayed as small dots in the plot are all proposed values, not only the accepted ones, so these points should not be used to reflect the whole distribution.

Particle Markov chain Monte Carlo methods

On-line parameter estimation in general statespace models

The pseudo-marginal approach for efficient Monte Carlo computations

Efficient real-time monitoring of an emerging influenza pandemic: How feasible?

Deciphering the impacts of vaccination and immunity on pertussis epidemiology in thailand

Estimating in real time the efficacy of measures to control emerging communicable diseases

An introduction to sequential Monte Carlo

A new framework and software to estimate time-varying reproduction numbers during epidemics

EpiEstim v2.2-3: A tool to estimate time varying instantaneous reproduction number during epidemics

Four key challenges in infectious disease modelling using data from multiple sources

Time-aggregated mobile phone mobility data are sufficient for modelling influenza spread: the case of Bangladesh

A theoretical singleparameter model for urbanisation to study infectious disease spread and interventions

spread: Infectious Disease Spread Models. Norwegian Institute of Public Health

Spatial modelling of the early-phase of the COVID-19 epidemic in Norway

Markov chain Monte Carlo, sufficient statistics, and particle filters

Quantifying SARS-CoV-2 transmission suggests epidemic control with digital contact tracing

Estimating parameters in stochastic compartmental models using Markov chain methods

Following a moving target-Monte Carlo inference for dynamic Bayesian models

The iterated auxiliary particle filter

Controlled sequential Monte Carlo

On particle methods for parameter estimation in state-space models

Modeling infectious diseases in humans and animals

Statistical Inference for Partially Observed Markov Processes via the R Package pomp

Sequential Monte Carlo methods for dynamic systems

Analyses of infectious disease data from household outbreaks by Markov chain Monte Carlo methods

Bayesian inference for partially observed stochastic epidemics

Model selection and parameter estimation for dynamic epidemic models via iterated filtering: application to rotavirus in germany

Particle filters for state-space models with the presence of unknown static parameters

Simonsen, A. Vespignani, and the RAPIDD Ebola Forecasting Challenge group

Trajectories of hospitalisation for patients infected with SARS-CoV-2 variant B. 1.1. 7 in Norway

We acknowledge the support from many colleagues at the Norwegian Institute of Public Health in the collection and preparation of the data. Funding to this project was provided by the centre BigInsight (NFR project 237718), the project "A real-time analytical pipeline for preparedness, planning and response during the Covid-19 pandemic in Norway" (NFR project 312721) and the Nordic project "Data streams and mathematical modelling pipelines to support preparedness and decision making for Covid-19 and future pandemics" (Nordforsk NordicMathCovid project).

A Derivations for the dynamic models for {R t } Here we define L t = log(R t ) for simplicity of notation and neglect the first time points (before March 7 2020) for which R t is constant. Note that for each model, the sufficient statistics involved are all easily updated with increasing t.

Assume a prior p(τ ) = Gamma(α τ , β τ ) with τ = 1/σ 2 R . Given {R t }, the posterior for (φ c , τ ) becomesSimulation of L t+1 can then be performed by first simulating τ from the Gamma distribution and thereafter R t+1 from the model.

Assume a prior p(φ c , τ ) = Beta(a φ , b φ )Gamma(α τ , β τ ) with τ = 1/σ 2 R . Given {R t } and changepoint indicators C t (= 1 if there is a changepoint at time t), the posterior for (φ c , τ ) becomeswhere n c t is the number of changepoints up to timepoint t. Simulation of R t+1 can then be performed by first simulating (φ c , τ ) from the Beta-Gamma distribution and thereafter R t+1 from the model.

We consider here only the case where µ = 0. Assume p(a, τ ) =N (a 0 , [κ 0 τ ] −1 ) × Gamma(α τ , β τ )