Nowcasting fatal COVID-19 infections on a regional level in Germany

Before we discuss our modelling approach in detail, we want to describe our major findings. First, Table 1 shows that age and gender both play a major role when estimating the daily death toll. As is generally known, elderly people exhibit a much higher death rate which is for the age group 80+ around 100 times higher than for people in the age group 35-59. A remarkable difference is also observed between genders, where the expected death rate of females is around 40% ($\approx 1-\exp(-0.503$)) lower than the death rate for males. Furthermore, we see that significantly less deaths are attributed to infections registered on Sundays compared to weekdays, due to the existing reporting delay during weekends.

Our model includes a global smooth time trend representing changes in the death rate since March 26th. This is visualized in Figure 1. The plotted death rate is scaled to give the expected number of deaths per 100.000 people in an average district for the reference group, i.e. males in the age group 35 - 59. Overall, we see a peak in the death rate on April 3rd and a downwards slope till end of April. However, our nowcast reveals that the rate remains constant since beginning of May. Note that this recent development cannot be seen by simply displaying the raw death counts of these days. The nowcasting step inevitably carries statistical uncertainty, which is taken into account in Figure 1 by including best and worst case scenarios. The latter are based on bootstrapped confidence intervals, where details are provided in Section 6.3 later in the paper.

Our aim is to investigate spatial variation and regional dynamics. To do so, we combine a global geographic trend for Germany with unstructured region-specific effects, where the latter uncover local behaviour. In Figure 2 we combine these different components and map the fitted nowcasted death counts related to Covid-19 for the different districts of Germany, cumulating over the last seven days before the day of analysis (here May 14, 2020). While in most districts of Germany the death rate is relatively low, some hotspots can be identified. Among those, Traunstein and Rosenheim (in the south-east part of Bavaria) are the most evident, but Greiz and Sonneberg (east and south part of Thuringia) stand out as well, to mention a few. A deeper investigation of the spatial structure is provided in Section 6, where we show the global geographic trend and provide maps that allow to detect new hotspot areas, after correcting for the overall spatial distribution of the infection.

On the day of analysis, we do not observe the total counts of deaths for recently registered infections, since not all patients with an ongoing fatal infections have died yet. We therefore nowcast those numbers, i.e. we predict the prospective deaths which can be attributed to all registration dates up to today. This is done on a national level, and the resulting nowcast of fatal infections for Germany is shown in Figure 3. For example, on May 14, 2020 there are 25 deaths reported where the infection was registered on May 5th (red line on May 5th). We expect this number to increase to about 50 when all deaths due to Covid-19 for this registration date will have been reported (blue line on May 5th). Naturally, the closer a date is to the present, the larger the uncertainty in the nowcast. This is shown by the shaded bands. Details on how the statistical uncertainty has been quantified are provided in Section 5 below. The fit of this model has been incorporated into the district model discussed before, but the nowcast results are interesting in their own right. The curve confirms that the number of fatal infections is decreasing since the beginning of April. Note that the curve also mirrors the ”weekend effect” in registration, as less infections are reported on Sundays. Further analyses and a detailed description of the model are given in the following sections.

The above model cannot be fitted directly to the available data, since we need to take the course of the disease into account. For a given registration date $t$, the number of deaths of patients registered as positive on that day, $Y_{t,r,g}$, may not yet be known, since not all patients with a fatal outcome of the disease have died yet. This requires the implementation of nowcasting. We do this on a national level, and cumulate the numbers over district/region $r$ and gender and age groups $g$. This allows to drop the corresponding subscripts in the following and we simply notate the cumulated number of deaths with registered infections at day $t$ with $Y_{t}$. Let $N_{t,d}$ denote the number of deaths reported on day $t+d$ for infections registered on day $t$. Assuming that the true date of death is at $t+d$, or at least close to it, we ignore any time delays between time of death and its notification to the health authorities. We call $d$ the duration between the registration date as a Covid-19 patient and the reported day of death, where $d=1,\ldots,d_{max}$. Here, $d_{max}$ is a fixed reasonable maximum duration, which we set to 30 days (see e.g. Wilson et al. 2020). The minimum delay is one day. In nowcasting we are interested in the cumulated number of deaths for infections registered on day $t$, which we define as

The total number of deaths with a registered infection at $t$ is apparently unknown at time point $t$ and becomes available only after $d_{max}$ days. In other words, only after $d_{max}$ days we know exactly how many deaths occurred due to an infection which was registered on day $t$. We define the partial cumulated sum of deaths as

so that by definition $C_{t,d_{max}}=Y_{t}$.

On day $t=T$, when the nowcasting is performed, we are faced with the following data constellation, where NA stands for not (yet) available:

We may consider the time span between registered infection and (reported) death as a discrete duration time taking values $d=1,\ldots,d_{max}$. Let $D$ be the random duration time, which by construction is a multinomial random variable. In principle, for each death we can consider the pairs $(D_{i},t_{i})$ as i.i.d. and we aim to find a suitable regression model for $D_{i}$ given $t_{i}$, including potential additional covariates $x_{t,d}$. We make use of the sequential multinomial model (see Agresti 2010) and define

Let $F_{t}(d)$ denote the corresponding cumulated distribution function of $D$ which relates to probabilities $\pi()$ through

for $d=1,\dots,d_{\max}-1$ and $F_{t}(d_{\max})=1$.

The available data on cumulated death counts allow us to estimate the conditional probabilities $\pi(d;)$ for $d=2,\dots,d_{\max}$. In fact, the sequential multinomial model allows to look at binary data such that

with

where

• •

  $s_{1}(t)$ is an overall smooth time trend over calendar days,

• •

  $s_{2}(d)$ is a smooth duration effects, capturing the course of the disease,

• •

  $x_{t,d}$ are covariates which may be time and duration specific.

Assuming that $D$, the duration between a registered fatal infection and its reported death, is independent of the number of fatal Covid-19 infections, we obtain the relationship

Note further that if we model $Y_{t}$ with a quasi-Poisson model as presented in the previous chapter, we have no available observation $Y_{t}$ for time points $t>T-d_{max}$. Instead, we have observed $C_{t,T-t}$, which relates to the mean of $Y_{t}$ through (7). Including therefore $\log F_{t}(T-t)$ as additional offset in model (2), allows to fit the model as before, but with nowcasted deaths included. That means, instead of $\lambda_{t,r,g}$ as in (2), the expected death rates are now parametrized by $\lambda_{t,r,g}^{\star}=\lambda_{t,r,g}\exp(\log F_{t}(T-t))$, where the latter multiplicative term is included as additional offset in the model.

We fit the nowcasting model (5) with parametrization (6). We include a weekday effect for the registration date of the infection with reference category ”Monday”. The estimates of the fixed linear effects are shown in Table 2. The fitted smooth effects are shown in Figure 4, where the top panel shows the effect over calendar time, which is very weak and confirms that the course of the disease hardly varies over time. This shows that the German health care system remained stable over the considered period, and hence survival did not depend on the date on which the infection was notified.

The bottom panel of Figure 4 shows the course of the disease as a smooth effect over the time between registration of the infection and death. We see that the probabilities $\pi(d;\cdot)$ decrease in $d$, where this effect is the strongest in the first days after registration. Thus, most of the Covid-19 patients with fatal infections are expected to die not long after their registration date. The effect of $d$ becomes easier to interpret by visualizing the resulting distribution function $F_{t}(d)$. This is shown in Figure 5 for two dates $t$, i.e.. April 14th and May 13th. The plot also shows how the course of the disease hardly varies over calendar time: In fact, the small differences between the two distribution functions is dominated by the weekday effect, since the red curve is related to a Tuesday while the blue one is from a Wednesday.

In Figure 3 above we have shown the nowcasting results along with uncertainty intervals shaded in grey. These were constructed using a bootstrap approach as follows. Given the fitted model, we simulate $n=10\text{ }000$ times from the asymptotic joint normal distribution of the estimated model parameters which results through (4). This leads to a set of bootstrapped distribution functions $\mathcal{F}=\{\widehat{F}_{t}^{(i)}(T-t),i=1,\dots,n;t=T-d_{\max}+1,\dots,T-1\}$. This set is used to compute the simulated nowcasts $\widehat{Y}_{t}^{(i)}=C_{T-t}/\widehat{F}_{t}^{(i)}(T-t)$ applying (7), where $C_{T-t}$ is the observed partial cumulated sum of deaths at time point $T-t$. The point-wise lower and upper bounds of the 95% prediction intervals for the nowcast for $Y_{t}$ are then given by the 2.5 and the 97.5 quantiles of the set $\{\widehat{Y}_{t}^{(i)},i=1,\dots,n\}$, respectively.

In Section 3 we presented the fitted death rate, which is the convolution of a smooth spatial effect as well as region specific effects. It is of general interest to disentangle these two spatial components. This is provided by the model. We visualize the fitted global geographic trend $m_{2}(\cdot)$ for Germany in Figure 6. The plot confirms that up to May 2020 the northern parts of the country are less affected by the disease in comparison to the southern states. The two plots in Figure 7 map the region specific effects, i.e. the predicted long term level of a district $u_{r0}$ (left hand side) and the predicted short term dynamics $u_{r1}$ (right hand side). Both plots uncover quite some region-specific variability. In particular, the short term dynamics captured in the right hand side plot ($u_{r1}$) pinpoint districts with unexpectedly high nowcasted death rates in the last two weeks, after correcting for the global geographic trend and the long term effect of the district. Some of the noticeable districts have already been highlighted in Section 3 above, but we can detect further districts, which are less pronounced in Figure 1. For instance, Steinfurt (in the north-west of North Rhine-Westphalia), Olpe (southern North Rhine-Westphalia) or Gotha (center of Thüringen) presently show a high rate of fatal infections.

A large number of the registered deaths related to Covid-19 stem from people in the age group 80+. Locally increased numbers are often caused by an outbreak in a retirement home. Such outbreaks apparently have a different effect on the spread of the disease, and the risk of an epidemic infection caused by outbreaks in this age group is limited. Thus, the death rate of people in the age group 80+ could vary differently across districts when compared to regional peaks in the death rate of the rest of the population. In order to respect this, we decompose the district-specific effects $\bm{u}_{r}$ in (2) into $\bm{u}_{r}^{80-}=(u_{r0}^{80-},u_{r1}^{80-})^{\top}$ for the age group 80- and $\bm{u}_{r}^{80+}=(u_{r0}^{80+},u_{r1}^{80+})^{\top}$ for the age group 80+, where the age group 80- consists of the aggregated age groups 15-34, 35-59 and 60-79. We put the same prior assumption on the random effects as we did in (3), but now the variance matrix that needs to be estimated from the data has dimension 4 by 4.

The fitted age group-specific random effects are shown in Figure 8, where the $\bm{u}_{r}^{80-}$ are shown in the top panel and the $\bm{u}_{r}^{80+}$ in the bottom panel. Most evidently, the variation of the random effects is much higher in the age group 80+ when compared to the younger age groups, as more districts occur which are coloured dark blue or dark red, respectively. When comparing the district-specific short term dynamics of the last 14 days ($u_{r1}$) in Figure 8 to those in Figure 7, we recognize that in most of the districts which recently experienced very high death intensities (with respect to the whole period of analysis), these stem from the age group 80+. As mentioned before, this can often be explained by outbreaks in retirement homes.

When fitting the mortality model (1) we included the fitted nowcast model as offset parameter. This apparently neglects the estimation variability in the nowcasting model, which we explored via bootstrap as explained in Section 5.3 and visualized in Figure 3. In order to also incorporate this uncertainty in the fit of the mortality model, we refitted the model using (a) the upper end and (b) the lower end of the prediction intervals shown in Figure 3. It appears that there is little (and hardly any visible) effect on the spatial components, which is therefore not shown here. But the time trend shown in Figure 1 does change, which is visualized by including the two fitted functions corresponding to the 2.5% and 97.5% quantile of the offset function. We can see that the estimated uncertainty of the nowcast model mostly affects the last ten days, with a strong potential increase in the death rate mirroring a possible worst case scenario.

In Figure 9 we show a normal QQ-plot of the Pearson residuals in the nowcasting model. Apart from some observations in the lower tail, the Pearson residuals are distributed very closely to a standard normal distribution when considering the estimate $\widehat{\phi}=1.766$ of the dispersion parameter in the quasi-poisson model (7). Overall, the model seems to fit to the available data quite well.

The paper presents a model to monitor the dynamic behaviour of Covid-19 infections based on death counts. It is important to highlight that the proposed model makes no use of new infection numbers, but only of observed deaths related to Covid-19. This in turn means that the results are less dependent on testing strategies. The nowcasting approach enables us to estimate the number of deaths following a registered infection today, even if the fatal outcome has not occurred yet. Moreover, the district level modelling uncovers hotspots, which are salient exclusively through increased death rates. A differential analysis of the number of current fatal infections on a regional level allows to draw conclusions on the current dynamics of the disease assuming a constant case fatality rate, i.e. a stable proportion of death compared to the true number of infections when adjusting for age and gender.

A natural next step would now be to consider the nowcasted deaths in relation to the number of newly registered infections, which is, in contrast, highly dependent on both testing strategy and capacity. We consider this as future research, and the proposed model allows us to explore data in this direction. This might ultimately help us in shedding light on the relationship between registered and undetected infections as well as on the effectiveness of different testing strategies.

There are several limitations to this study which we want to address as well. First and utmost, even though death counts are, with respect to cases counts, less dependent on testing strategies, they are not completely independent from them. This applies in particular to the handling of post-mortem tests. We therefore do not claim that our analysis of death counts is completely unaffected by testing strategies. Secondly, a fundamental assumption in the model is the independence between the course of the disease and the number of infections. Overall, if the local health systems have sufficient capacity and triage can be avoided, this assumption seems plausible, but it is difficult or even impossible to prove the assumption formally. Finally, the nowcasting itself is not carried out on a regional level, though the model focuses on regional aspects of the pandemic. While it would be desirable to fit the nowcast model regionally, the limited amount of data simply prevents us from extending the model in this direction.