Lizenz: CC BY 4.0
Source code: Github Gist
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
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))
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 |
# 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).
Geht nur in R-Studio
library(manipulate)
manipulate(
make_and_plot_sim(Rt),
Rt=slider(1,2, initial=1.2, step=0.1))