Auf geht es in eine neue Runde COVID-19 Datenanalyse mit R. Thema heute: Wie hängen Zahl der Tests und Zahl der berichteten Infektionen zusammen?
Die Behauptung, steigende Infektionszahlen seien dadurch begründet, dass mehr getestet wird, dürfte jedem schon begegnet sein. Erfreulicherweise finden sich beim RKI nicht nur die Infektionszahlen sondern auch Daten zu Tests:
- Infektionen (Rohdatei, csv): https://npgeo-corona-npgeo-de.hub.arcgis.com/datasets/dd4580c810204019a7b8eb3e0b329dd6_0
- Testzahlen: https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Testzahl.html
Dass die Zahl der Tests und die Zahl der Infektionen zusammenhängt ist offensichtlich und trivial: Ohne Tests, kann es keine berichteten Infektionen geben. Die berichtete Infektionszahl ist abhängig von der Testzahl, es geht gar nicht anders.
Ob Veränderungen in er Teststrategie oder das tatsächliche Infektionsgeschehen ursächlich für die veränderten berichteten Infektionszahlen sind ist eine komplett andere Frage. Sie lässt sich anhand der Zeitreihendaten analysieren.
Ich habe die Daten normalisiert (auf eine Skala gebracht) und geplottet. Infektionszahlen, werden täglich von Gesundheitsämtern an das RKI übermittelt, Testzahlen wöchentlich von Laboren. Deswegen sind Testzahl und Positivrate „eckig“.
Diese Hypothese lässt sich mit den Daten prüfen. Prüfen heißt im wissenschaftlichen Sinn, die Hypothese versuchen mit Daten zu widerlegen (Falsifikationsprinzip). Dabei helfen die Daten von Juli 2020 bis Ende September 2020:
Veränderung der Testzahl ohne Effekt auf Infektionszahl
Es ist von immenser Bedeutung zu verstehen, dass mit dieser Prüfung nur widerlegt wurde, dass die Testzahl der einzige Faktor ist, der die Infektionszahl beeinflusst. Er zeigt, dass es noch andere Faktoren geben muss, die die berichtete Infektionszahl beeinflussen, allen voran (aber nicht ausschließlich!) das Infektionsgeschehen. (Für Statistik-Interessierte: Stichwort aufgeklärte Varianz)
Es gibt klare Effekte von Testzahl und berichteter Infektionszahl. Diese zeigen sich in den Daten von Oktober 2020 bis Februar 2021.
Effekte von Testzahlen auf Infektionszahlen
- Im Oktober steigt die Zahl der Tests stark, gleichzeitig steigt die Zahl der Infektionen sehr schnell -> Die Steigerung der Testzahl führt zu einer (leichten) Überschätzung der tatsächlichen Infektionsdynamik
- Anfang November sinkt die Zahl der Tests, die Infektionszahlen stagnieren -> Die Reduzierung der Testzahl führt zu einer Unterschätzung der tatsächlichen Infektionsdynamik
- Ab Dezember steigen die Testzahlen wieder, die Infektionszahlen steigen auch-> Die Verschärfung des Lockdowns ist die logische Konsequenz, weil im gesamten November die tatsächliche Infektionsdynamik unterschätzt und damit die Effektivität des „Lockdown Light“ überschätzt wurde
Interessant ist, dass die Schwankungen bei der Positivrate deutlich geringer sind. Sie steigt im gesamten Zeitraum einfach an.
Der grünlich markierte Bereich um Weihnachten zeigt einen ähnlichen Effekt. Aufgrund der Feiertage wurde weniger getestet (auch Labore machen Urlaub). Die berichtete Infektionszahl ging sprunghaft zurück. Gleichzeitig erreicht die Positivrate ihr absolutes Maximum.
Gesamthaft betrachtet legen diese Daten nahe, dass die Infektionszahl (und damit der Inzidenzwert) anfällig für Schwankungen der Testzahlen sind. Inzidenz und Infektionszahl (alleine) sind damit nur begrenzt geeignet, um die Infektionsdynamik zu überwachen und steuernde Maßnahmen zu ergreifen.
Die Positivrate scheint anhand dieser Daten ein besserer Indikator für den Infektionsverlauf zu sein. Sie ist von den Schwankungen der Testzahlen im Beispiel relativ unabhängig. Aber auch die Positivrate ist beeinflussbar, insbesondere sobald man nicht auf den Bundesschnitt, sondern auch auf das lokale Geschehen schaut. Daher ist es am sinnvollsten Positivrate und Infektionszahl immer gemeinsam zu betrachten und zu kommunizieren, um einen möglichst klaren Blick auf das Infektionsgeschehen zu erhalten.
Kennzahlen zur weiteren Steuerung
Zusätzlich zur Verwendung von Positivrate und Infektionszahl, ist es sinnvoll noch weitere Variablen zur Bewertung der Lage heranzuziehen. Um ein möglichst umfangreiches Bild, mit möglichst wenig Aufwand zu erzeugen, eignet sich aus meiner Sicht eine Kombination von drei Zahlen die auch verschiedenen Ursprungs sind:
- Positivrate der Tests geliefert von Laboren (unter der Annahme eine halbwegs stabilen Teststrategie)
- Die Zahl der Patienten auf Intensivstationen geliefert von den Krankenhäusern
- Die Zahl der Tode geliefert von den Gesundheitsämtern
Das entsprechende Bild sieht so aus:
Spannend ist der Verlauf der Daten ab Februar 2021. Die Linien laufen auseinander: Die Todeszahl bleibt gering, sinkt sogar, während Positivrate und COVID-19 Patientenzahl steigt. Meine wackelige, auf Hoffnung bauende Hypothese ist, dass dieser Effekt durch Impfungen ermöglicht wird: „Impfungen schützen vor schweren COVID-19 Verläufen und führen so zu weniger Toten, trotz weiterhin hoher Infektionsdynamik.“
Dieser Hypothese widerspricht jedoch der Anstieg der COVID-19 Patienten in Intensivbehandlung. Bisher habe ich keine stichfeste Erklärung, warum die Zahl der Intensivpatienten steigt, aber die Zahl der Toten vergleichsweise niedrig ist.
Die wahrscheinlichsten Faktoren, die eine Rolle spielen könnten, sind die Verbreitung der Mutante B.1.1.7 („B.1.1.7 ist gefährlicher“) und eine möglicherweise veränderte Altersstruktur von Intensivstation-Patienten („Jüngere Patienten liegen länger auf Intensivstationen und sterben nicht“). Da ich aber weder Daten zum Anteil der Mutanten auf Intensivstationen noch zur Altersstruktur der Patienten habe, kann ich diese Hypothesen nicht prüfen.
Alternativ wäre denkbar, dass die meisten Menschen nicht auf Intensivstationen sterben, sondern zu Hause oder in Altenheimen. Dann könnte die Zahl der Toten sinken („zu Hause und im Altenheim wird dank Impfung nicht mehr an COVID-19 gestorben“), während die Zahl der Patienten auf Intensivstationen steigt. Auch hier fehlen mir verlässliche Daten um diese Vermutung zu prüfen.
Die Aufgabe von Politikern und Managern ist, anhand all dieser Informationen zu steuern. Ich glaube, dass die Politik das auch tut. Zumindest kann ich anhand dieser Daten politische Entscheidungen nachvollziehen (s.o. Verschärfung „Lockdown Light“). Leider ist die Kommunikation und Berichterstattung zu diesen Entscheidungen mit viel unnützer Information und wenig Daten bespickt. Es bleibt nötig sich selbst ein Bild zu machen – die Datengrundlage dafür ist frei verfügbar.
R Skript
„Glaube keiner Statistik, die du nicht selbst gefälscht hast“. Deswegen hier wieder die Skripte zur Bearbeitung der Daten und Erstellung der Grafiken. Die Excel mit den Testdaten musste ich lokal speichern und den ersten Reiter rauswerfen. D.h. wer das Skript ohne diesen Schritt ausführt wird Fehlermeldungen bekommen. Obwohl nicht betrachtet, sind hier auch die DIVI Daten wieder drin.
- Download R: https://ftp.fau.de/cran/
- Download R Studio: https://www.rstudio.com/products/rstudio/download/#download
#########Libraries ############
library(dplyr)
library(pracma)
library(BBmisc)
library(plyr)
library(readxl)
library(tidyr)
library(ggplot2)
library(scales)
########### Daten aus DIVI laden und nur die Daten für Deutschland auswählen ##########
#Einlesen
bundesland.zeitreihe <- read.csv("https://diviexchange.blob.core.windows.net/%24web/bundesland-zeitreihe.csv")
#Filtern
divi <- filter(bundesland.zeitreihe, Bundesland == "DEUTSCHLAND")
divi_s <- divi[,-(2:3)]
max(divi$Anzahl_Meldebereiche_Erwachsene)
plot(divi$Anzahl_Meldebereiche_Erwachsene, type = "l")
#Aufräumen
names(divi_s) <- c("Datum", "covCount", "bedsUsed", "bedsFree", "reserve", "ivFree", "ivFreeCov", "sitGood", "sitMedium", "sitBad", "sitNA")
divi_s$Datum <- as.Date(divi_s$Datum)
#Moving Averages berechnen
divi_s$bedsUsed_7 <- movavg(divi_s$bedsUsed,7,"s")
############### Daten vom RKI laden, aggregieren und moving averages der letzten 7 Tage berechnen ##############
#Einlesen
rki <- read.csv("https://opendata.arcgis.com/datasets/dd4580c810204019a7b8eb3e0b329dd6_0.csv")
rki_s <- data.frame(rki$Meldedatum, rki$AnzahlFall, rki$AnzahlTodesfall, rki$AnzahlGenesen)
#Aggregieren
names(rki_s) <- c("Datum", "Infekt", "Tod", "Gesund")
rki_a <- aggregate(rki_s[,2:4], list(rki_s$Datum), FUN = sum)
names(rki_a) <- c("Datum","Infekt", "Tod", "Gesund")
#Moving Averages berechnen
Infekt_7 <- movavg(rki_a$Infekt,7,"s")
Tod_7 <- movavg(rki_a$Tod,7,"s")
Gesund_7 <- movavg(rki_a$Gesund,7,"s")
#Aufräumen
rki_c <- data.frame(rki_a, Infekt_7, Tod_7, Gesund_7)
rki_c$Datum <- as.Date(rki_c$Datum)
####### RKI Fallzahlen und DIVI Daten Mergen und Normalisieren #########
#RKi Fallzahlen und DIVI zahlen mergen und weitere Variablen bereechnen
d <- merge(rki_c, divi_s, by=c("Datum"), all.x = TRUE)
d$beds <- d$bedsFree + d$bedsUsed
d$bedsT <- d$bedsFree + d$bedsUsed + d$reserve
d$ivFree
#Normalisierung von Variablen
z <- as.data.frame(apply(d[,c("Infekt_7", "Tod_7", "covCount", "ivFree", "bedsUsed_7", "beds", "bedsT", "sitGood", "ivFreeCov")], 2, FUN = normalize, method = "range", range = c(0,1)))
z$Datum <- d$Datum
z$ivFree[z$ivFree == 0] <- NA
z$ivFreeCov[z$ivFreeCov == 0] <- NA
########## RKI Testdaten aus lokaler Excel einlesen ##########
#Einlesen und unnötige Zeilen rauswerfen
tests <- read_excel("Documents/Corona Date/Testzahlen-gesamt.xlsx")
tests <- tests[-c((nrow(tests)-(nrow(tests)-1)),nrow(tests)),]
names(tests) <- c("KW", "Zahl", "ZahlPos", "PosAnteil", "Labore")
#Datensatz bereinigen, erste und letzte Zeile entfernen, Daten normalisieren
tz <- as.data.frame(apply(tests[,c("Zahl", "ZahlPos", "PosAnteil", "Labore")], 2, FUN = normalize, method = "range", range = c(0,1)))
tz$KW <- tests$KW
######### RKI Testdaten von Wochenbericht auf Tage strecken ###################
#Vorbereitung des Datensatzes
ndays = nrow(z) # Zahl der Tage, die schon im RKI Bericht sind
td <- matrix(0, nrow = ndays, ncol = 4)
td <- as.data.frame(td)
names(td) <- c("Datum", "Zahl", "ZahlPos", "PosAnteil")
td$Datum <- z$Datum
# Wochenzahlen der Tests für jeden Tag reinschreiben
for(w in 0:(nrow(tests)-1)) {
for(p in 0:6){
td$Zahl[31 + (w*7) + p] <- tz$Zahl[w+1]
td$PosAnteil[31 + (w*7) + p] <- tz$PosAnteil[w+1]
td$ZahlPos[31 + (w*7) + p] <- tz$ZahlPos[w+1]
}
}
#Die Tage, für die es am Ende noch keine Daten gibt großzügig rauswerfen
s <- ndays - 5
for(x in s:ndays){
td$Zahl[x] <- NA
td$PosAnteil[x] <- NA
td$ZahlPos[x] <- NA
}
############ Mergen der DIVI, RKI, und Testdaten #############
total <- merge(z, td, by=c("Datum"), all.x = TRUE)
total$PosAnteil_7 <- movavg(total$PosAnteil, 7, "s")
########### Plots für Artikel zu Tests ##################
plot(total$Datum,total$Infekt_7, type="l",col="purple", main = "SARS-CoV 2 Infektionen, Tests, Positivrate", xlab = "Datum", ylab = "Normalisierte Werte - Wertebereich 0 bis 1")
lines(total$Datum,total$PosAnteil, col="magenta")
lines(total$Datum,total$Zahl, col="cyan")
legend("topleft",
c("Infektionen", "Tests", "Positivrate"),
fill=c("purple", "cyan", "magenta"))
#Zeitbereich komplett
t_plot <- ggplot(data= subset(total, total$Datum > "2020-03-07" & total$Datum < "2021-04-18"), aes(x=Datum)) + geom_line(aes(y = Zahl, color = "Tests")) + geom_line(aes(y = Infekt_7, color = "Infektionen")) + geom_line(aes(y = PosAnteil, color = "Positivrate"))
t_plot <- t_plot + theme_classic() + ylab("Normalisierte Werte - Wertebereich 0 bis 1") + ggtitle("Infektionen, Positivrate, Tests März 2020 bis April 2021") + theme(plot.title = element_text(hjust = 0.5), legend.title = element_blank())
t_plot <- t_plot + scale_color_manual(values=c("purple", "magenta", "cyan"))
t_plot
#Zeitbereich Juli bis Oktober
t_plot2 <- ggplot(data= subset(total, total$Datum > "2020-07-01" & total$Datum < "2020-09-30"), aes(x=Datum)) + geom_line(aes(y = Zahl, color = "Tests")) + geom_line(aes(y = Infekt_7, color = "Infektionen")) + geom_line(aes(y = PosAnteil, color = "Positivrate"))
t_plot2 <- t_plot2 + theme_classic() + ylab("Normalisierte Werte - Wertebereich 0 bis 1") + ggtitle("Infektionen, Positivrate, Tests Juli 2020 bis September 2020") + theme(plot.title = element_text(hjust = 0.5), legend.title = element_blank())
t_plot2 <- t_plot2 + scale_color_manual(values=c("purple", "magenta", "cyan"))
t_plot2
#Zeitbereich Oktober bis Februar
rect_1 <- geom_rect(aes(xmin=as.Date("2020-11-08"), xmax=as.Date("2020-12-02"), ymin=0, ymax=Inf), fill = "floralwhite", alpha = 0.1)
rect_2 <- geom_rect(aes(xmin=as.Date("2020-12-23"), xmax=as.Date("2021-01-09"), ymin=0, ymax=Inf), fill = "honeydew", alpha = 0.1)
t_plot3 <- ggplot(data= subset(total, total$Datum > "2020-10-01" & total$Datum < "2021-02-01"), aes(x=Datum)) + rect_1 + rect_2 + geom_line(aes(y = Zahl, color = "Tests")) + geom_line(aes(y = Infekt_7, color = "Infektionen")) + geom_line(aes(y = PosAnteil, color = "Positivrate"))
t_plot3 <- t_plot3 + theme_classic() + ylab("Normalisierte Werte - Wertebereich 0 bis 1") + ggtitle("Infektionen, Positivrate, Tests Oktober 2020 bis Februar 2021") + theme(plot.title = element_text(hjust = 0.5), legend.title = element_blank())
t_plot3 <- t_plot3 + scale_color_manual(values=c("purple", "magenta", "cyan"))
t_plot3
#Tode, Positivanteil, COVID-19 Patienten etc.
d_plot <- ggplot(data= subset(total, total$Datum > "2020-03-07" & total$Datum < "2021-04-18"), aes(x=Datum)) + rect_1 + rect_2 + geom_line(aes(y = PosAnteil, color = "Positivrate")) + geom_line(aes(y = Tod_7, color = "Tode")) + geom_line(aes(y = covCount, color = "COVID-19 Patienten"))
d_plot <- d_plot + theme_classic() + ylab("Normalisierte Werte - Wertebereich 0 bis 1") + ggtitle("COVID-19 Patienten, Positivrate, Tode März 2020 bis April 2021") + theme(plot.title = element_text(hjust = 0.5), legend.title = element_blank())
d_plot <- d_plot + scale_color_manual(values=c("red", "magenta", "black"))
d_plot 
Denkbär 