Corona-Zahlen auswerten und visualisieren ​

Inhaltsverzeichnis

Nachrichten, Sondersendungen, Extra-Beiträge – das Corona-Virus ist das beherrschende Thema in den Medien. Dort kommen auch Forscher zu Wort, die zwar die Gefahren des Virus detailreich beschreiben, der Frage nach einer Lösung und der Zeitspanne bis zur Lösung aber ausweichen. Glücklicherweise können Sie die Lücke mit ein paar Zeilen Python-Code einfach selbst schließen.

Das zugrunde liegende Modell ist sehr grob, die Vorhersagen daher vage. Für eine grobe Einschätzung, wie lange die Politik gezwungen sein wird, die Ausgangs­beschränkungen in Deutschland bestehen zu lassen, reicht es aber aus.

Die Anzahl der Mitmenschen, die ein Infizierter im Schnitt ansteckt – die sogenannte Basisreproduktionszahl R₀ des Virus –, liegt in einer Gesellschaft wie in Deutschland (vor der Einführung der aktuellen Einschränkungen) bei ungefähr 3. Dadurch wächst die Anzahl der Erkrankten exponentiell.

Damit sich das Virus nicht mehr weiterverbreiten kann, muss R₀ unter 1 sinken. Dann würde es langsam aussterben. Denn wenn keiner immun ist, stecken sich 3 Leute an. Wenn 2 von 3 immun sind, steckt sich nur noch einer an.

Impfstoff das beste Mittel

Zwei Drittel Immune in der Bevölkerung erreicht am schnellsten ein Impfstoff. Es wird zwar an mehreren geforscht, keines der Teams glaubt aber an eine Massenproduktion vor Frühjahr 2021. Alternativ gibt es irgendwann genügend Menschen, die eine Covid-19-Infektion überstanden haben. Deren Zahl wächst aber bei einer Infektionsrate, die für das Gesundheitssystem verkraftbar ist, nur sehr langsam.

Eine direkte Möglichkeit, um R₀ unter 1 zu drücken, besteht darin, die Kontakte von Infizierten mit Mitmenschen so weit zu beschränken, dass sich nur noch wenige anstecken können. Das versucht die Politik mit den Kontaktsperren zu erreichen. Dabei nimmt man in Kauf, dass die Zahl der genesenen und damit immunen Menschen deutlich langsamer wächst, als für eine ausreichende „Herdenimmunität“ nötig wäre. Letztlich dienen die derzeitigen Maßnahmen also dazu, auf einen Impfstoff zu warten. Würde man die Kontaktsperren nämlich lockern, würde R₀ sofort wieder auf etwa 3 steigen und es gäbe eine zweite Infektionswelle.

Das hier vorgestellte Rechenmodell benutzt aktuelle Infektionszahlen, um eine Prognose über den weiteren Verlauf der Epidemie zu berechnen und zu visualisieren. Sie können es mit jeder aktuellen Python-Distribution (Beispielsweise Anaconda unter Windows) nachvollziehen. Den Code (ein Jupyter-Notebook) bekommen Sie aus unserem GitHub-Repository. Alle nötigen Bibliotheken wie SciPy, Pandas und Altair installieren Sie mit pip. Für macOS und Linux reichen die folgenden Befehle auf der Konsole:

git clone https://github.com/pinae/Covid-Predict.git
cd Covid-Predict
python3 -m venv env
source env/bin/activate
pip install -U pip wheel
pip install -r requirements.txt

Unter Windows aktivieren Sie das Vir­tualenv grafisch über Anaconda, öffnen über das Menü ein Konsolenfenster und tippen dort die letzten beiden Befehle ein. Den lokalen Webserver für das Jupyter-Notebook startet jupyter notebook.

Daten sammeln

Für eine Abschätzung, wie sich die Infektionszahlen weiterentwickeln, kann man die bisherigen Werte extrapolieren. Dafür haben wir Zahlen für die Corona-Infektionen in Deutschland von worldometers.info abgetippt. Pandas lädt die Daten mit einer Zeile, eine weitere konvertiert die Datumsangaben:

data = pd.read_csv("corona_infections.csv", header=0, names=["day", "cases"])
data["day"] = [dt.datetime.strptime(d + ' 2020', "%b %d %Y") for d in data["day"]]

Die CSV-Datei enthält lediglich in der ersten Spalte eine Datumsangabe wie "Mar 20" und nach einem Komma in der zweiten Spalte die Zahl der Infizierten. Der zweite Befehl des obigen Codes konvertiert die kurzen Datumsangaben in datetime-Objekte, womit später die X-Achse im Diagramm sinnvoll beschriftet wird.

Wachstum modellieren

Leider enthält der Datensatz bisher nur wenige verwertbare Datenpunkte. Daher ist es sinnvoll, zunächst mit einem einfachen Modell anzufangen: Ein Virus trifft zunächst auf viele Wirte, wodurch es sich exponentiell verbreitet. Je mehr Wirte jedoch bereits infiziert sind, desto weniger neue Wirte findet es, sodass sich die Verbreitung später einem Maximum an­nähert. Ein solches Wachstum modelliert die logistische Funktion:

def corona_curve(x, b0, x0, k, s):
return s * 1 / (1 + np.exp(-1 * k * s * (x - x0)) * (s / b0 - 1))

Sie steigt vom Basiswert (ca. 20.000 Infizierte) b0 zunächst exponentiell an, erreicht aber einen Wendepunkt, ab dem sie abflacht und sich an ein Maximum s annähert. Eine sinnvolle obere Grenze für s ist die Bevölkerungszahl Deutschlands (ca. 83 Mio.). In der Praxis ergeben sich aber kleinere Werte, da die Kurve bei geringerer Reproduktionszahl des Virus schneller eine Sättigung erreicht.

Maßgeblich für die Steilheit der größten Steigung ist die Wachstumskonstante k. Die Bedeutung der Konstante wird klarer, wenn man ein paar Berechnungen mit ihr anstellt, die in der Formel nicht vorkommen müssen: Rechnet man e^k, kommt man von der Konstante auf den Wachstumsfaktor b. 1-b/100 gibt an, wie viele Prozentpunkte die Kurve pro Zeitschritt wächst.

Der Parameter x0 verschiebt die Kurve lediglich entlang der X-Achse, damit die Werte zum richtigen Datum passen.

Um ein Modell zu entwerfen, das beschreibt, wie sich die Infektionszahlen entwickeln werden, muss man die vier Parameter der Funktion corona_curve() so wählen, dass sie möglichst gut zu den Werten aus dem Datensatz passen. Um diese Aufgabe kümmert sich curve_fit() aus dem Python-Modul scipy.optimize. Um sie zu nutzen, brauchen Sie Daten für X und Y als Numpy-Array.

Kontaktsperren in Betracht ziehen

Seit dem 20. März gelten in allen Bundesländern weitgehende Kontaktsperren, die die Wachstumskonstante beeinflussen. Wegen der Inkubationszeit sind diese Maßnahmen in der Kurve aber durchschnittlich erst ab dem 28. März sichtbar (Variable fqd). Eine Prognose sollte sich deshalb nur auf die Zahl der nachgewiesenermaßen Erkrankten seit diesem Tag stützen. Die folgenden Zeilen extrahieren daher die Werte seit dem 28. März als Numpy-Arrays und speichern die Tag-Nummern in days_since_quarantine und die Infektionszahlen in cases_since_quarantine:

fqd = dt.datetime(year=2020, month=3, day=28)
cases_since_quarantine = np.array(data[data["day"] >= fqd]["cases"])
day_no_since_quarantine = np.array([d.toordinal() for d in data[data["day"] >= fqd]["day"]])

Die vier Parameter von corona_curve() soll curve_fit() nun automatisch so wählen, dass sie möglichst genau zum vorbereiteten Datensatz passen. Die Funktion startet dafür standardmäßig mit einem Wert von 1 für alle Parameter und arbeitet sich schrittweise an ein optimales Ergebnis heran. Falls die Exponentialfunktion bei diesem Prozess zu große Werte liefert, kommt es zu Fehlern.

Kurvenanpassung stabilisieren

Deswegen übergibt der Code im Jupyter-Notebook über den Parameter p0 zusätzlich eine Liste mit besseren Startwerten als 1, mit denen die Optimierung stabiler läuft. Eine ähnliche Hilfestellung bieten die bounds, die Listen mit Minimal- und Maximalwerten für die Parameter enthalten. xdata und ydata enthalten die Eingaben und erwarteten Ausgaben der Funktion:

params, _ = curve_fit(
corona_curve,
xdata=day_no_since_quarantine,
ydata=cases_since_quarantine,
p0=[
cases_since_quarantine[0],
dt.datetime(year=2020, month=3, day=22).toordinal(),
8e-9,
5.6e7],
bounds=(
[0, day_no_since_quarantine[0],
1e-11, cases_since_quarantine[-1]],
[cases_since_quarantine[-1], dt.datetime(year=2021, month=6, day=1).toordinal(),
1e-8, 8.35e7])
)

Die vier Parameter b0, x0, k und s gibt curve_fit() anschließend als ersten Rückgabewert params zurück (der zweite enthält Informationen zum Fehler). In dem 4-Tupel params interessieren danach vor allem der dritte und vierte Wert: Je kleiner k, desto weniger schnell steigen die Infektionszahlen. Der vierte Wert s ist mit Vorsicht zu genießen, da das Maximum der Infizierten (wer genesen ist, wird in dieser Statistik weiter mitgezählt) mit den wenigen bisher zur Verfügung stehenden Werten stark schwanken kann. Ein verlässliches s liefert SciPy wohl erst, wenn sich die Kurve dem Wendepunkt nähert.

Visualisiert

Die Kurve der tatsächlichen und der erwarteten Infektionszahlen zeigt ein Diagramm viel anschaulicher als eine Tabelle. Mittel der Wahl ist da Altair, dessen Diagramme Jupyter-Notebooks praktischerweise direkt anzeigen. Das Jupyter-Notebook in unserem Repository enthält entsprechenden Beispielcode.

Immunität statt Impfstoff

Statt auf einen Impfstoff zu warten, könnten einfach so viele Menschen Covid-19 überleben, dass irgendwann zwei Drittel auf natürlichem Weg immun werden. Experten vermuten aber, dass die Corona-Immunität wie bei der Grippe nicht ewig hält. Schätzungen reichen von sechs Monaten bis anderthalb Jahren Immunität. Um rechtzeitig zwei Drittel zu erreichen, müssten viel mehr Menschen gleichzeitig krank werden als jetzt. Für die Berechnung reicht ein Dreisatz:

bevoelkerung_de = 83019213
infizierte_pro_tag = ["{:.0f}".format(bevoelkerung_de*0.7/(365*dauer))
for dauer in [0.5, 1, 1.5]]

Ergebnis: Damit der Plan aufgeht, müssten sich im besten Fall, wenn die Immunität 18 Monate hält, täglich über hunderttausend Menschen neu anstecken.

Die hier gezeigten Modelle sind zu grob, um die Vorhersagen von Experten zu ersetzen. Wenn die aber Antworten schuldig bleiben, kann man sich jedoch mit Python selber behelfen.

Zusätzliche Links

• Zahlen zu den Corona-Infektionen in Deutschland (und anderen Ländern).
• Quelle zur Zahl der Intensivbetten für schwere COVID-19-Fälle.
• Quelle zur Anzahl der schweren Fälle, die Intensivbetten mit Beatmung benötigen.
• Epidemic Calculator. Schöne browserbasierte Modellierung anhand des Standardmodells der Epidemologie.
• Erklärung zu einem einfachen epidemologischen Modell aufbauend auf Differentialgleichungen bei YouTube.
• Epidemologisches Modell mit Daten aus Deutschland.
• Modellierung zu COVID-19 mit veränderbaren Parametern.
• Browserbasiertes Modell mit einer einfachen Exponentialfunktion.

Dieser Artikel stammt aus c't 9/2020. (pmk)