Halbwegs gute Nachrichten

Zeitreihen

Wie habe ich das gerechnet? Nun, die Korrelation von Zeitreihen ist ein ganzer Satz von Wissenschaften, bei denen es meist darum geht, ungleiche Abdeckungen, unregelmäßig gesetzte Messpunkte sowie allerlei Rauschen und Schmutz weggefummelt zu kriegen, ohne allzu viele Informationen zu verlieren oder, vielleicht schlimmer, damit Artefakte einzubauen.

Die so gereinigten Daten lassen sich dann zum Beispiel geeignet skaliert und verschoben übereinanderlegen. Tatsächlich sind die Parameter dieser Transformationen im Regelfall (und gewissermaßen auch hier) viel interessanter als die Zeitreihen selbst[1]. Es gibt daher zahlreiche mathematische Verfahren, die das von einer Augenmaß-Übung in etwas verwandeln, das reproduzierbar und auch quantifizierbar ist. Vorsicht: ich bin da kein Experte und gehe hier nur mit nicht allzu schwer erkranktem Menschenverstand ran, verwende also (zumindest in Summe) gerade kein wirklich wohldurchdachtes Verfahren. Wer sowas wie das hier für Hausaufgaben oder Hausarbeiten verwendet, tut das auf eigene Gefahr.

Andererseits bin ich recht zuversichtlich, dass der Kram insgesamt schon stimmt.

Ausgangsdaten sind meine aus RKI-Berichten gescrapten Intensivbelegungen und ein RKI-Sheet (in, ach weh, „Office Open XML“ a.k.a. XSLX – muss das sein?) mit den Inzidenzen. Dabei ist das erste Problem, dass ich aus verschiedenen Gründen nicht für jeden Tag Belegungszahlen habe, und so ist mein erster Schritt, fehlende Punkte zu interpolieren. Das Scipy-Paket macht es leicht, aus einem Satz von Zeit/Wert-Paaren (die Zeit wird in den Lesefunktionen auf Sekunden seit einer Epoche gewandelt) ein handliches Array zu rechnen:

grid_points = numpy.arange(
  raw_crits[0][0],
  raw_crits[-1][0],
  86400)  # ein Tag in Sekunden
critnums = interpolate.griddata(
  raw_crits[:,0], raw_crits[:,1], grid_points)

Ich wollte diese Interpolation eigentlich visualisieren, habe aber keine gute Stelle gefunden, an der sie einen sichtbaren Unterschied gemacht hätte, und entscheidend ist eigentlich nur, dass ich ab diesem Schritt blind mit Arrays arbeiten kann; deren Index ist zunächst die Zahl der Tage seit dem ersten Tag mit Daten, hier speziell dem 11.8.2020, denn damals habe ich mit dem Screenscrapen der DIVI-Daten angefangen.

Übergeplottete Kurven

Wenn ich diese interpolierten Kurven übereinanderlege, kommt das hier heraus:

Weil das, wie gesagt, erst im August 2020 anfängt, fehlt die erste Welle.

Das bloße Auge reicht für die Bestätigung der Erwartung, dass die Intensivbelegung der Inzidenzkurve meist etwas hinterherläuft, wenn auch nicht so, dass mensch hoffen könnte, die beiden durch etwas Schieben global übereinanderzubekommen. Die großen Zacken in den Inzidenzen durch Weihnachten und Ostern finden sich in der Intensivbelegung gar nicht wieder, was ein klares Zeichen ist, dass sie weitgehend Erfassungsartefakte sind. In der Hinsicht wären die „großen“ RKI-Zahlen mit Referenzdaten (vgl. die Film-Geschichte) bestimmt besser, aber ich wollte für dieses Ding nicht die 200 Megabyte durchkämmen, zumal die „kleinen“ RKI-Daten besser auf die Meldedaten aus den Tagesberichten passen, um die es mir hier ja geht.

Vor allem fällt auf, dass meine Jammerei darüber, dass der Impffortschritt die Intensivantwort nicht wesentlich abgeflacht hat, unzutreffend ist: Mit der Skalierung aus dem Plot liegen Inzidenz und Belegung in der zweiten und dritten Welle ziemlich übereinander, während in der vierten Welle doch ein knapper Faktor zwei dazwischenliegt; da scheint die Impfung doch ein wenig gegenüber Delta zu gewinnen (aber, klar, nirgendwo hinreichend).

Ich hatte erwartet, dass die abfallenden Flanken der Intensivkurven deutlich flacher sind als die der Inzidenzkurven, weil Leute unter Umständen lang auf Intensiv liegen und es entsprechend lang dauern sollte, bis sich die Stationen wieder leeren. Das sieht im Abfall der zweiten Welle auch ein wenig so aus, nicht jedoch bei dem der dritten Welle. Bei ihr fällt die Intensivbelegung sehr treu mit der Inzidenz. So viele LangzeitpatientInnen gibt es glücklicherweise wohl doch nicht.

Nach Glätten differenzierbar

Wie kann ich jetzt den Verzug zu quantifizieren? Mein (wie gesagt eher intuitiv gefasster) Plan ist, über die Ableitung der jeweiligen Kurven zu gehen, und zwar aus der Überlegung heraus, dass mich ja Veränderungen viel mehr als Pegel interessieren. Nun habe ich aber keine differenzierbaren Funktionen, sondern lediglich Arrays, bei denen ich Ableitungen allenfalls durch Subtrahieren benachbarter Elemente simulieren kann. Solche numerischen „Ableitungen“ reagieren ziemlich empfindlich auf das Gewackel („Rauschen“), das es in realen Daten immer gibt („Subtraktion ist in der Regel numerisch schlecht konditioniert“). Deshalb will ich meine rohen Daten glätten, bevor ich die „Ableitung“ ausrechne, sprich: das Rauschen rausnehmen, ohne das Signal wesentlich zu verzerren.

Glättung heißt eigentlich immer, Kurvenpunkte in einem Zelle für Zelle über die Daten laufenden Fenster zu mitteln, also z.B., indem mensch je fünf Nachbarpunkte rechts und links auf den aktuellen Punkt draufaddiert und ihn dann durch die durch elf geteilte Summe ersetzt. Mit diesem ganz naiven Rezept bekommen relativ weit entfernte Werte aber genauso viel Einfluss auf den aktuellen Punkt wie die unmittelbaren Nachbarn. Das modelliert meist die sachlichen Grundlagen des Rauschens nur schlecht. Es ist auch aus theoretischeren Gründen normalerweise eher ungünstig.

Der eher theoretische Hintergrund ist grob, dass bei der Glättung letztlich zwei Funktionen gefaltet werden. Das ist äquivalent dazu, die Spektren der beiden Signale (ihre Fouriertransformierten) zu multiplizieren und dann wieder zurückzutransformieren. Bei einem einfachen Fenster des oben skizierten Typs („Rechteckfunktion“) gibts nun sehr scharfe Kanten, die wiederum zu einem sehr breiten Spektrum führen, so dass auch das Spektrum der geglätteten Funktion allerlei unwillkommene Features bekommen kann (manchmal ist das aber auch genau das, was mensch haben will – hier nicht).

Wie auch immer: scipy macht es einfach, mit einer lärmarmen Funktion – die sieht ein wenig gaußig aus und ist im nächsten Plot im kleinen Inset zu sehen – zu glätten, nämlich etwa so:

import numpy
from scipy import signal

smoothing_kernel = numpy.kaiser(20, smoothing_width)
smoothing_kernel = smoothing_kernel/sum(smoothing_kernel)
convolved = signal.convolve(arr, smoothing_kernel, mode="same"
  )[smoothing_width:-smoothing_width]

Die Division durch die Summe von smoothing_kernel in der zweiten Zeile sorgt dafür, dass sich an der „Höhe“ der geglätteten Funktion insgesamt nichts ändert: Sozusagen ausgeklammert ist die Faltung eine Multiplikation mit eins. Das Wegschneiden der Ränder in der letzten Zeile wiederum entfernt Punkte, bei denen fehlende Werte am Anfang und Ende der Zeitreihe im Fenster waren. Scipy ersetzt die durch Nullen, so dass die geglätteten Werte am Rand steil abfallen. Was ich hier mache, entspricht technisch dem valid-Mode der convolve-Funktion. Nur weiß ich hier zuverlässig, wie lang das Array am Schluss ist.

Der Effekt, hier auf Indzidenzdaten der vierten Welle:

Damit kann ich jetzt meine numerische „Ableitung“ bilden, ohne dass mir das Ergebnis furchtbar rauscht:

diff = convolved[1:]-convolved[:-1]

Die nächsten beiden Grafiken zeigen diese Pseudo-Ableitungen für Inzidenz und Intensivbelegung. Ich habe jeweils in blau reingemalt, wie es ohne Glättung aussehen würde, um deutlich zu machen, warum diese eine gute Idee ist und was ich mit „furchtbar rauschen“ meine. Weiterverwendet werden natürlich die orangen Verläufe:

In der Inzidenzkurve fallen wieder Weihnachten und Ostern besonders auf, weil sie selbst in den geglätteten Graphen noch wilde Ausschläge verursachen.

In der Zeit verschieben

Schon der optische Eindruck aus dem Rohdaten-Plot legt nahe, die einzelnen Wellen getrennt zu untersuchen, und das ist angesichts von über die Zeit stark veränderlichen Viren, Testverhältnissen, demographischen Gegebenheiten und nichtpharamzeutischen Maßnahmen sicher auch sachlich geboten.

Deshalb hier zunächst der Verlauf der Ableitungen von Inzidenzen und Intensivbelegung in der aktuellen vierten Welle (das ist eine Kombination der rechten Enden der orangen Graphen der letzten beiden Plots):

In dem Graphen steckt noch eine weitere Normalisierung. Und zwar habe ich für beide Arrays etwas wie:

incs /= sum(abs(incs))

laufen lassen. Damit ist die Fläche zwischen den Kurven und der Abszisse in beien Fällen (etwa) gleich groß (nämlich eins), obwohl die absoluten Änderungen bei den Intensivbelegungen etwa eine Größenordnung größer sind als bei den Inzidenzpunkten. Ziel ist, dass ich nur die Absolutwerte der Differenzen der jeweiligen Komponenten zusammenaddieren muss, um ein relativ robustes Maß für den Unterschiedlichkeit der beiden Graphen zu haben (nämlich ungefähr die Fläche zwischen den beiden Linien). Ich hätte auch z.B. auf gleiche Minimal- und Maximalwerte normieren können, aber das ist meist weniger stabil gegen noch vorhandenes Rauschen.

Ach so: In Wirklichkeit rechne ich hier mit Quadrat statt Betrag („euklidische Distanz“), was wiederum eher theoretische Vorteile (z.B. Differenzierbarkeit bei Null) hat, aber dafür nicht ganz so direkt als Fläche interpretierbar ist.

Jetzt kommt jedenfalls der eigentliche Trick. Ich verschiebe nämlich die Kurven zeitlich („in x-Richtung“) gegeneinander und berechne dann mein Unterschiedsmaß. Das mache ich für Verschiebungen („offsets“) zwischen einem und dreißig Tagen und sehe nach, wann die Differenz minimal wird. Übereinanderlegen und Residuenbildung („Residuum“ ist Jargon für „was halt so an Unterschied übrigbleibt“) sehen in Code (konzeptionell; in Wirklichkeit passiert das in einer List Comprehension) so aus:

diffs = incs[:-offset]-crits[offset:]
resid = sum(diffs**2)/len(diffs)

– durch die Länge muss ich dividieren, weil der Überlapp mit größerem Offset kleiner wird und weniger Punkte in diesem Spiel sonst ein Vorteil wären.

Auswerten

Und damit sind wir fertig. Plottet mensch diese Residuen über den Offsets, kommt die Kurve ganz am Anfang des Artikels raus:

Das Minimum dieser Kurve ist ein recht robustes Maß für den Verzug von Inzidenzsignal und Intensivantwort, und das lässt sich hier als zwischen acht und neun Tagen ablesen. Am Freitag lagen bei einer Inzidenz von 170 laut DIVI 2332 SARS-2-PatientInnen auf Intensivstationen. Neun Tage davor war die Inzidenz 118, ist also in der Zwischenzeit um 44% angestiegen. Diesen Anstieg müssten wir bis in neun Tagen auch auf den Intensivstationen sehen, ganz egal, was jetzt noch passiert, denn die Leute, die dort landen werden, sind schon krank. Die Vorhersage wäre nur dann falsch, wenn sich plötzlich massiv etwas an der Demographie der Kranken und damit ihrer mittleren Krankheitsschwere geändert hätte, was mitten in einer Welle zunächst nicht zu erwarten ist.

Also: ich würde relativ hohe Wetten eingehen, dass um den 15.11. herum knapp 3500 SARS-2-PatientInnen auf deutschen Intensivstationen liegen werden.

Das sind im Hinblick auf Möglichkeiten, noch etwas zu machen, gute Nachrichten, denn wenn die Antwort wirklich 21 Tage später käme, würden wir auf den 15.10. zurückschauen. Seitdem hat sich die Inzidenz um 250% erhöht, und wir müssten in 21 Tagen fest mit knapp 6000 SARS-2-PatientInnen auf Intensiv rechnen, was schon knapp würde. So, wie es in Wirklichkeit ist, können wir das noch einbremsen (wenn wir wollen).

Ich habe die Rechnung auch für die beiden Wellen vorher machen lassen. In der dritten Welle sahen die normierten Ableitungen so aus:

Die Osterfeiertage sehen hier ganz schlimm aus. Entsprechend sind die Residuen absolut eine Ecke größer, auch wenn die generelle Form der Kurve der der vierten Welle erstaunlich ähnlich sieht:

Hier ist der Intensivantwort eher nur sieben Tage verzögert, und zumindest von der Anmutung der Kurve her würde ich das trotz des Osterzackens für eine recht realitätsnahe Schätzung halten.

In der zweiten Welle sieht es fast bilderbuchartig aus, weil ich vor den Weihnachtsferien aufhöre:

Die Residuen hier:

Die Intensivantwort kam hier also nochmal etwas schneller, braucht nämlich nur gut fünf Tage.

Ob die scheinbar immer langsamer werdende Intensivantwort viel zu bedeuten hat, kann ich überhaupt nicht einschätzen; dass sie ein echter Effekt in den Daten ist, ist aber glaube ich nicht bestreitbar. Da ich hier ja mit Meldedaten arbeite, könnte der Grund ganz schlicht sein, dass die Meldeketten im Laufe der Pandemie schneller geworden sind. Wenn nämlich Menschen im Schnitt 14 Tage nach der Infektion auf Intensiv kämen, wäre der beobachtete Effekt schlicht erklärbar, wenn in der zweiten Welle im Schnitt neun Tage nach Infektion gemeldet worden wäre (wegen Testverzug, Faxstau usf) und inzwischen so schnell getestet und berichtet wird, dass zwischen Infektion und Meldung nur noch fünf oder sechs Tage vergehen.

Es kann aber auch sein, dass ein demographischer Effekt dabei ist: Vielleicht dauert es bei jüngeren Menschen einfach etwas länger, bis sie intensivpflichtig werden? Meine Daten lassen da wohl keine Schlüsse zu.

Ein wenig Python

Das Programm, das das alles rechnet und malt, ist nicht wirklich spektakulär und nicht sehr flexibel geschrieben. Es hat insbesondere einen dramatischen Tiefpunkt beim Parsen des XLSX vom RKI – gibt es eine nicht-durchgeknallte Art, solche Spreadsheets maschinell weiterzuverarbeiten? Das ist durchaus ernstgemeint: Ich verstehe nicht viel von Office-Software, und vielleicht übersehe ich ja, wie in dem Geschäft ordentlicher Markup simuliert wird.

Nein: ich erwarte das nicht wirklich. Aber wenn die halbe Welt mit dem Zeug rummurkst, bin ich ja sicher nicht der erste, der da recht sinnlos über die Zeilen iteriert und freihändig rät, was wohl wo sein könnte.

Gut (wenn ich das selbst sagen darf) gefällt mir hingegen dieser Context Manager, der bei der großen Plotorgie m.E. sehr hilft:

@contextlib.contextmanager
def fig_writeto(dest_name):
  """starts a new pyplot figure and saves it to dest_name once
  the controlled section is left.
  """
  pyplot.figure(figsize=(9,3))
  pyplot.axes()
  yield
  pyplot.legend()
  pyplot.savefig(dest_name)
  pyplot.close()

Damit ist ein Plot im Code mit Speichern und allem so kompakt:

with fig_writeto("crit-and-inc.svg"):
  pyplot.plot(critnums/30, label="Intensivbelegung/30")
  pyplot.plot(incs, label="Inzidenz")

Auf dieser Basis ließe sich glatt eine nicht furchtbare Approximation der ja nach heutigen Maßstäben doch etwas arg angestaubten pyplot-API basteln…

[1]	Eines der spannenderen Beispiele für Information, die aus solchen Zeitverschiebungen kommt, ist die Bestimmung der Hubblekonstante mithilfe von gravitationsgelinsten Quasaren, ganz schön illustriert z.B. im Foliensatz 10.5281/zenodo.4062126.