mhoehle/notbremse.Rmd

## notbremse.Rmd
---
title: "50 von 100,000 Notbremse innherhalb von 7 Tagen"
author: "Michael Höhle ([https://www.math.su.se/~hoehle](https://www.math.su.se/~hoehle))"
date: "`r Sys.Date()`"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
# Lade Pakete
suppressPackageStartupMessages(library(tidyverse))
```

Lizenz: [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/)<br>
Source code: [Github Gist](https://gist.github.com/hoehleatsu/9e93402949c291b030a5e56c626ccf23)

# 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).

```{r}


# 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
```

## 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:

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


## Vergleich zwischen Approximation und Simulation
```{r}
inner_join(sims_zuviel %>% select(Rt, zuviel_faktor), approx, by="Rt", suffix=c(".sim",".approx")) %>% knitr::kable(digits=2)
```


## Plots
```{r}
# 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)
```

# Interaktive Visualisierung

Geht nur in R-Studio

```{r, eval=FALSE}
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.
	---
	title: "50 von 100,000 Notbremse innherhalb von 7 Tagen"
	author: "Michael Höhle ([https://www.math.su.se/~hoehle](https://www.math.su.se/~hoehle))"
	date: "`r Sys.Date()`"
	output: html_document
	---

	```{r setup, include=FALSE}
	knitr::opts_chunk$set(echo = TRUE)
	# Lade Pakete
	suppressPackageStartupMessages(library(tidyverse))
	```

	Lizenz: [CC BY 4.0](https://creativecommons.org/licenses/by/4.0/)<br>
	Source code: [Github Gist](https://gist.github.com/hoehleatsu/9e93402949c291b030a5e56c626ccf23)

	# 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).

	```{r}


	# 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
	```

	## 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:

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


	## Vergleich zwischen Approximation und Simulation
	```{r}
	inner_join(sims_zuviel %>% select(Rt, zuviel_faktor), approx, by="Rt", suffix=c(".sim",".approx")) %>% knitr::kable(digits=2)
	```


	## Plots
	```{r}
	# 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)
	```

	# Interaktive Visualisierung

	Geht nur in R-Studio

	```{r, eval=FALSE}
	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.
No results found