Dies ist ein Diskussionsbeitrag - denn die Ergebnisse sind sehr merkwürdig
Als Data-Scientist liegt es für mich nahe, die Daten zur Corona/COVID-19 Krise zu analysieren. Ich habe versucht, verschiedene Modelle zu fitten, um den möglichen Verlauf der Krankheit in Deutschland vorherzusagen.
Ich kenne mich aber mit Virologie überhaupt nicht aus und es kann sehr gut sein, dass die folgende Analyse dumme Fehler beinhaltet. Wäre alles richtig und würde sich die Wirklichkeit tatsächlich an das Modell halten, dann wäre der Peak der Ausbreitung in Deutschland schon in wenigen Tagen erreicht.
Die Daten
Ich verwende die Daten des Johns Hopkins University Center for Systems Science and Engineering (JHU CSSE). Die Daten werden täglich auf GitHub veröffentlicht, unter: https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv.
Exponentiell oder Sigmoid?
Häufig wird über das exponentielle Wachstum von COVID-19 gesprochen. Auch die Stellungnahme der Deutschen Gesellschaft für Epidemiologie (DGEpi) spricht davon, dass sich das Virus "anfangs mit exponentiell steigender Geschwindigkeit" (https://www.dgepi.de/assets/Stellungnahmen/Stellungnahme2020Corona_DGEpi-20200319.pdf) ausbreitet. Ich bin wie gesagt kein Experte, verstehe die Stellungnahme aber so, dass in den Modellen von einer gleichbleibenden Geschwindigkeit (gemessen in der Basisreproduktionszahl R0) ausgegangen wird. Ich gehe davon aus, dass die verwendeten Modelle dem Stand der Forschung entsprechen und theoretisch sehr gut begründet sind. Aus der Stellungnahme selbst und ohne Fachwissen lässt sich da auch gar nicht viel zu sagen.Aus der Data Science Perspektive sieht das schon anders aus: Anders als die klassische Statistik sind wir bemüht, wenige theoretische Annahmen zu machen und stattdessen induktiv möglichst viel Informationen aus echten Daten abzuleiten. Dabei zeigt sich häufig, dass Prozesse, von denen man zunächst annimmt, sie seien exponentiell, häufig einer Sigmoid-Funktion folgen. Nach dem exponentiellen Anstieg nimmt der Zuwachs genauso schnell wieder ab und das System pendelt sich auf einem neuen stabilen Wert ein. Gerade am Anfang verhalten sich Exponentielle Funktion und Sigmoid allerdings sehr ähnlich, was es schwierig macht, die richtige Funktion auszuwählen.
Entwicklung in Deutschland
Die Fallzahlen in Deutschland verändern sich jeden Tag rapide nach oben. Um Vorhersagen zu machen, braucht man ein mathematisches Modell, das die tatsächlichen Fällen einfach anhand der vergangenen Tage möglichst exakt schätzt. Als Modell verwende ich eine Exponentialfunktion und als Alternativen eine Sigmoid-Funktion mit drei und mit vier Parametern.
Es ist ganz wichtig, dass man sich klar macht, dass die Daten nicht die echte Anzahl der Infizierten anzeigen. Es sind die offiziellen Fallzahlen. Diese hängen zum Beispiel davon ab, wie viele Leute getestet werden. Bei COVID-19 hat sich die Teststrategie sehr verändert: Zunächst wurde (verkürzt gesagt) nur getestet, wer aus einem Krisengebiet gekommen ist, heute wer deutliche Symptome zeigt (und beim Arzt auftaucht). Die DGEpi sagt: "Veröffentlichte Daten sprechen
dafür, dass Infektionen mit SARS-CoV-2 in den meisten Fällen milde bis moderat verlaufen, während 2–5% der Fälle eine Behandlung auf der Intensivstation benötigen." D. h., unabhängig von der Teststrategie werden vermutlich viele Fälle gar nicht getestet. Für unser vorgehen hier ist das eigentlich gut. Wir gehen davon aus, dass die Daten durch alle möglichen Einflussfaktoren verändert sind. Gleichzeitig ist aber auch davon auszugehen, dass die schweren Fälle, um die es ja geht, irgendwann in den Daten auftauchen. Wie viele leichte Fälle es in die Statistik schaffen, wird aber vermutlich je nach Land und Lage unterschiedlich sein.
Wenn ich die drei genannten Funktionen so anpasse, dass sie die Daten möglichst gut abbilden, entsteht folgendes Bild:
Man sieht, die Sigmoidfunktionen bilden die Entwicklung in den Daten besser ab. Das liegt zunächst einfach daran, dass diese Funktionen mit mehr Parametern flexibler sind. Um Overfitting zu vermeiden, versucht man einen Kompromiss zwischen Genauigkeit und Einfachheit der Funktion zu finden. Im vorliegenden Fall spricht meiner Meinung nach sehr viel für die Sigmoid 3. Wie in den Codes am Ende zu sehen, hat sie sogar die geringste Abweichung.
Da wir jetzt Modelle mit den entsprechenden geschätzten Parametern haben können wir uns anschauen, wie sich diese Modelle in der Zukunft entwickeln. (Wohlgemerkt: Die Ergebnisse entspringen jetzt dem Modell, es ist keine Wahrsagerei, sondern die Fortschreibung der getroffenen Modellannahmen.)
Je nachdem, welches Modell wir bevorzugen fällt die zukünftige Entwicklung sehr unterschiedlich aus.
Dieses Bild enthält jetzt zwei interessante Informationen. Gehen wir davon aus, dass die echte Entwicklung durch eine Sigmoidfunktion abgebildet werden können, dann pendelt sich die Ausbreitung ziemlich schnell ein. Je nachdem, welches Sigmoidmodell wir nehmen, liegt die Schätzung der final Infizierten zwischen 100.000 und 350.000 Personen. Bei der Exponentialfunktion steigt die Zahl natürlich ins unendliche (weswegen in der echten Welt kein Wachstumsprozess exponentiell sein kann).
Wir können jetzt auch berechnen, zu welchem Zeitpunkt welches Modell einen Peak erreicht. Für die Sigmoidfunktionen ist klar, dass nach dem höchsten Zuwachs die Wachstumsrate rückläufig ist.
Bei der Exponentialfunktion steigt die Wachstumsrate beständig an. Allerdings können ja nur gesunde Leute angesteckt werden. Daher kann man davon ausgehen, dass sich auch unter exponentiellen Annahmen das Wachstum umdreht, spätestens wenn zweidrittel der Bevölkerung infiziert sind. Dieser Wert (rote gestrichelte Linie) ist aber extrem ungenau, weil ja längst nicht alle Infizierten in den Daten sind.
Entscheiden wir uns dafür, dass die Sigmoidfunktion den echten Prozess besser widerspiegelt, dann heißt das, dass sich die Fallzahlen relativ schnell stabilisieren würden, nämlich bereits in den nächsten zwei Wochen.
Ob das das Resultat der politische Intervention wäre, oder durch einen rein biologischen Vorgang erklärt werden könnte, entzieht sich meiner Kenntnis.
Warum Sigmoid: China!
Richtig spannend wird es, wenn man die Entwicklung in Deutschland jetzt mit der chinesischen Provinz Hubei vergleicht.
Hier sehen wir in den echten Daten die Form der Sigmoid-Kurve! Und natürlich passen dann die Modelle auch am besten.
Das krasse ist, dass die Parameter der Sigmoidfunktionen sich gar nicht groß vom deutschen Fall unterscheiden. Und auch der Breakpoint, wenn die Infektionen wieder zurückgehen, liegt bei 62 Tagen, also in dem Bereich, den auch die Modelle auf den deutschen Daten vorhersagen (mit größerer Unsicherheit, weil die Entwicklung noch nicht soweit fortgeschritten ist).
Ist die Krise bald vorbei?
Ich habe keine Ahnung! Ich kann nur sagen, dass das stumpfe Schauen auf die Daten diesen Eindruck vermittelt. Noch weniger weiß ich, woran das liegen könnte, falls es denn stimmt.
Mit den vorliegenden Daten kann ich auch nicht einschätzen, ob ein Anstieg der bestätigten Fälle auf 350.000 einen Kollaps des Gesundheitssystems bedeutet. Vermutlich schon, denn wenn wir davon ausgehen, dass die Krankheit in 2-5% der Fälle behandelt werden muss, dass vermutlich mehr oder weniger alle schweren Fälle auch bei den offiziell gemeldeten auftauchen, also die Zahl der zu behandelnden Fälle auch bei 10% der 350.000 liegen könnte und 30.000 Krankenhausplätze vorhanden sind (die nicht frei sind!), dann klingt die Situation immer noch dramatisch. Allerdings sagt das beste Modell eine offizielle Fallzahl von 100.000 voraus. Auch 10% davon wären vielleicht zu verkraften.
Wenn diese Analyse einen Wert hat, dann müssten wir schon in der nächsten Woche einen Rückgang der Wachstumsrate sehen. Es bleibt also spannend.
Wenn diese Analyse einen Wert hat, dann müssten wir schon in der nächsten Woche einen Rückgang der Wachstumsrate sehen. Es bleibt also spannend.
Ich freue mich sehr, wenn Leute, die sich auskennen, diesen Ansatz widerlegen! Ich kann nur wiederholen: Ich schaue auf die Daten! Von Viren und Seuchenpolitik verstehe ich nichts.
Code
# get data from John Hopkins University
df <- read.csv("https://github.com/CSSEGISandData/COVID-19/raw/master/csse_covid_19_data/csse_covid_19_time_series/time_series_19-covid-Confirmed.csv", stringsAsFactors = F)
ge <- as.numeric(df[df$Country.Region=="Germany",5:dim(df)[2]])
day <- c(1:length(ge))
dfGE <- data.frame(day,ge)
dfGE <- dfGE[dfGE$ge>0,]
fit1 <- lm(log(ge) ~ day, data = dfGE)
summary(fit1)##
## Call:
## lm(formula = log(ge) ~ day, data = dfGE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8997 -0.6296 0.2460 0.7800 1.3887
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.616371 0.298312 -2.066 0.0439 *
## day 0.152837 0.008411 18.172 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9366 on 51 degrees of freedom
## Multiple R-squared: 0.8662, Adjusted R-squared: 0.8636
## F-statistic: 330.2 on 1 and 51 DF, p-value: < 2.2e-16timevalues <- seq(1, max(dfGE$day), 0.1)
predCurve <- exp(predict(fit1,list(day=timevalues)))
plot(dfGE$ge~dfGE$day,pch=16,
main = "COVID-Models - Germany",
xlab = "Days", ylab="Number infected")
lines(timevalues, predCurve,lwd=2, col = "red", xlab = "Day", ylab = "Counts")
MAE <- function(x, p){
z <- mean(abs(p - x), na.rm = TRUE)
return(z)
}
pred1 <- exp(predict(fit1))
MAE1 <- MAE(dfGE$ge, pred1)
# Verbreitung wird langsamer, wenn die 66% der Bevölkerung infiziert sind
oneYear <- exp(predict(fit1, list(day=c(1:365))))
breakpoint1 <- as.numeric(which(oneYear>=82000000*0.66)[1])
library(drc)## Loading required package: MASS##
## 'drc' has been loaded.## Please cite R and 'drc' if used for a publication,## for references type 'citation()' and 'citation('drc')'.##
## Attaching package: 'drc'## The following objects are masked from 'package:stats':
##
## gaussian, getInitialfit2 <- drm(ge ~ day, fct = L.3(), data = dfGE)
summary(fit2)##
## Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -2.7099e-01 7.4787e-03 -36.2353 < 2.2e-16 ***
## d:(Intercept) 1.1542e+05 4.2155e+04 2.7381 0.008539 **
## e:(Intercept) 6.4911e+01 1.7206e+00 37.7265 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 108.583 (50 degrees of freedom)predCurve <- predict(fit2, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
lines(timevalues, predCurve,lwd=2, col = "green")
MAE2 <- MAE(dfGE$ge, predict(fit2))
timevalues <- c(1:365)
oneYear <- predict(fit2, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
oneYearDif <- oneYear[-1]-oneYear[-365]
breakpoint2 <- which.max(oneYearDif)
fit3 <- drm(ge ~ day, data = dfGE, fct = LL.4(c(fixed = NA, NA, NA, NA), names = c("Slope", "Lower Limit", "Upper Limit", "ED50")))
summary(fit3)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## Slope:(Intercept) -1.4100e+01 1.7458e-01 -80.7656 < 2.2e-16 ***
## Lower Limit:(Intercept) 2.5670e+01 1.8432e+01 1.3927 0.169995
## Upper Limit:(Intercept) 3.6170e+05 1.1408e+05 3.1705 0.002624 **
## ED50:(Intercept) 7.2369e+01 1.7963e+00 40.2888 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 114.6457 (49 degrees of freedom)predCurve <- predict(fit3, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
lines(timevalues, predCurve,lwd=2, col = "blue")
legend("topleft", legend = c("exponential", "sigmoid 3", "sigmoid 4"), col = c("red", "green", "blue"), lty = 1)
MAE3 <- MAE(dfGE$ge, predict(fit3))
oneYear <- predict(fit3, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
oneYearDif <- oneYear[-1]-oneYear[-365]
breakpoint3 <- which.max(oneYearDif)
cat("MAE fit1: ", MAE1)## MAE fit1: 851.5581cat("MAE fit2: ", MAE2)## MAE fit2: 53.07257cat("MAE fit3: ", MAE3)## MAE fit3: 55.95039# the future
timevalues <- seq(1, 130, 0.1)
predCurve <- exp(predict(fit1,list(day=timevalues)))
plot(dfGE$ge~dfGE$day,pch=16, xlim=c(0,130), ylim=c(0,500000), main = "COVID Pseudo-Forecast - Germany",
sub= "breakpoints in dashed lines", xlab = "Days", ylab="Number infected")
lines(timevalues, predCurve,lwd=2, col = "red", xlab = "Day", ylab = "Counts")
predCurve <- predict(fit2, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
lines(timevalues, predCurve,lwd=2, col = "green")
predCurve <- predict(fit3, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
lines(timevalues, predCurve,lwd=2, col = "blue")
abline(v=breakpoint1, col="red",lty=2)
abline(v=breakpoint2, col="green",lty=2)
abline(v=breakpoint3, col="blue",lty=2)
legend("topleft", legend = c("exponential", "sigmoid 3", "sigmoid 4"), col = c("red", "green", "blue"), lty = 1)
hb <- as.numeric(df[df$Province.State=="Hubei",5:dim(df)[2]])
# ge <- as.numeric(df[df$Country.Region=="Germany",5:dim(df)[2]])
day <- c(44:(length(hb)+43))
dfHB <- data.frame(day,hb)
dfHB <- dfHB[dfHB$hb>0,]
fit1hb <- lm(log(hb) ~ day, data = dfHB)
summary(fit1)##
## Call:
## lm(formula = log(ge) ~ day, data = dfGE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8997 -0.6296 0.2460 0.7800 1.3887
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.616371 0.298312 -2.066 0.0439 *
## day 0.152837 0.008411 18.172 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9366 on 51 degrees of freedom
## Multiple R-squared: 0.8662, Adjusted R-squared: 0.8636
## F-statistic: 330.2 on 1 and 51 DF, p-value: < 2.2e-16timevalues <- seq(1, max(dfHB$day), 0.1)
predCurve <- exp(predict(fit1hb,list(day=timevalues)))
plot(dfHB$hb~dfHB$day,pch=16,
main = "COVID-Models - Hubei",
xlab = "Days", ylab="Number infected")
lines(timevalues, predCurve,lwd=2, col = "red", xlab = "Day", ylab = "Counts")
pred1hb <- exp(predict(fit1hb))
MAE1hb <- MAE(dfHB$hb, pred1hb)
# Verbreitung wird langsamer, wenn die 66% der Bevölkerung infiziert sind
oneYear <- exp(predict(fit1hb, list(day=c(1:365))))
breakpoint1hb <- as.numeric(which(oneYear>=82000000*0.66)[1])
fit2hb <- drm(hb ~ day, fct = L.3(), data = dfHB)
summary(fit2hb)##
## Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -2.3491e-01 8.1363e-03 -28.872 < 2.2e-16 ***
## d:(Intercept) 6.7714e+04 4.2569e+02 159.068 < 2.2e-16 ***
## e:(Intercept) 6.2689e+01 1.7350e-01 361.315 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 2073.187 (55 degrees of freedom)predCurve <- predict(fit2hb, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
lines(timevalues, predCurve,lwd=2, col = "green")
MAE2hb <- MAE(dfHB$hb, predict(fit2hb))
timevalues <- c(1:365)
oneYear <- predict(fit2hb, newdata = data.frame(ge=rep(NA,length(timevalues)),day = timevalues))
oneYearDif <- oneYear[-1]-oneYear[-365]
breakpoint2hb <- which.max(oneYearDif)
fit3hb <- drm(hb ~ day, data = dfHB, fct = LL.4(c(fixed = NA, NA, NA, NA), names = c("Slope", "Lower Limit", "Upper Limit", "ED50")))
summary(fit3hb)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## Slope:(Intercept) -14.82929 0.74486 -19.9088 <2e-16 ***
## Lower Limit:(Intercept) 768.36095 923.40740 0.8321 0.409
## Upper Limit:(Intercept) 68099.49071 487.93534 139.5666 <2e-16 ***
## ED50:(Intercept) 62.71719 0.24551 255.4615 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 2189.899 (54 degrees of freedom)predCurve <- predict(fit3hb, newdata = data.frame(per=rep(NA,length(timevalues)),day = timevalues))
lines(timevalues, predCurve,lwd=2, col = "blue")
legend("topleft", legend = c("exponential", "sigmoid 3", "sigmoid 4"), col = c("red", "green", "blue"), lty = 1)
MAE3hb <- MAE(dfHB$hb, predict(fit3hb))
oneYear <- predict(fit3hb, newdata = data.frame(per=rep(NA,length(timevalues)),day = timevalues))
oneYearDif <- oneYear[-1]-oneYear[-365]
breakpoint3hb <- which.max(oneYearDif)
abline(v=breakpoint1hb, col="red",lty=2)
abline(v=breakpoint2hb, col="green",lty=2)
abline(v=breakpoint3hb, col="blue",lty=2)
summary(fit2)##
## Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -2.7099e-01 7.4787e-03 -36.2353 < 2.2e-16 ***
## d:(Intercept) 1.1542e+05 4.2155e+04 2.7381 0.008539 **
## e:(Intercept) 6.4911e+01 1.7206e+00 37.7265 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 108.583 (50 degrees of freedom)summary(fit2hb)##
## Model fitted: Logistic (ED50 as parameter) with lower limit fixed at 0 (3 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## b:(Intercept) -2.3491e-01 8.1363e-03 -28.872 < 2.2e-16 ***
## d:(Intercept) 6.7714e+04 4.2569e+02 159.068 < 2.2e-16 ***
## e:(Intercept) 6.2689e+01 1.7350e-01 361.315 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 2073.187 (55 degrees of freedom)summary(fit3)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## Slope:(Intercept) -1.4100e+01 1.7458e-01 -80.7656 < 2.2e-16 ***
## Lower Limit:(Intercept) 2.5670e+01 1.8432e+01 1.3927 0.169995
## Upper Limit:(Intercept) 3.6170e+05 1.1408e+05 3.1705 0.002624 **
## ED50:(Intercept) 7.2369e+01 1.7963e+00 40.2888 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 114.6457 (49 degrees of freedom)summary(fit3hb)##
## Model fitted: Log-logistic (ED50 as parameter) (4 parms)
##
## Parameter estimates:
##
## Estimate Std. Error t-value p-value
## Slope:(Intercept) -14.82929 0.74486 -19.9088 <2e-16 ***
## Lower Limit:(Intercept) 768.36095 923.40740 0.8321 0.409
## Upper Limit:(Intercept) 68099.49071 487.93534 139.5666 <2e-16 ***
## ED50:(Intercept) 62.71719 0.24551 255.4615 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error:
##
## 2189.899 (54 degrees of freedom)cat("Day Peak reached, model 2: ", breakpoint2)## Day Peak reached, model 2: 64cat("Day Peak reached, model 2, Hubei: ", breakpoint2hb)## Day Peak reached, model 2, Hubei: 62cat("Day Peak reached, model 3: ", breakpoint3)## Day Peak reached, model 3: 71cat("Day Peak reached, model 3, Hubei: ", breakpoint3hb)## Day Peak reached, model 3, Hubei: 62
Kommentare
Kommentar veröffentlichen