Lizenz: CC BY 4.0
Source code: Github Gist

Indikatorregel:

50 Neuinfektionen pro 100.000 Einwohnern in einem Landkreis innerhalb von sieben Tagen.

Annahme: Operationalisert wird das durch die Anzahl an gemeldeten Neuinfektionen innerhalb einer Woche (Zeitpunkt der zählt: Tag an dem der Fall beim Gesundheitsamt ankommt).

# Verzug zwischen Neuinfektion und Ansteckung
dauer_neuinf_ansteckend <- 4 #Generationsintervall
# Verzug zwischen Neuinfektion und Meldung (Inkubationszeit + Dauer Erkrankungsbeginn bis Meldung am GA)
dauer_neuinf_meldung    <- 5 + 5 


##################################################################################
# Vereinfachte Simulation in einer Bevölkerung mit 100,000 Einwohnern.
#
# @param Rt Zeitlich variierende Reproduktionszahl
# @param n_equilibrium Anzahl Fälle in einer stabilen Situation, um den es zunächst mit R(t)=1 schwankt
# @param n_steps Anzahl Zeit Schritte in der Simulation
#
# @return Situation wenn es mehr als 50 gemeldete Fälle gibt, d.h. tatsächliche
##        Anzahl Neuinfektionen
##################################################################################

simulate <- function(Rt = 1.2,   n_equilibrium =5, n_steps = 50) {

  # Tägliche Anzahl Neuinfektionen
  neuinfektionen <- c(rep(n_equilibrium, 7), rep(NA, n_steps-7))
  # Wochensumme der Neuinfektionen
  neuinfektionen_summe <- c(rep(7*n_equilibrium, 7), rep(NA, n_steps-7))
  
  # Gemeldete Neuinfektionen
  neuinfektionen_gemeldet <- c(rep(n_equilibrium,7+dauer_neuinf_meldung), rep(NA, n_steps-7-dauer_neuinf_meldung))
  # Summe für die Notbremse
  neuinfektionen_gemeldet_summe <- rep(NA, n_steps)
  # Alarm
  notbremse <- rep(NA, n_steps)
  
  # Monitoring
  for (t in 8:(n_steps-dauer_neuinf_meldung)) {
    # Infektionsprozess
    neuinfektionen[t] <- Rt*neuinfektionen[t-dauer_neuinf_ansteckend]
    neuinfektionen_summe[t] <- sum(neuinfektionen[(t-6):t])
    
    # Meldeprozess (Wann werden die Neuinfektionen gemeldet?
    neuinfektionen_gemeldet[t+dauer_neuinf_meldung] <- neuinfektionen[t]
    
    # Indikator - summe über die letzten 7 Tage
    neuinfektionen_gemeldet_summe[t] <- sum( neuinfektionen_gemeldet[(t-6):t] )
    
    #Notbremse ziehen?
    notbremse[t] <- neuinfektionen_gemeldet_summe[t] > 50
  }
  
  # Ergbnis als data.frame
  ts <- data.frame(t=seq_len(n_steps), neuinfektionen=neuinfektionen, neuinfektionen_summe=neuinfektionen_summe, neuinfektionen_gemeldet=neuinfektionen_gemeldet, neuinfektionen_gemeldet_summe=neuinfektionen_gemeldet_summe, notbremse=notbremse) %>% slice(1:(n_steps - dauer_neuinf_meldung))
  ts
  
  # Situation, falls "notbremse" gerissen wird, sowie gesamter Verlauf.
  res <- ts %>% filter(notbremse) %>% slice(1) %>% mutate(Rt=Rt) %>% select(Rt, everything()) %>%  
    mutate(zuviel_faktor =  neuinfektionen_summe/50 )
  attr(res, "thesim") <- ts
  return(res)
}

# Durchspielen von drei R(t) Szenarien. Angezeigt wird jeweils der erste Zeitpunkt
# wo die Notbremse auslöst.
Rt_vec <- c(1.1, 1.2, 1.3)

sims <- map(Rt_vec, simulate) 
sims_zuviel <- sims %>% bind_rows() %>% mutate(zuviel_faktor =  neuinfektionen_summe/50 )
sims_zuviel
##    Rt  t neuinfektionen neuinfektionen_summe neuinfektionen_gemeldet
## 1 1.1 34       9.743586             64.66198                 8.05255
## 2 1.2 27      12.441600             80.87040                 8.64000
## 3 1.3 24      18.564650             97.65665                 8.45000
##   neuinfektionen_gemeldet_summe notbremse zuviel_faktor
## 1                      50.64455      TRUE      1.293240
## 2                      52.08000      TRUE      1.617408
## 3                      51.35000      TRUE      1.953133

Grobe Schätzung des zu-viel Faktors

Der Faktor lässt sich grob abschätzen als:

\[ R(t)^{\left\{\frac{\text{Mittlere Dauer von Exposition bis Meldung beim GA}}{\text{Mittlere Generationszeit}}\right\}} \] Herleitung: Aus Wallinga und Lipsitch (2007) ist bekannt, dass der Zusammenhang zwischen exponentiellen Anstieg und der Reproduktionszahl \(R\) bei einer konstanten Generationszeit von \(G\) Tagen gleich \(R = \exp(r \cdot G)\) ist, wobei \(r\) die pro Capita Rate des Anstiegs ist. Somit ist der multiplikative Ansiteg pro Tag bei einer Reproduktionszahl von \(R\) gleich \[ \exp(r) = \exp\left( \log(R) \cdot \frac{1}{G}\right) = R^{\frac{1}{G}}. \] Die Anzahl an gemeldeten Fällen zum Tag \(t\) entspricht der Anzahl an Neuinfektionen zum Zeitpunkt \[t-\text{Mittlere Dauer von Exposition bis Meldung beim GA}.\] Diese mittlere Dauer setzt sich aus mittlere Dauer von Ansteckung bis Symptombeginn sowie mittlere Dauer von Symptombeginn bis Meldung des Falles beim GA zusammen. Die beiden Dauern zusammen betragen im Mittel ungefährt 5+5=10 Tage. Legt man also den Schwellenwert auf der Anzahl an gemeldeten Fällen zum Zeitpunkt \(t\) an, ist die tatsächliche Anzahl an Neuinfektionen zum gleichem Zeitpunkt unter obigen Annahmen bereits um den Faktor \(\exp(r)^{10} = R^{\frac{10}{G}}\) höher als die Meldezahlen.

In R-Code:

approx <- data.frame(Rt=Rt_vec, zuviel_faktor=Rt_vec^(dauer_neuinf_meldung / dauer_neuinf_ansteckend))

Vergleich zwischen Approximation und Simulation

inner_join(sims_zuviel %>% select(Rt, zuviel_faktor), approx, by="Rt", suffix=c(".sim",".approx")) %>% knitr::kable(digits=2)
Rt zuviel_faktor.sim zuviel_faktor.approx
1.1 1.29 1.27
1.2 1.62 1.58
1.3 1.95 1.93

Plots

# Funktion um eine Simulation für einen Rt-Wert zu erstellen und zu plotten.
make_and_plot_sim <- function(Rt) {
  #Simulate
  sims <- simulate(Rt)
  zuviel_faktor <- sims %>% pull(neuinfektionen_summe) / 50
  attr(sims, "thesim")  %>% select(t,neuinfektionen, neuinfektionen_summe) 
  
  t_notbremse <- sims %>% filter(notbremse) %>% pull(t)
  
  # Long format
  df <- attr(sims, "thesim") %>% pivot_longer(cols=-t, names_to = "Zeitreihe", values_to = "Anzahl")
  
  ## Plot für den R(t) Wert
  df_relevant <- df %>% filter(t <= t_notbremse+1)
  p1 <- ggplot(df_relevant %>% filter(Zeitreihe %in% c("neuinfektionen", "neuinfektionen_gemeldet")), aes(x=t, y=Anzahl, color=Zeitreihe)) + geom_step(direction="vh") + 
    theme(legend.position = 'bottom') +
    scale_color_brewer(type="qual") +
    ggtitle(sprintf("R(t) = %.2f", Rt))  
   
  p2 <- ggplot(df_relevant %>% filter(Zeitreihe %in% c("neuinfektionen_summe", "neuinfektionen_gemeldet_summe")), aes(x=t, y=Anzahl, color=Zeitreihe)) + geom_step(direction="vh") + 
    theme(legend.position = 'bottom') +
    geom_hline(yintercept= 50, col="salmon2", lty=2) +
    geom_vline(xintercept=t_notbremse, col="salmon2", lty=2) +
    scale_color_brewer(type="qual") + 
    ggtitle(sprintf("Zuviel bei Notbremse: %.2f Mal", zuviel_faktor)) 
   
  gridExtra::grid.arrange(p1, p2, nrow=2, ncol=1)
  invisible()
}

make_and_plot_sim(Rt=1.2)
## Warning: Removed 7 row(s) containing missing values (geom_path).

Interaktive Visualisierung

Geht nur in R-Studio

library(manipulate)
manipulate(
  make_and_plot_sim(Rt),
  Rt=slider(1,2, initial=1.2, step=0.1))

Literatur

  • Wallinga, J. und Lipsitch, M. (2007), How generation intervals shape the relationship between growth rates and reproductive numbers, Proc. R. Soc. B., 274:599-604. https://doi.org/10.1098/rspb.2006.3754.