key: cord-0151318-zjnr7vwm
authors: Altmejd, Adam; Rocklov, Joacim; Wallin, Jonas
title: Nowcasting Covid-19 statistics reported withdelay: a case-study of Sweden
date: 2020-06-11
journal: nan
DOI: nan
sha: 48b2dd6196779cc186362d01b75dbd5084ca89ab
doc_id: 151318
cord_uid: zjnr7vwm

The new corona virus disease -- COVID-2019 -- is rapidly spreading through the world. The availability of unbiased timely statistics of trends in disease events are a key to effective responses. But due to reporting delays, the most recently reported numbers are frequently underestimating of the total number of infections, hospitalizations and deaths creating an illusion of a downward trend. Here we describe a statistical methodology for predicting true daily quantities and their uncertainty, estimated using historical reporting delays. The methodology takes into account the observed distribution pattern of the lag. It is derived from the removal method, a well-established estimation framework in the field of ecology.

The new corona virus pandemic is affecting societies all around the world. As countries are challenged to control and fight back, they are in need of timely, unbiased, data for monitoring trends and making fast and well-informed decisions (Nature, 2020) . Official statistics are usually reported with long delay after thorough verification, but in the midst of a deadly pandemic, real time data is of critical importance for policymakers (Jajosky and Groseclose, 2004) . The latest data are often not finalized, but change as new information is reported. In fact, reporting delays make the most recent days have the least cases accounted for, producing a dangerous illusion of an always improving outlook.

Still, these unfinished statistics offer crucial information. If the pandemic is indeed slowing, we should not wait for the data to be finalized before using it. Rather, we argue that actual case counts and deaths should be nowcasted to account for reporting delay, thus allowing policymakers to use the latest numbers availiable without beinig misled by reporting bias.

Such predictions provide an additional feature that is perhaps even more important. They explicitly model the uncertainty about these unknown quantities, ensuring that all users of these data have the same view of the current state of the epidemic.

In this paper we describe a statistical methodology for nowcasting the epidemic statistics, such as hospitalizations or deaths, and their degrees of uncertainty, based on the daily reported event frequency and the observed distribution pattern of reporting delays. The prediction model is building on methods developed in ecology, referred to as the "removal method" (Pollock, 1991) .

To help motivate why such forecasting is needed, we now turn to the case of Sweden. The model is flexible by design, however, and could easily be applied to other countries as well.

The Swedish Public Health Agency updates the COVID-19 statistics daily 1 . During a press conference, they present updates on the number of deaths, admissions to hospitals and intensive care, as well as case counts.

One of the reasons for following these indicators is to enable public health professionals and the public to observe the evolving patterns of the epidemic (Anderson et al., 2020) . In relation to policy, it is of specific interest to understand if the growth rates changes, which could indicate the need for a policy response. However, in each daily report only a proportion of the number of recent deaths is yet known, and this bias produces the illusion of a downward trend.

The death counts suffer from the longest reporting delay. In their daily press conference, the Swedish Public Health Agency warns for this by stopping the reported 7-day moving average trend line 10 days before the latest date. But not only are deaths often reported far further back than 10 days, a bar plot still shows the latest information, creating a sense of a downward trend. In fact, this might be the reason why the number of daily deaths have been underestimated repeatedly. At the peak, deaths were initially believed to level out at around 60 per day, but after all cases had been reported more than two weeks later, the actual number was close to 120 (Öhman and Gagliano, 2020) .

We propose to use the removal method, developed in animal management (Pollock, 1991) , to present an estimate of the actual frequencies at a given day and their uncertainty. The method has a long history dating back at least to the 1930s (Leslie and Davis, 1939) . However, the first refined mathematical treatment of the method is credited to Moran (1951) , more modern derivatives exits today (Matechou et al., 2016) . It is a commonly applied method today when analyzing age cohorts in fishery and wildlife management.

The removal method that has three major advantages over simply reporting moving averages:

• it does not relay any previous trend in the data,

• we can generate prediction intervals for the uncertainty about daily true frequencies,

1 The data is published on https://www.folkhalsomyndigheten. se/smittskydd-beredskap/utbrott/aktuella-utbrott/covid-19/ bekraftade-fall-i-sverige/.

• the uncertainty estimates can be carried over to epidemiological models to help create more realistic models.

A classic example where the method proposed to solve this problem has been used is in estimating statistics of trapping a closed population of animals (Pollock, 1991) . Each day the trapped animals are collected, and kept, and if there is no immigration the number of trapped animals the following days will, on average, decline. This pattern of declining number of trapped animals allows one to draw inference of the underlying population size. Here we replace the animal population with the true number of deaths on a given day. Instead of traps we have the new reports of COVID-19 events. As the number of new reported deaths for a given day declines, we can draw inference on how many actually died that day. If we assume that the reporting structure is constant over time we can after a while quickly get good estimate of the actual number.

Suppose for example that on day one, 4 individuals are reported dead for that day. On the second day, 10 deaths are recorded for day two. Then, with no further information, it is reasonable to assume that more people died on day two. If the proportion reported on the first day is 3%, the actual number of deaths would be 133 for day one and 333 for day two.

If additionally, 60 deaths are reported during the second day to have happened during day one, and on the third day, only 40 are reported for day two, we now have conflicting information. From the first-day reports it seemed like more people had died during day two, but the second day-reports gave the opposite indication. The model we propose systematically deals with such data, and handles many other sources of systematic variation in reporting delay. In fact, the Swedish reporting lag follows a calendar pattern. The number of events reported during weekends is much smaller. To account for this, we allow the estimated proportions of daily reported cases to follow a probability distribution taking into consideration what type of day it is.

We propose a Bayesian version of the removal model that assumes an overdispersed binomial distribution for the daily observations of deaths in Sweden in COVID-19. We then calculate the posterior distribution, prediction median and 95% prediction intervals of the expected deaths from the reported deaths on each specific day. The method and algorithm is thoroughly described in the Supplementary Information.

To get accurate estimates we apply two institution-specific corrections. First, we only count workdays as constituting reporting delay, as very few deaths are reported during weekends. Second, we apply a constant bias correction to account for the fact that Swedish deaths come from two distinct populations with different trends: deaths in hospitals, and in elderly care.

In Figure 1 we apply the model to the latest statistics from Sweden. The graph shows reported and predicted deaths (with uncertainty intervals) as bars, and a dashed line plots the 7-day (centered) moving average. A version without predictions is used in the Public Health Agencyś daily press briefings. As expected, the model provides estimates of actual deaths considerably above the reported number of deaths. Not how the model predicts additional deaths above the moving average line. 

To judge whether or not the model is accurate we need to compare it to a benchmark. The moving average of reported deaths is not useful, since it is biased for deaths that occurred within the last week. Instead, we create a benchmark prediction by a Normal distribution where the mean and standard deviation is taken from the historical lags from the last two weeks to the reported numbers 2 . Figure 2 depicts four randomly chosen dates where the model is compared to the benchmark. The model and the benchmark are tasked with predicting the total number of individuals who have died at a given date and have been reported within 14 days of that date. As time progresses, more deaths are reported and the dashed grey line approaches the horizontal line. Meanwhile model uncertainty decreases. Figure 3 shows model performance compared to the benchmark for three difference performance metrics. All three graphs are based on predictions of reported deaths within 14 days, and show how performance increases as more data has been reported. Each data point is the average of all dates where predictions can be evaluated. SCRPS is a measure of accuracy that rewards precision, it is a proper scoring rule like the continuous probability rank score or the Brier score (see definition in Appendix) (Bolin and Wallin, 2019 ). The central plot shows the width of the prediction intervals, and the rightmost one the proportion of PIś that cover the true value.

Benchmark and model point estimates are similarly close to the truth. The model produces tighter prediction intervals. For 8-5 days of reporting lag (see Figure 3 ), the intervals are too tight. This is likely because the Public Health Agency queries the Swedish death registry for Covid-19 deaths only once or twice a week. Since we do not know the process, it has not been explicitly modeled. 

The model proposed here can estimate the trends in surveillance data with reporting delays, such as the daily COVID-19 reports in Sweden. To generate accurate estimates of the actual event frequencies based on these reports is highly relevant and can have large implications for interpretations of the trends and evolution of disease outbreaks. In Sweden, delays are considerable and exhibit a weekday and holiday pattern that need to be accounted for to draw conclusions from the data. The method and algorithm proposed overcomes major shortcomings in the daily interpretation and practice analyzing and controlling the novel Corona virus pandemic. It also provides valuable measures of uncertainty around these estimates, showing users how large the range of possible outcomes can be. Whenever case statistics are collected from multiple sources and attributed to its actual event date in the middle of a public health emergency, similar reporting delays to the ones in Sweden will necessarily occur. The method described thus has implications and value beyond Sweden, for any situation where nowcasts of disease event frequencies are of relevance to public health.

Nevertheless, the method also has its limitations. As presented, the model assumes that all deaths are reported in the same manner. Given there exists many regions in Sweden this is unlikely to be the case. For example, it is easy to see that the Swedish region Västra Götland follows a different reporting structure than Stockholm. Building a model for each region separately would most likely give better results and make the assumptions more reasonable. Unfortunately we do not currently have access to the high resolution data required to do so.

Moreover, deaths are reported from two distinct populations that seem to follow different trends. At the time of writing, the daily deaths in elderly care, reported with a longer delay, seem to be decreasing slower than hospital deaths. But statistics offer only aggregate numbers, prohibiting us from modeling two distinct processes. However, we have noted a clear decline in proportions of deaths reported the two first working days. For example the number of deaths occurring at the second of April ≈ 30% of deaths where reported within the first two working days whereas for the eighteens of May only ≈ 10% where reported during the two first working days. We address this by assuming that the deaths reported during the two first working days comes from a different population then the remainder of days.

Another limitation is that the model assumes that the number of new reported deaths for a given day cannot be negative, which is not actually true, due to miscount or misclassification of days. The number of such cases is very small, however, and its removal should not make much difference. The central assumption of the model is that the proportions deaths reported each day is fixed (up to the known covariates). If actual reporting standards change over time, the model will not be able to account for this. But reporting likely becomes faster as the crisis infrastructure improves. One can imagine that after a while the reporting improves, or is changed, if this is not accounted for by a covariate in the model, it will report incorrect numbers. Of course, there might be unknown variables that we have failed to incorporate, but at the least the model is an improvement from the estimates using moving averages. When the covariates to the reporting delay pattern are known, the model can incorporate them and provide more accurate predictions.

In this paper, we provide a method to accurately nowcast daily Covid-19 statistics that are reported with delay. By systematically modelling the delay, policy makers can avoid dangerous illusory downward trends. Our model also gives precise uncertainty intervals, making sure users of these statistics are aware of the fast-paced changes that are possible during this pandemic.

Death date r 11 r 12 · · · · · · r 1T r 22 · · · · · · r 2T r 33 · · · r 3T . . . . . . p ij , i.e.

Typically in removal sampling one would set the probability of reporting uniform, i.e. p i,j := p. However for this data this is clearly not realistic given weekly patterns in reporting -very little reporting during the weekends. Instead we assume that we have k different probabilities. Further, to account for overdispertion, we assume that each probability rather being a fixed scalar is a random variable with a Beta distribution. The Beta distribution has two parameters α and β. This resulting the following distribution for the probabilities

). Here, if j ∈ H then day j is a holidays or weekends, and the parameters above are

else.

These extra parameters are created to account for the under-reporting that occurs during weekend and holidays. Finally we add an extra mixture component that allows for very low reporting.

For the α and β parameters we use an (improper) uniform prior. For the deaths, d, one could imagine several different prior ideally some sort of epidemiological model. However, here we just assume a log-Gaussian Cox processes (Møller et al., 1998) , but instead of Poisson distribution we use a negative binomial to handle possible over dispersion. The latent Gaussian processes has a intrinsic random walk distribution (Rue and Held, 2005) i.e.

This model is created to create a temporal smoothing between the reported deaths. For the hyperparameter σ 2 we impose a inverse Gamma distribution, this prior is suitable here because it guarantees that the process is not constant (σ 2 = 0) which we know is not the case.

Putting the likelihood and priors together we get the following hierarchical Bayesian model

where where and j ≤ i and i = 1, . . . , T .

As the main goal to generate inference of the number of death d is through the posterior distribution of number of deaths d given the observations r. In order to generate samples from this distribution we use a Markov Chain Monte Carlo method (Brooks et al., 2011) . In more detail we use a blocked Gibbs sampler, which generates samples in the following sequence:

• We sample α, β, α H , β H |d, r using the fact that one can integrate out p in the model, and then d|α, β, α H , β H , r, λ follows a Beta-Binomial distribution. Here to we use an adaptive MALA (Atchadé, 2006) to sample from these parameters.

• To sample d|α, β, α H , β H , r, λ, that each death, d i is conditionally independent, and we just use a Metropolis Hastings random walk to sample each one.

• To sample λ|d, σ 2 we again use an adaptive MALA.

• Finally We sample σ 2 |d,and p 0 , π directly since this distribution is explicit, and φ using a MH-RW.

In this section, we present additional comparison of the model to the benchmark. We first describe the benchmark model in detail. The benchmark model simply takes the sum of average historical reporting lags for the preceding 14 days. As before r ij is the number of deaths that happened on day i and were recorded on day j. To predict the number of people that died on a given day, we first calculate lag averages:

wherer i,i+L is the average number of deaths reported with a lag of L days, based on the 14 reports closest preceding day i. If we are looking at data released 2020 − 04 − 28 and call this day 0, the latest death date that we have 10-day (L = 10) reporting lag observation for is r −10,0 . The average for Lag(0, 10) is therefore taken over the 14 days between r −24,−14 and r −10,0 (2020-04-04 and 2020-04-18). For this reason, some of the earlier predictions will not have data from 14 days. The average is then taken over all available reports.

In the comparisons we aim at predicting the total number of deaths that will have been reported within 14 days of the death date. To do so, we sum over the average lag that has yet to be reported. If we are predicting the number of people that have yet to be reported dead for day -3, we already know the true values for r −3,−3 , r −3,−2 , r −3,−1 , and r −3,0 so we only need to predict r −3,1 . . . r −3,10 . The prediction is then Benchmark(i, j) = j l=i r i,l + 14 l=jr i,l .

(2)

As confidence interval we simply use a Normal assumption with standard deviations of the reporting lags, assuming independence, i.e. this is just the square root of the sum of V ar(r). 

How will country-based mitigation measures influence the course of the COVID-19 epidemic?

An Adaptive Version for the Metropolis Adjusted Langevin Algorithm with a Truncated Drift

Scale dependence: Why the average crps often is inappropriate for ranking probabilistic forecasts

Handbook of Markov Chain Monte Carlo

Evaluation of reporting timeliness of public health surveillance systems for infectious diseases

An Attempt to Determine the Absolute Number of Rats on a Given Area

Open models for removal data

Log Gaussian Cox Processes

A Mathematical Theory of Animal Trapping

Coronavirus: Three things all governments and their science advisers must do now

Antalet virusdöda har underskattats

Review Papers: Modeling Capture, Recapture, and Removal Statistics for Estimation of Demographic Parameters for Fish and Wildlife Populations: Past, Present, and Future

Gaussian Markov Random Fields: Theory and Applications

A Appendix

Before presenting the model we describe some notation used through out the appendix. For a m × n matrix r we use the following broadcasting notation r k,j:l = [r k,j , r k,j+1 , . . . , r k,l ]. Further x|y ∼ π(.) implies that the random variable x if we conditioning on y follows distribution π(.). The relevant variables in the model are the following:Variable name Dimension DescriptionLatent prior parameter for p α H 2 × 1 parameter for the probability, p for holiday adjustment. β H 2 × 1 parameter for the probability, p for holiday adjustment. µ T × 1 µ i is the intensity of the expected number of deaths at day i. σ 2 1 × 1 Variation of the random walk prior for the log intensity. φ 1 × 1 overdispersion parameter for negative binomial distribution. p 0 1 × 1 probability of reporting for a low reporting event. pi 1 × 1 probability of a low reporting event.

The most complex part of our model is the likelihood, i.e. the density of the observations given the parameters. Here the data consist the daily report of recorded deaths for the past days. This can conveniently be represented upper triangular matrix, r, where r i,j represents number of new reported deaths for day i reported at day j. This matrix is displayed on the left in Table 1 . We assume that given the true number of deaths at day i, d i , that each reported day j the remaining death d i − j−1 k=1 r i,k each recored with probability