von Thomas Viehmann tv@lernapparat.de

Eine der aktuell meistdiskutierten Kennzahlen der Epidemie ist die Reproduktionszahl R und deren Schätzung durch das Robert-Koch-Institut. In diesem Artikel vollziehen wir - so gut wie uns möglich - die zugrundeliegende Schätzung der Neuerkrankungszeitreihe. Das erlaubt es uns auch, analoge Kennzahlen für Bundesländer zu berechnen.

Das Kleingedruckte vorab: Ich bitte darum, meinen Codes oder aus diesem weiterentwickelte Programme oder der Berechnungsresultate von solchen nur unter Angabe dieses Artikels als Quelle zu verwenden. Ich stelle diese Informationen nur zu Informationszwecken zur Verfügung, jede Haftung für die richtigkeit der Ergebnisse schließe ich aus. Die erzeugten Daten unterliegen wie die zugrundeliegenden Daten des Robert-Koch-Instituts der Open Data Datenlizenz Deutschland – Namensnennung – Version 2.0 zur Verfügung.

Und noch eine Vorbemerkung: Wie man spätestens beim Nowcasting merkt, ist dies ein Arbeitsstand, der noch aktualisiert werden wird. Es schien mir für die Transparenzdiskussion nützlich, auch diesen Arbeitsstand zugänglich zu machen.

Neuinfektionen, Maßnahmen und Reproduktionszahl

Für die Erfassung der Entwicklung der Epidemie ist der Blick auf die Neuinfektionen maßgeblich. Die absolute Höhe ist wichtig - und anhand ihr, bezogen auf die Bevölkerungsgröße, wurde ja auch die "Notfall-Grenze" gezogen Wie es zu der Obergrenze kam, ist ja unklar. Das RKI hatte in der Vorwoche von ca. 400 Fällen pro Tag in Deutschland gesprochen, die sie für Verfolgbar halten. Das sind ca. 4-5 pro 1 Mio. Einwohner und Tag. Jetzt ist es 50 pro Woche und 100000 Einwohner, also ca. 70 pro 1 Mio. Einwohner und Tag und niemand konnte so recht sagen, wie man darauf kommt. Wenn das mal gut geht. Ein Schelm, wer denkt, es wäre ungefähr so: Ministerpräsident: Wann können wir endlich Lockerungen einführen? Experten: Das wäre möglich bei ca. 4-5 pro 1 Mio. Einwohner und Tag. Ministerpräsident: Oh, und wieviel haben wir jetzt? Experte: das Variiert von Landkreis zu Landkreis bis ca. 70. Ministerpräsident: Ich will aber jetzt Lockerungen, also darf die Obergrenze nicht niedriger als 70 sein..

Die Entwicklung der Fallzahlen, wird maßgeblich dadurch gesteuert, wieviele weitere Personen jeder Infizierte im Mittel ansteckt. Dies wird als die ReproduktionszahlDie größe heißt Zahl, weil es eine Anzahl ist, nicht, etwa, weil man "nur" zählt. Genauer wird die Reproduktionszahl unter der Annahme bestimmt, dass alle anderen Personen noch nicht immun sind. Wenn jetzt ein Anteil $\alpha$ der Bevölkerung schon immun ist (und die Bevölkerung homogen ist, was, wie z.B. in den Diskussionen über die Schlussfolgerungen aus der Gangelt-Studie immer wieder betont wurde, eher nicht der Fall ist), wäre die durchschnittliche Zahl der Ansteckungen $R * (1 - \alpha)$. So kommt ja auch die Herdenimmunitätsschwelle zustande als Schwelle $\alpha_{HI}$, so dass für die Basisreproduktionszahl $R_0$ gilt, dass die Ansteckungen unter eins liegt: Diese liegt daher bei $R_0 * (1 - \alpha_{HI}) = 1$. Wenn wir maßnahmen ergreifen, um $R$ von $R_0$, dass offenbar für SARS-CoV-2 bei 2-4 liegt, zu senken, haben wir einen Rückgang bei $R * (1 - \alpha)$. Da man aber davon ausgehen muss, dass $\alpha$ nahe 0 liegt (außer vielleicht in New York), nimmt man fast immer $\alpha=0$ an, wenn man jetzt über $R$ spricht und es misst, die Pedanten mögen sich überall $R / (1 - \alpha)$ denken. $R$ bezeichnet.

Wenn man man die Neuinfektionen auf in einem Diagramm mit logarithmischer vertikaler Achse aufträgt, gibt die Reproduktionszahl $R$ gerade eine (Tangenten-) Steigung anEine Steigung muss ja als Einheit die Einheiten der vertikalen Achse (Anzahl) geteilt durch die der horizontalen Achse (Zeit in Tagen) haben, daher kann es so einfach nicht sein. Genauer muss man die Reproduktionszahl durch das serielle Intervall, also die Zeit zwischen aufeinanderfolgenden Ansteckungen, teilen. Das serielle Intervall wird vom RKI auf ca. 4 Tage geschätzt.. Wir hatten in Den Rückgang der Epidemie messen die gemeldeten Fallzahlen genommen, geglättet, und die Steigung am Ende des Graphen extrapoliert.

Ein Nachteil unseres Verfahrens ist der Zeitbezug. Wenn das Verhalten, mithin die ergriffenen Maßnahmen, die Reproduktionszahl $R$ bestimmen (das muss man sich klar machen, dass $R$ durch das Verhalten gesetzt und nicht etwa mit einer Tendenz versehen wird), wollen wir natürlich möglichst zeitnah wissen, welches $R$ wir aktuell erreichen.

Um zu verstehen, wie es um das Aktuell dabei steht, muss man überlegen, was wir mit den Fallzahlen machen (im oben verlinkten Blogpost haben wir zweimal über 7 Tage gemittelt, das ergibt einen Verzug), aber auch, wie sie zustandekommen.

Wie die Fallzahlen zustandekommen

Werfen wir einen Blick auf das Meldewesen. Im Infektionsschutzgesetz ist die Kette der Meldungen bei meldepflichtige Krankheiten, zu denen CoViD-19 zählt, definiert. Das Robert-Koch-Institut empfielt eine Übersicht aus der Apotheken-Umschau. In der folgenden Abbildung sehen wir uns die Meldekette und Ihren Vorlauf von Infektion zu Diagnose an.

Am Anfang sind die Infektionen. Wie viele weitere Personen jeder Infizierte im Mittel ansteckt, ist die Reproduktionszahl $R$, die uns interessiert. Wir sehen sie jedoch nicht direkt. Eine infizierte Person durchläuft zunächst die Inkubationszeit (im Mittel laut RKI ca. 5 Tage), bevor sie Symptome entwickelt.Eines der Probleme ist gerade, dass die oft Symptome nicht oder erst nach Beginn der Infektiösität auftritt, sogar das serielle Intervall (mittlere Zeit, nicht Zeit bis zu erster Folgeinfektion) gilt als kürzer als die durchschnittliche Inkubationszeit. Es gibt auch Fälle, die entdeckt werden, bevor sie symptomatisch werden und solche, die entdeckt werden, obwohl sie asymptomatisch verlaufen. Dies kann zum Beispiel durch vorsorgliche Tests oder Kontaktpersonenverfolgung sein. Dann wird eine Probe vereinbart und genommen. Diese muss vom Labor ausgewertet werden. Dann beginnt - bei einem positiven Befund die Meldekette. Das lokale Gesundheitsamt wird informiert, dieses gibt die Daten an die zuständigen Ämter des Landes weiter. Von dort werden sie an das Robert-Koch-Institut als maßgebliche Bundesbehörde weitergeleitet. Das veröffentlicht dann die Situationsberichte und (mittelbar?) die berühmte Tabelle mit den Fallzahlen (RKI_COVID19.csv).

Wenn man jetzt von dort die Fallzahlen nimmt, wie im RKI-Dashboard angezeigt, und die Entwicklung verfolgt ist man am Ende der Meldekette, ganz unten. Das ist ganz schön weit weg von der eigentlich interessanten Ansteckung! Einige Zeitungen (z.B. die Zeit oder die Berliner Morgenpost) eliminieren den Meldeverzug, indem sie die Informationen direkt von den lokalen Gesundheitsämtern "einsammeln".

Das RKI geht in seiner Analyse noch ein stück Weiter und schaut für seine täglichen Situationsberichte auf das Erkrankungsdatum. Dazu verwendet es statistische Verfahren um von den ihm gemeldeten Fällen auf Erkrankungen einschließlich der ihm noch nicht gemeldeten Fälle zu schließen. Dieses Verfahren, genauer die Imputation und das Nowcasting wollen wir uns genauer ansehen. Der Weg vom Nowcasting zur Berechnung der Reproduktionszahl ist vom RKI sehr detailliert dokumentiert. Wie Björn Schwentker vom NDR beschreibt, ist aber der Weg von den gemeldeten Fällen zur Schätzung der Erkrankungszahlen nicht so leicht nachzuvollziehen.

Exkurs: Ist es realistisch zu wissen, was man noch nicht erfahren hat?

Man kann sich fragen - und es scheint Menschen zu geben, die sich da sehr empören können - ob es solche Schätzungen funktionieren können. Die Grundidee, von einem bekannten Teil auf ein größeres Ganzes zu schließen, ist ja nicht vollkommen abwägig, die Demoskopen haben sie für Umfragen perfektioniert.

Aber es gibt auch andere Gelegenheiten, bei denen noch ähnlichere Verfahren weitgehend geräuschlos verwendet werden. Wenn wir zum Beispiel den Geschäftsbericht 2019 der Münchener Rück ansehen, finden wir auf Seite 101 die Verbindlichkeiten in der Bilanz wie in der nächsten Abbildung gezeigt.

Die Rückstellung für noch nicht abgewickelte Versicherungsfälle, also Schadensfälle, die schon eingetreten sind, aber bei denen die Versicherungsleistung noch nicht gezahlt wurde, beträgt 70,9 Milliarden Euro - fast ein Viertel der gesamten Bilanzsumme. Auf Seite 151 ist diese dann aufgegliedert, wie wir in in folgender Abbildung.

In der Schadenrückstellung entfallen 36,9 Mrd. Euro auf die Spätschadenrückstellung. Auf Englisch heißt sie incurred but not reported, kurz IBNR, d.h. eingetreten aber noch nicht gemeldet. Diese Schäden werden - wie beim Nowcasting die noch nicht gemeldeten Infektionsfälle - statistisch geschätzt. Gut ein Achtel der Verbindlichkeiten in der Bilanz der Münchener Rückversicherung beruht also auf einer Art NowcastingEs ist aber nicht alles "Geld der Münchener Rück", da ein Teil (2,9 Mrd EUR) im Rahmen der Retrozession bei anderen Rückversicherern weiterversichert wurde..

Rekonstruktion der Imputation und des Nowcasting des RKI

Schauen wir also genauer darauf, wie das Robert-Koch-Institut auf Erkrankungszahlen kommt, die seiner Schätzung der Reproduktionszahl zugrunde liegen. Primärquelle ist der Artikel Schätzung der aktuellen Entwicklung der SARS-CoV-2-Epidemie in Deutschland - Nowcasting von an der Heiden und Hamouda im Epidemiologischen Bulletin 17/2020. Die Autoren machen einen Kompromiss - der Artikel ist für eine interessierte Öffentlichkeit verständlich gehalten, dafür sind manche Details, die eine Rekonstruktion vereinfachten, weggelassen. Die andere Text-Quelle sind die Beschreibungen in den täglichen Situationsberichten.

Nachdem ich - im Hinblick auf die NDR-Berichterstattung und meinen Stand bei der Rekonstruktion - geäußert hatte, dass ich die Methode für im Wesentlichen nachvollziehbar halte, hat sich ein interessanter Austausch mit Björn Schwendtker und auch Peter Steinbach ergeben. Vielen Dank!

Im Sinne der Aufklärung beschreibe ich hier also mein Verständnis des Verfahrens und den Code, den ich dazu geschrieben habe, die Methode nachzubauen und eigene Nowcasts, zum Beispiel für Bundesländer, zu machen.

Ich sollte an dieser Stelle sagen, dass ich weder Epidemiologe noch Journalist binOscar Wilde hat ja "Social Media" beschrieben, auch wenn er über Bücher spricht: In old days books were written by men of letters and read by the public. Nowadays books are written by the public and read by nobody., sondern Mathematiker. Wenn ich für meine Kunden mathematische Modellierung für Machine Learning betreibe, sind aber auch die Daten das A und O. Vor einiger Zeit habe ich auch mal für eine Europäische Aufsichtsbehörde die Übereinstimmung der Berechnungen mit der Dokumentation nachvollzogen. Insofern bin ich mir bewusst, wie schwierig das selbst mit Code ist.

Überblick

Die Methode des RKI besteht im Wesentlichen aus zwei Teilen:

• Zunächst werden fehlende Angaben in den bekannten Fällen ergänzt. Konkret liegt für ungefähr zwei Drittel der Fälle ein Erkrankungsdatum vor, für das übrige Drittel nicht. Das kann verschiedene Ursachen habe, z.B. asymptomatische Fälle (d.h. ohne Erkrankung), vor der entwicklung Symptomen festgestellte Erkrankungen (obwohl die auch nachgemeldet werden), oder ein unbekanntes Erkrankungsdatum. Da wir ja aber die Erkrankungen letztlich als Hilfsgröße für die Schätzung der Infektionszahlen benutzen wollen, tun wir so, als würde sei ein vorhandenes Erkrankungsdatum nur nicht gemeldet worden. Wir füllen diesen fehlenden Wert aufgrund der unserer Beobachtungen für die bekannten Erkrankungsdaten. Dieses Vorgehen ist die Imputation. Hier geht es also um die Vervollkommnung der Aufzeichnungen für die schon gemeldeten Fälle.

• Um den Meldeverzug zu behandeln, schätzen wir - wie die Münchener Rück für Ihre Bilanz - was noch auf uns zu kommt. Für Erkrankungsdaten, für die schon einige Fälle vorliegen, kann man aus der Beobachtung des Meldeverzugs in der Vergangenheit schließen, wie viele Fälle uns nach unserem Datenstand noch erreichen werden. Dies ist das Nowcasting.

Damit erhalten wir eine Zeitreihe der Neuerkrankungszahlen. Aus dieser wird dann ein gleitendes 4-Tages- (und neuerdings 7-Tages-) Mittel gebildet. Dann wird schließlich die Reproduktionszahl $R$ als Quotient von vier Tage (dem seriellen Intervall) auseinanderliegenden Neuerkrankungswerten geschätzt.

Eingangsdaten

Wir haben natürlich nur öffentliche Daten, das heißt, wir verwenden die oben Erwähnten RKI_COVID19.csv. Die Datei kommt mit einer ausführlichen Datensatzbeschreibung. Das RKI hat aber mehr Informationen: Zum einen hat das RKI Informationen auf Einzelfallebene, den es vermutlich auch über die Zeit verfolgen kann. Diesen "Luxus" haben wir nicht. Inwiefern es Datenschutzrechtlich möglich wäre, die Einzelfälle (vermutlich wäre das pseudonymisiert) herauszugeben, weiß ich nicht. In jedem Fall fehlt uns das "Eingangsdatum" beim RKI (in der Abbildung oben am Ende des Meldeverzugs).

Was wir aber machen können - und nach meinen Erkenntnissen auch müssen, um die Daten anzunähern - ist, verschiedene Datenstände in der Gesamtsicht zu sehen. Das ist im Vergleich zur Behandlung von Einzelfällen über die Zeit ungenauer. Es ist gut möglich, dass man mit aktuell öffentlich verfügbaren Daten keine vollständige Übereinstimmung mit den RKI-Werten erreichen kann. Warum, sehen wir, wenn wir die Imputation im Detail angucken.

Vergleichsdaten

Es ist versuchend, einfach zu sehen, ob man das $R$ des RKI reproduzieren kann. Das ist ja aber das Ergebnis eines mehrstufigen, relativ komplexen Prozesses:

Ausgehend auf den gemeldeten Daten,

• werden die fehlenden Erkrankungsdaten imputiert,
• dann wird das Nowcasting gemacht,
• es werden gleitende Mittelwerte der Ergebnisse berechnet,
• daraus wird $R$ geschätzt.

Das jetzt in einem Rutsch alles nicht nur "richtig", sondern auch noch genauso wie das RKI zu machen, ist relativ unwahrscheinlich. Deshalb sind Zwischenergebnisse für uns sehr wichtig. Das RKI stellt die Nowcasting-Ergebnisse als Zahlen zur Verfügung.

Aber damit sind die beiden schwierigen Schritte - die Imputation und das Nowcasting - immer noch zusammen. Wie können wir die trennen?

Als Approximation können wir versuchen, die Fallzahlen aus dem Graphen im Situationsbericht des RKI zu benutzen. Dann bekommen wir die Zahlen zu folgender Abbildung:

Das ist per Hand mühsam. Automatisiert auch, ich habe ein Auswertungsskript in Python geschrieben, dass eine SVG-Datei ausliest.Vermutlich ist das aber nicht sehr "stabil", man muss also damit rechnen, dass man das Programm ergänzen muss, wenn im PDF oder der Übersetzung in SVG etwas komisches passiert. Direkt mit dem PDF zu arbeiten, wäre natürlich einfacher in der Bedienung, aber offene Dateiformate sind so viel schöner... Mit einer UI kann man da bestimmt ein ganzes Data-Startup draus machen. Äh. Vermutlich hat das schon jemand gemacht..

Unser Nowcasting

Also nun an die Arbeit. Wir wollen die Methode in Python nachverfolgen.

Vorbereitung: Daten laden

Als erstes müssen wir die Daten laden. Ich benutze als Ausgangsmaterial eine Sammlung alter RKI_COVID19.csv. Weil der Datenstand ein sehr deutsches Format hat, übersetze ich den per Hand.

Zuerst machen wir einen großen Datenframe aus den verschiedenen Datenstände. Relevant für uns sind jene mit IstErkrankungsbeginn und Refdatum.

dfs = []
for date in pandas.date_range('2020-04-29', '2020-05-17', freq='D'):
    dstr = date.date().isoformat()
    adf = pandas.read_csv(f'data/rki_covid19.{dstr}.csv')
    adf['Datenstand'] = pandas.to_datetime(adf.Datenstand.apply(lambda s: '-'.join(s.split(',')[0].split('.')[::-1])))
    assert (date==adf['Datenstand'].min()) and (date==adf['Datenstand'].max())
    dfs.append(adf)
large_df = pandas.concat(dfs, axis=0, sort=True)

Wir machen die Datumsspalten zu echten Daten.

for k in ('Meldedatum', 'Datenstand', 'Refdatum'):
    large_df[k] = pandas.to_datetime(large_df[k])

Schließlich führen wir noch drei Hilfsspalten ein: - Die Zeit zwischen Erkrankung und Meldedatum (wenn der Erkrankungsbeginn angegeben ist). - Eingangsdatum und Verzögerung bis zum Eingang. Diese sind nur für als NeuerFall markierte Daten gültig. Wir benötigen das Eingangsdatum nachher summarisch auch für andere Fälle, das kommt aber später.

large_df['delay'] = (large_df.Meldedatum - large_df.Refdatum).dt.days # nur gültig für IstErkrankungsbeginn==1
large_df['Eingangsdatum'] = large_df.Datenstand - pandas.Timedelta(days=1) # nur gültig für NeuerFall==1
large_df['delay_eingang'] = (large_df.Eingangsdatum - large_df.Refdatum).dt.days # nur gültig für NeuerFall==1

Damit haben wir unsere Daten beisammen und wir können anfangen, die drei Fall-Zeitreihen des RKI zu rekonstruieren.

Fälle mit angegebenem Erkrankungsbeginn

Im vom RKI veröffentlichten Datensatz gibt es (seit 29. April) Spalten IstErkrankungsbeginn und Refdatum. Ist in ersterer eine 1 hinterlegt, gibt letztere das angegebene Datum des Erkrankungsbeginns an.

Soweit so gut. An der Heiden und HamoudaAlle unmarkierten Zitate im folgenden sind aus an der Heiden M, Hamouda O: Schätzung der aktuellen Entwicklung der SARS-CoV-2-Epidemie in Deutschland – Nowcasting. Epid Bull 2020;17:10 – 16; DOI 10.25646/6692.4 schreiben

In 522 Fällen war der zeitliche Abstand zwischen dem Datum der Übermittlung an das RKI und dem Erkrankungsbeginn negativ oder lag über 30 Tage. Diese Fälle wurden bei der nachfolgen- den Analyse und der Imputation des Erkrankungsbeginns nicht miteinbezogen.

Wir müssen also um Fälle bereinigen, deren Erkrankungsbeginn vor dem Übermittlungsdatum (Eingangsdatum am RKI) liegt. Für den aktuellen Tag ist das möglich, da die Spalte NeuerFall mit 1 anzeigt, dass der Fall in dem Tag bis zum Datenstand (der ja 0 Uhr des Folgetages ist) eingegangen ist. (Ein negativer Wert schint da ohne Hellsehen nicht möglich...) Für zurückliegende Tage ist das nicht so leicht möglich: Zwar könnten wir durch rekurrieren auf ältere Datenstände den damaligen Stand abfragen (ab dem 29.4.), jedoch werden ja auch einige Datensätze wieder "ausgelistet", vermutlich wegen Korrekturen, so dass wir nicht wissen, ob die Fälle im neuesten Datenstand noch auftauchen.

Hier scheint es aber kein so großes Problem zu sein, so dass wir für unsere Zahlen von Fällen mit angegebem (oder auch bekanntem) Erkrankungsbeginn die Verzögerung zum Meldedatum als Hilfsmittel nehmenMan könnte natürlich auch eine Approximation durchführen, wie wir sie gleich für die Fälle mit unbekanntem Erkrankungsbeginn beschreiben..

Mit dieser Vereinfachung haben wir folgende Methode, um in Python die Zeitreihe der Fälle mit angegebenem Erkrankungsbeginn zu bekommen. Wir benutzen nur Daten aus der neuesten Datenstands-Meldung und Filtern auch die "gelöschten Fälle" (mit NeuerFall als -1 markiert) weg.

def faelle_mit_erkrankungsbeginn_nach_erkrankungsbeginn(large_df):
    # eigentlich müsste es Delay zum jeweiligen Eingang sein. Kennen wir i.A. aber nicht
    known = large_df[(large_df.Datenstand == large_df.Datenstand.max())
                          & (large_df.IstErkrankungsbeginn == 1) 
                          & (large_df.NeuerFall >= 0)
                          & (large_df.delay >= 0)
                          & (large_df.delay <= 30)].groupby(['Refdatum']).AnzahlFall.sum()
    known = known.reindex(pandas.date_range('2020-01-01', large_df.Eingangsdatum.max()))
    return known

Um mit den aus den Graphen der Situationsberichte abgelesenen Daten zu vergleichen, müssen wir noch die gleitenden 4-Tage-Mittel benutzen. Da uns das öfter begegnet, definieren wir eine Hilfsfunktion und stellen die Mittelung als Faltung mit einem konstanten Glättungskern dar.Das haben Sie davon, wenn es ein Mathematiker aufschreibt.

def smooth(s, size=4):
    return pandas.Series(numpy.convolve(s, numpy.full((size,), 1/size))[:-(size-1)], index=s.index, name=s.name)

Wir bekommen eine Zeitreihe, die sehr nah an jener des RKI in seiner Grafik ist, wie die folgende Abbildung zeigt.

Damit haben wir die Fälle mit angegebenem Erkrankungsbeginn behandelt.

Imputation

Die Schwierigkeit des fehlenden Eingangsdatums begegnet uns für die Imputation wieder, aber in einer verschärften Form.

Als wichtigste Angabe zur Ermittlung des fehlenden Erkrankungsbeginns wurde dabei das Datum des Eingangs der Fallmeldung am RKI verwendet, die fehlenden Werte wurden getrennt nach Geschlecht und Altersgruppe geschätzt.

Wenn man die Datei so vor sich hat, ist man erst mal versucht, das Meldedatums statt des Übermittlungsdatums zu benutzen. Die folgende Grafik zeigt die Fallzahlen für beide Varianten.

Beide Varianten - Imputation nach Meldedatum oder Eingangsdatum - sind legitim, führen aber naturgemäß zu voneinander abweichenden Ergebnissen. Wenn wir aber wie das RKI das Eingangsdatum benutzen wollen, benötigen wir einen Trick, um die Eingangsdaten der Fälle in der aktuellen Fall-Tabelle zu approximieren.

Übermittlungsdaten approximieren

Wie also das Übermittlungsdatum approximieren? Für die neuesten Fälle mit $NeuerFall=1$ haben wir die Information. Wir können aber wegen der Korrekturen ($NeuerFall = -1$ am Ende des Korrekturtages) nicht einfach davon ausgehen, dass die Fälle, die in alten Datenständen so markiert sind, auch tatsächlich noch ohne Erkrankungsdatum sind.

Um dennoch approximativ die Übermittlungsdaten zu schätzen, müssen wir um diese Korrekturen bereinigen. Wir haben also folgende Situation: Im Datenstand zu Tag $d+1$ werden (in einer der Rasterkategorien, z.B. Alter und GeschlechtNach der Beschreibung zu schließen scheint das RKI eventuell andere Altersgruppen (20-Jahres-Bänder) zu benutzen als im öffentlichen Datensatz. Man könnte das Raster für diese Suche nach Übermittlungsdaten auch noch durch Hinzunahme der Landkreise verfeinern. Das gibt ggf. genauere Ergebnisse., $N_{d+1, \textrm{alt}}$ "alte Fälle" ($NeuerFall = 0$ im $d+1$-Datensatz) ausgewiesen. Am Tag $d$ werden $N_{d,\textrm{neu}} + N_{d,\textrm{alt}} \geq N_{d+1,\textrm{alt}}$ ausgewiesen. Wir müssen uns jetzt entscheiden, wie wir $N_{d,\textrm{neu}}$ und $N_{d,\textrm{alt}}$ anpassen, damit die Summe genau $N_{d+1,\textrm{alt}}$ wird. Ich habe mich dazu entschlossen, in erster Linie die alten Fälle des Tages $d$ zu vermindern, um die Differenz auszugleichen.Es gibt andere Möglichkeiten, diese Differenz zu verarbeiten. Naheliegend wäre zum Beispiel eine ratierliche Kürzung. Mit dieser habe ich jedoch zum Datenstand 17. Mai weniger gute Resultate bei der Imputation beobachtet. Eigentlich sollte ich - solange das RKI keinen entsprechenden Datensatz rausrückt - die Übermittlungsdaten so rekonstruieren und eine angepasste RKI_COVID19.csv rausrücken. Dann könnte man nachfolgend auch weitgehend ohne Historie arbeiten. Wir setzen also für den letzten Tag (pro Rasterkategorie) $\hat N_{d_{max}, \textrm{alt/neu}} = N_{d_{max}, \textrm{alt/neu}}$ und dann rückwärts über die Datenstände laufend $\hat N_{d, \textrm{alt}} = \max(0, \hat N_{d+1, \textrm{alt}} -N_{d, \textrm{neu}})$ sowie $\hat N_{d, \textrm{neu}} = \hat N_{d + 1, \textrm{alt}} - \hat N_{d, \textrm{alt}}$.

Für die Datenstände vor dem 29. April setzen wir einfach das Meldedatum als Eingangsdatum. Das führt zu erheblichen ungenauigkeiten, jedoch ist dieser Zeitraum für aktuelle R-Schätzungen weniger interessant. Man könnte zum Beispiel alle Fälle der vorherigen Daten wie solche ohne Erkrankungsdatum behandeln. Damit sind bestimmt noch Verbesserungen zu erzielen.

Wir machen das nur für die Fälle mit unbekanntem Erkrankungsdatum. Ich habe mich entschlossen, als Datum das Übermittlungsdatum (also Datenstand minus 1 Tag) zu verwenden. Im Code sieht das dann so aus:

def faelle_ohne_erkrankungsbeginn_nach_rekonstruiertem_eingangsdatum(large_df, missing_tol=2):
    unbekannte_roh = large_df[(large_df.IstErkrankungsbeginn == 0)
                              & (large_df.NeuerFall >= 0)].groupby(['Datenstand', 'NeuerFall',  'Meldedatum', 'Altersgruppe', 'Geschlecht']).AnzahlFall.sum()
    unbekannte_roh = unbekannte_roh.reindex(pandas.MultiIndex.from_product(unbekannte_roh.index.levels, names=unbekannte_roh.index.names),
                                            fill_value=0)
    new_cases_melde = {}
    old_cases = None
    newer_sum = 0
    for d in unbekannte_roh.index.levels[0][::-1]:
        if old_cases is None:
            old_cases = unbekannte_roh[d, 0]
            new_cases = unbekannte_roh[d, 1]
        else:
            last_old_cases = old_cases

            old_cases = (old_cases - unbekannte_roh.loc[d, 1]).clip(lower=0)
            new_cases = last_old_cases - old_cases
        newer_sum += new_cases.sum()
        # Eingangsdatum = Datum vor Datenstand
        new_cases_melde[d - pandas.Timedelta(days=1)] = new_cases

    new_cases = pandas.DataFrame(new_cases_melde).sort_index(axis=1)
    new_new = new_cases.sum(level=[1,2])
    unknown_new_cases_per_day = pandas.concat([old_cases[:d-pandas.Timedelta(days=2)].unstack(level=0, fill_value=0), new_new], axis=1)

    missing = old_cases[d-pandas.Timedelta(days=1):]
    missing = missing[missing > 0]
    # Wir verlieren für D nach Alter/Geschlecht mind. 2 Fälle, die offenbar nie als neuer Fall markiert sind.
    # Für Teil-Daten ggf. mehr.
    assert missing.sum() <= missing_tol
    return unknown_new_cases_per_day

Jetzt wollen wir aber endlich imputieren.

Vom Übermittlungsdatum zur Imputation

Mit dieser Vorbereitung ist die Imputation ganz einfach. Für die Historie benutzen wir (ungenau) die neuen Fälle mit zu den jeweiligen Datenständen und Gruppieren nach unserer Verzögerungs-Spalte (und dem Kategorienraster). Wieder nehmen wir für die Daten vor dem 29. April einfach das Meldedatum.

Im Code:

def verzug_von_erkrankung_bis_eingang(large_df):
    bekannte_roh = large_df[(large_df.IstErkrankungsbeginn == 1) 
                             & (large_df.NeuerFall == 1)
                             & (large_df.delay_eingang >= 0)
                             & (large_df.delay_eingang <= 30)].groupby(['Eingangsdatum',  'Altersgruppe', 'Geschlecht', 'delay_eingang']).AnzahlFall.sum()
    bekannte_roh = bekannte_roh.unstack(level=0, fill_value=0)
    bekannte_roh = bekannte_roh.reindex(pandas.MultiIndex.from_product(bekannte_roh.index.levels[:-1]+[pandas.RangeIndex(0, 31)], names=bekannte_roh.index.names),
                                            fill_value=0)
    bekannte_roh_old = large_df[(large_df.IstErkrankungsbeginn == 1) 
                             & (large_df.Datenstand == large_df.Datenstand.min())
                             & (large_df.NeuerFall == 0)
                             & (large_df.delay >= 0)
                             & (large_df.delay <= 30)].groupby(['Meldedatum', 'Altersgruppe', 'Geschlecht', 'delay']).AnzahlFall.sum()
    bekannte_roh_old = bekannte_roh_old.reindex(pandas.MultiIndex.from_product(bekannte_roh_old.index.levels, names=bekannte_roh_old.index.names),
                                            fill_value=0).unstack(level=0)
    # wir verlieren einen Fall
    assert bekannte_roh_old.loc[:, bekannte_roh.columns.min():].sum().sum() <= 1
    delay_known_new_cases = pandas.concat([bekannte_roh_old.loc[:, :bekannte_roh.columns.min()-pandas.Timedelta(days=1)], bekannte_roh], axis=1)
    delay_known_new_cases = delay_known_new_cases.reindex(pandas.date_range('2020-01-01',
                                                     delay_known_new_cases.columns.max()), axis=1, fill_value=0)
    delay_known_new_cases.index.names = ['Altersgruppe', 'Geschlecht', 'delay']
    delay_known_new_cases.columns.name = 'Eingangsdatum'
    return delay_known_new_cases

Wenn wir Verteilungen für die Zeitbereiche mit detaillierten Daten und ohne ansehen, erhalten wir die Verteilungen der jeweiligen Abbildungen.

Wir sehen insbesondere, dass die Verzögerungen nach Meldetag klarerweise kleiner sind als die echt nach Übermittlungsdatum. Das ist ja aber nicht völlig unerwartet.Fun fact: In dem oben verlinkten Apotheken-Umschau-Bild wird 2021 als Zieltermin für die digitalisierung des Meldewesens genannt.

Wie im Epidemiologischen Bulletin beschrieben, wird nach Altersgruppe und Geschlecht getrennt imputiert, also auch die Verteilungen getrennt ermittelt. Das ist relativ detailliert - in manchen Kategorien gibt es dann an einem Tag nur wenige Fälle um die Verteilung zu schätzen (zum Glück meist auch in denen mit wenig zu schätzenden Fällen). Weitere Dimension, mit, denen man Arbeiten könnte wären Bundesländer (falls der Verzug zwischen Meldung und Übermittlung unterschiedlich ist) oder die Zeit (man könnte zum Beispiel die Verzugsverteilung anhand eines Zeitfensters schätzen). Mein Eindruck aus dem Vergleich meiner Ergebnisse mit denen des RKI ist, dass das RKI nur die neuen Fälle des jeweiligen Tages für Schätzung der die Imputations-Verteilungen heranzieht. Bundesländer betrachtet das RKI nicht getrennt. Wenn man nach Bundesländern trennt, ist es gegebenenfalls ratsam, die andere Segmentierung zu vergröbern oder zum Schätzen der Verteilung auf den gesamten Datensatz (also ganz Deutschland) zu gehen, um eine vernünftige Grundlage für die Verzugsverteilung zu haben. Ich habe mich für die Bundesländer-Nowcasts unten für letzteres entschieden.

Mit der Verteilung und den zu imputierenden Fällen haben wir haben nun verschiedene Möglichkeiten, die Imputation technisch zu bewerkstelligen. Eine ist, wie das RKI stochastische Pfade zu produzieren, also für jeden Fall $K=200$ Ziehungen aus der Verzugsverteilung vorzunehmen und das resultierende Erkrankungsdatum zu berechnen. Aggregiert man jeweils die $i$-te Ziehung zu einer Zeitreihe, bekommt man $K$ mögliche Zeitreihen von Fällen pro Erkrankungsdatum. Der Mittelwert gibt dann den Punktschätzer und man kann quantile für die Vorhersageintervalle benutzen.

Da wir erst mal vor allem am Mittelwert interessiert sind, gehen wir einen leicht anderen Weg: Wir teilen den Fall entlang der Verzugsverteilung auf und aggregieren so zur Zeitreihe. Wir schreiben mit, wo wir keine Fälle für die Schätzung der Imputationsverteilung im Raster hatten. Meist sind dies Kategorien mit unbekanntem Geschlecht o.ä. mit sehr wenigen zu imputierenden Fällen, diese werden dann implizit auf die letzte Kategorie mit Fällen abgebildet (also z.B. weibliches Geschlecht).

Der Code hat als Kern eine dreifache Schleife über Datum, Altersgruppe und Geschlecht. Da ich zuvor mehrere Tage gesamthaft betrachtet habe, benutze ich wieder eine Faltung, das könnte man natürlich auch anders machen, wenn man Tagesweise vorgeht. Aus dem gleichen Grund bestimme ich hier die Verteilung aus den Fallzahlen "on the fly".

def imputiere_erkrankungsdaten(unknown_new_cases_per_day, delay_known_new_cases, reuse_tol=200, completely_missing_tol=0):
    dist = None
    melde_range = pandas.date_range('2020-01-01', unknown_new_cases_per_day.columns.max(), name='Refdatum')
    reuse_count = 0
    completely_missing_count = 0
    imputed_by_day_and_cat = {}
    for date in unknown_new_cases_per_day.columns:
        ##['A00-A04', 'A05-A14', 'A15-A34', 'A35-A59', 'A60-A79', 'A80+', 'unbekannt']:
        for ag in unknown_new_cases_per_day.index.levels[0]:
            for g in unknown_new_cases_per_day.index.levels[1]:
                section = unknown_new_cases_per_day.loc[ag, g]
                section = section[date:date] # note: This is including the upper boundary!
                section = section.reindex(melde_range, fill_value=0)
                count = section.sum()
                unnormed_dist = delay_known_new_cases.loc[ag, g].loc[:, date]
                unnormed_count = unnormed_dist.sum()
                if unnormed_count == 0 and count > 0:
                    # print("### Warning", date, ag, g, count)
                    reuse_count += count
                elif unnormed_count > 0:
                    dist = unnormed_dist / unnormed_count
                if count > 0 and dist is not None:
                    imputed_section = numpy.convolve(section, dist[::-1])[-len(section):] # we only want temporal padding
                    imputed_by_day_and_cat[(date, ag, g)] = pandas.Series(imputed_section, index=melde_range)
                    assert count - imputed_section.sum() < 2, "loosing too many cases in imputation"
                elif count > 0:
                    completely_missing_count += count
    if INTERACTIVE:
        print("reused:", reuse_count, "completely missing:", completely_missing_count)
    assert completely_missing_count <= completely_missing_tol # in D eigentlich nicht.
    assert reuse_count <= reuse_tol  # Stand 17. Mai: 179 Einzelfälle ohne andere Fälle in gleicher Kategorie in D
    imputed =  pandas.DataFrame(imputed_by_day_and_cat).T
    imputed.index.names = ['Eingangsdatum', 'Altersgruppe', 'Geschlecht']
    #imputed.columns.name = 'Eingangsdatum'
    imputed = imputed.reindex(pandas.date_range('2020-01-01', large_df.Eingangsdatum.max(), name=imputed.columns.names[0]), axis=1)
    return imputed

Die nächste Abbildung vergleicht die so gewonnene Zeitreihe imputierter Erkrankungsbeginne (ohne Bekannte Fälle) mit den Zahlen des RKI. Wir sehen, dass die Übereinstimmung im späteren Bereich mit ca. \pm 10 Fällen sehr gut ist.

Damit haben wir die Imputation, nun kommen wir zum eigentlichen Nowcasting.

Nowcasting

Für das Nowcasting können wir wieder die dieselbe Verzögerungsverteilung benutzen, wir aggregieren aber alle (neuen) Fälle eines Tages mit bekanntem Erkrankungsbeginn.

Wir werden hier eine einfachere Methode, als sie das RKI (mutmaßlich) benutzt, vorstellen, das Verfahren von Lawless, Adjustments for Reporting Delays and the Prediction of Occurred but Not Reported Events, Can J Stat 20 (1994). Wieder wollen wir die Verzögerungsverteilung schätzen, aber diesmal aus den 8 letzten Zeitschritten (Datenstände $d_{max}$ bis $d_{max}-7$).

Die übliche Darstellung erfolgt Anhand von Erkrankungsbeginn und VerzögerungWenn Sie die Dreiecke in den Arbeiten on Höhle und an der Heiden und Lawless verwirrt haben, sind sie nicht allein. Selbst im R-Paket surveillance, in dem das Nowcasting von M. Höhle programmiert wurde, findet sich an einer Stelle in der Implementierung des Lawless-Verfahrens eine Index-Berechnung mit einem Verweis audf eine Formel in der Arbeit von Lawless und "knifflig" im Kommentar: x = T-t + seq_len(min(T1-t,D)-(T-t)) #(2.13) with modifications. knifflig. Ich selbst habe es auf Anhieb auch prompt nicht richtig gemacht.. Dann erhält man ein Dreieck, aus denen man bestimmte Felder zu rate zieht. Das machen wir jetzt auchIch habe in einer früheren Version an dieser Stelle versucht, einen einfacheren Weg zu gehen, dieser führt dann aber, worauf M. Höhle mich hinwies, nicht zur Lawless-Schätzung, so dass ich das umgestellt habe. Vielen Dank!.

Für das Dreieck müssen wir die Verzugsdaten nach Erkrankungsdatum statt wie oben nach Übermittlungsdatum betrachten. Dazu wandeln wir den Code ab:

def verzug_von_erkrankung_bis_eingang_nach_erkrankung(large_df):
    bekannte_roh = large_df[(large_df.IstErkrankungsbeginn == 1) 
                             & (large_df.NeuerFall == 1)
                             & (large_df.delay_eingang >= 0)
                             & (large_df.delay_eingang <= 30)].groupby(['Refdatum',  'Altersgruppe', 'Geschlecht', 'delay_eingang']).AnzahlFall.sum()
    bekannte_roh_old = large_df[(large_df.IstErkrankungsbeginn == 1) 
                             & (large_df.Datenstand == large_df.Datenstand.min())
                             & (large_df.NeuerFall == 0)
                             & (large_df.delay >= 0)
                             & (large_df.delay <= 30)].groupby(['Refdatum', 'Altersgruppe', 'Geschlecht', 'delay']).AnzahlFall.sum()

    bekannte_roh.index.names = bekannte_roh.index.names[:-1] + ['delay']
    bekannte_roh = bekannte_roh.unstack(level=0, fill_value=0)
    bekannte_roh_old = bekannte_roh_old.unstack(level=0, fill_value=0)

    bekannte_roh = bekannte_roh.reindex(pandas.MultiIndex.from_product(bekannte_roh.index.levels[:-1]+[pandas.RangeIndex(0, 31)], 
                                                                       names=bekannte_roh.index.names),
                                            fill_value=0)
    bekannte_roh_old = bekannte_roh_old.reindex(pandas.MultiIndex.from_product(bekannte_roh_old.index.levels[:-1]+[pandas.RangeIndex(0, 31)], 
                                                                       names=bekannte_roh_old.index.names),
                                            fill_value=0)

    return bekannte_roh.add(bekannte_roh_old, fill_value=0)

Tragen wir also in den Zeilen die Erkrankungsbeginne ab und in den Spalten den Verzug bis zum Eingang. Zum Datenstand 17. Mai (bzw. jeweils als neuen markierten Fällen der Datenstände bis zum 17. Mai) erhalten wir die folgende Tabelle. Jede Zeile enthält Eingangszahlen zu einem Erkrankungsdatum, in den Spalten die Dauer bis zum Eingang beim Robert-Koch-Instituts (ohne nachträgliche Korrekturen, da wir jeweils den Datenstand nehmen). In der Zeile zu Erkrankungsbeginn 12. Mai gingen also die Fälle aus der Spalte für 0 Tage Verzögerung am selben Tag beim RKI ein, aus der Spalte 1 Tag am 13., aus der Spalte 2 am 14. und so weiter. Klarerweise ist am 16. Schluss, daher die Dreiecksgestalt.

Das RKI schätzt nun

die Wahrscheinlichkeit dass ein Fall, der mit einer Verzögerung, $D$, von höchstens $d$ Tagen übermittelt wird, genau nach $d$ Tagen übermittelt wird, $g(d) = P(D = d | D ≤ d)$Hier ist $D$ die Zufallsvariable, die die Verzögerung beschreibt und $d$ eine Grenze, die wir festlegen.. Um $g$ für einen Verzug von $d$ Tagen zu schätzen, verwenden wir alle Fälle mit Erkrankungsbeginn im Bereich $T-d$ und $T-d-7$.

Um diese 8-Tages-Bereiche abzudecken, Tragen wir für jeden letzten Erkrankungsbeginn die Summe aus 8 Werten nach unten summiert ein. Die 8 rot markierten Werte oben ergeben zum Beispiel den rot markierten Wert in der folgenden Tabelle.

Von Interesse sind jetzt die jeweils letzten eingetragenen Werte in einer Zeile und die Summe der Zeile, bei Lawless heißen sie jeweils $N_x$ und $n_x$, wobei $x$ die Zeile bzw. die Verzögerung angibt, die bei uns $d$ heißt. Der letzte Wert $n_d$ ist die Anzahl der Fälle, die in den letzten 8 Tagen mit einer Verzögerung von $d$ Tagen eingegangen sind. $N_d$ ist die Gesamtzahl der Fälle zu den zu $n_d$ gehörenden Erkrankungstagen, die mit einer Verzögerung von maximal $d$ Tagen eingegangen sind. Damit ist der bei Lawless analysierte Quotient $$\hat g(d) = \frac{n_d}{N_d}$$ ein Schätzer für die obige bedingte Wahrscheinlichkeit $g(d)$.

Nun ist die Wahrscheinlichkeit, dass ein Fall, der mit höchstens $d$ Tagen Verzögerung übermittelt wird, mit genau $d$ Tagen Verzögerung übermittelt wird, noch nicht so interessant. Wir können aber einen Schritt machen, indem wir das komplementäre Ereignis betrachten, also die Wahrscheinlichkeit $$P(D \leq d -1| D \leq d) = 1 - g(d)$$ berechnen, dass ein Fall, der nach höchstens $d$ Tagen übermittelt wird, schon nach $d-1$ Tagen ankommt. So kommen wir von $d$ zu $d-1$!

Wenn wir jetzt annehmen, dass alle Fälle nach maximal $30$ Tagen übermittelt werden, können wir uns mit dieser Formel sukkzessive nach vorn hangeln, denn (weil die Bedingung in den bedingten Wahrscheinlichkeiten immer notwendig für das Ereignis ist) ist $$ F(x) = P(D \leq x) = P(D \leq x | D \leq 30) = \prod_{d = x+1}^{30} P(D \leq d-1 | D \leq d) = \prod_{d=x+1} (1 - g(d)). $$

Wenn wir für die Wahrscheinlichkeiten $g(d)$ jeweils die Schätzer $\hat g(d)$ einsetzen, erhalten wir einen Schätzer $\hat F(x)$ dafür, dass ein Fall mit $x$ Tage nach Erkrankungsbeginn gemeldet wird.

Beispielsweise schätzen wir also, dass $38.7\%$ aller Fälle mit maximal 3 Tagen Verzögerung gemeldet werden. Wenn wir also nach 3 Tagen für den 13. Mai nach Imputation(!) $264$ Fälle übermittelt wurden, schätzen wir, dass insgesamt $264 / 38.7\% \approx 682$ Fälle aufgetreten sind und im Laufe der Zeit gemeldet werden.

Diese Zahl liegt rund ein Zehntel höher als die vom RKI geschätzte Zahl von $620$. Dies lässt sich vor allem dadurch erklären, dass das RKI nach Imputation nur ca. 240 Fälle mit Erkrankungsbeginn 13.5. hatWir haben ja keine ungeglätteten Imputationswerte des RKI besprochen, aber mit den ungeglätteten Werten für den Nowcast aus dem Nowcast-Excel des RKI kann man diese zurückrechnen: Für mehr als 30 Tage zurückliegende Daten gibt es keinen Nowcast-Aufschlag mehr, also sind die Imputierten Werte da gleich den Nowcast-Ergebnissen. Die Änderung des geglätteten Wertes von einem Tag auf den Anderen beträgt gerade ein viertel der Differenz der vier Tage auseinanderliegenden Werte (da ja einer ins gleitende Mittel neu dazukommt und der andere herausfällt). Dann bekommt man eine - grobe, denn mit mit einem Faktor 4 skalierten Fehlern - Schätzung der ungeglätteten imputierten Werte.. Der (implizite, weil mit anderen Methoden berechnete) Schätzer $\hat F(3) \approx 38\%$ des RKI scheint also vergleichbar. Das Nowcasting verstärkt aber die Abweichung in der Imputation, die in den gemittelten Imputations-Werten noch überschaubar war.

Eine Frage, die man sich stellen kann, ist, ob es günstiger wäre, die imputierten Fälle mit zuberücksichtigen. Dies würde zu einer Umgewichtung der "Rastergruppen" führen. Zumindest auf dem Datenstand 17. Mai scheint uns das aber nicht nennenswert näher an die RKI-Ergebnisse zu bringen, obgleich ich aufgrund der Beschreibung, dem Umstand, dass die verschiedenen Imputationspfade dann natürlich leicht zu berücksichtigen sind, und dem R-Paket surveillance für Nowcasting und viele andere infektionsepidemiologische BerechnungenDem R-Paket tun wir unrecht, wenn wir nur für das Nowcasting auf es verweisen. Tatsächlich handelt es sich um einen umfangreichen Werkzeugkasten für die Überwachung von Ausbrüchen von Infektionskrankheiten., das unter anderem von M. Höhle geschrieben wurde, davon ausgehe, dass das RKI das (implizit) macht. Daher habe ich das auch vorgesehen.

Im Code ist diese Anpassung durch Bestimmen der Umgewichtungsfaktoren und Skalierung der Verzögerungen hinterlegt. Dies hat auch den Vorteil, dass man beim Nowcasting für einzelne Bundesländer die Basis-Verzögerungsverteilung von Deutschland benutzen kann und diese dann an die demographische Verteilung der gemeldeten Fälle des Bundeslands anpasst.

factors = (imputed_x.sum(level=(1,2)) / delay_known_new_cases_per_refdatum.sum(level=(0,1))).fillna(0)+1
delay_imputed_cases = delay_known_new_cases_per_refdatum * factors

Das eigentliche Nowcasting ist dann, wie man erwarten würde, in der Methode von Lawless sehr einfach.

def nowcast_lawless(known_cases, imputed, delay_cases):
    assert delay_cases.columns.name == "Refdatum"
    counts = delay_cases.iloc[:, -8-30:].sum(level=2).T[::-1].values

    cumulative_over_days_per_delay = numpy.concatenate([numpy.zeros_like(counts[:1]), counts.cumsum(0)], axis=0)
    counts8 = numpy.tril(cumulative_over_days_per_delay[8:]-cumulative_over_days_per_delay[:-8]) # sum of 8 preceeding days
    n_x = numpy.diagonal(counts8)
    N_x = counts8.sum(1)    

    g_hat = n_x / N_x

    F_hat = (1-g_hat[::-1]).cumprod()[::-1]
    imputed_total = known_cases.add(imputed, fill_value=0)
    nowcast = imputed_total.copy()
    for i in range(1, len(F_hat)):
        nowcast.iloc[-i] /= F_hat[i]
    return nowcast

Die folgende Abbildung zeigt die Ergebnisse unseres Nowcasts verglichen mit jenem des RKI, wie zuvor auf den geglätteten Ergebnissen.

Wie erwartet stimmen unsere Ergebnisse nicht mit jenen des RKI überein. Die könnte an den verwendeten Daten, wahrscheinlicher aber an dem anderen Verfahren liegen, dass wir benutzt haben.

Da die Änderung des Verfahrens in Arbeit ist, verzichten wir für den Moment auf Vorhersageintervalle (diese kommen mit dem neuen Verfahren ohnehin).

Fazit

Wie man sieht, scheint es möglich, die RKI-Berechnungen näherungsweise Anhand der Beschreibung und der Daten nachzuvollziehen. Dass das RKI keinen Code veröffentlicht, kann ich gut verstehen, zumal dieser möglicherweise (wahrscheinlich?) auf Einzelfall-Listen arbeitet. Ich habe auch eine Weile mit mir gerungen, ob ich meinen Code veröffentlichen möchte. Gleichzeitig ist das Erstellen von Codes nicht unbedingt die Aufgabe des RKI (wenn Sie einen teuren Berater dafür benötigen, sagen Sie bescheid).

Zur Datenlage: Natürlich wären weitere Zwischenergebnisse sehr schön. Wenn ich mir nur eine sache wünschen dürfte, wäre es - es wird die aufmerksame Leserin nicht weiter überraschen - ein Feld "Übermittlungsdatum" / "Eingangsdatum", dass das Übermittlungsdatum an das RKI hinterlegt.

Eine Leseempfehlung für zur Berechnung der effektiven Reproduktionszahl habe ich auch noch: M. Höhle: Effective reproduction number estimation vergleicht verschiedene Methoden, die Reproduktionszahl aus Infektionszahlen zu bestimmen und insbesondere die Rolle der Annahme, dass das serielle Intervall nicht fix ist, sondern die Dauer zwischen zwei Erkrankungen in der Infektionskette zufällig entlang einer Verteilung.

Zurück zur Datenlage

Klarerweise wäre es alles viel einfacher, wenn wir ein Übermittlungsdatum in der RKI_COVID19.csv hätten. Haben wir aber nicht. Es sei denn...

...es sei denn, wir bauen es selbst ein. Aus dem obigen haben wir schon gesehen, dass man dafür eine Historie von RKI_COVID19.csv-Abzügen benötigt. Damit kann man es aber - nicht nur wie für die unbekannten Fälle - angehen. Wir versuchen also, Fälle zu durch die verschiedenen Abzüge zu verfolgen. Dafür ist es nützlich, sich zu überlegen, wie man die Beobachtungen günstig aufteit. Das RKI verwendet ja AnzahlFall für alle Fälle und AnzahlTodesfall und AnzahlGenesen für jene Fälle, bei denen die betreffende Person verstorben oder (mutmaßlich) genesen ist. Das klingt plausibel, für uns ist es aber günstiger, in drei disjunkte Kollektive aufzuteilen: Aktive Fälle, Todesfälle und Genesene. Zwischen den Kollektiven sind Übergänge möglich, die natürlichlichen von aktiv zu versorben oder genesen sowie korrekturbedingt auch alle übrigen. Es können neue Fälle zu allen Kollektiven zugeordnet werden und, wieder durch Korrekturen, auch Fälle gelöscht werden. Diese möglichen VorfälleDie Kollektive und Vorfälle lassen den Aktuar durchscheinen. zusammen mit den "nichts passiert"-Fortschreibungen von einem Tag auf den nächsten lassen sich durch die AnzahlFall-Spalte zusammen mit den Merkmalen NeuerFall, NeuerTodesfall sowie Neugenesen rekonstruierenMithin können also die Spalten AnzahlTodesfall und AnzahlGenesen aus den anderen Rekonstruiert werden.:

In jedem Abzug sind also zwei Datenstände und die Übergänge codiert: Jener des Vortages und der aktuelle sowie die Übergänge zwischen Ihnen. Wollen wir nun die Brücke zwischen zwei aufeinanderfolgenden Abzügen schlagen, müssen wir den finalen Stand des Vortagesabzugs mit dem im aktuellen Abzug codierten Vortagesstand verknüpfen. Wenn wir das können, können wir am ersten Tag, den ein Fall auftritt (NeuerFall = 1), das Übermittlungsdatum vermerken und den Fall schritt für schritt über die Zeit verfolgen bis wir beim aktuellen Datenstand ankommen.

Für die verknüpfung Teilen wir die Fälle nach den uns zur Verfügung stehenden Merkmalen auf: Landkreis, Altersgruppe, Geschlecht, IstErkrankungsbeginn und Refdatum (Altersgruppe2 ignorieren wir, das sie nur sporadisch vorliegt). Wie so oft bei Daten ist der große Teil leicht und der Teufel steckt im Detail der letzten paar Fälle. Leicht zu verknüpfenAuch hier werden aber Annahmen getroffen, die nicht unbedingt zutreffen müssen, es könnten durch Korrekturen "Tausche" stattfinden, die sich in den Gesamtzahlen aufheben. Außerdem beachte man, dass hier die gelöschten Fälle des aktuellen Tages (die ja am Vortag enthalten waren) mitabgebildet werden. sind jene Fälle, bei denen der Fall in derselben Kategorie bleibt. Das sind auch die allermeisten, beim Übergang vom 3. zum 4. Juni beispielsweise 181980 der 182370 Fälle (die bis zum 3. Juni 0 Uhr bekannt waren).

Eine sehr plausible Änderung der Kategorie ist das Nachtragen eines Erkrankungsbeginns - wenn zuvor kein Erkrankungsbeginn bekannt war, z.B. weil noch keine Symptome aufgetreten waren oder er schlicht noch nicht erhoben wurde. Von den 390 verbleibenden Fällen in obigem Beispiel lassen sich 285 Fälle mit hinterlegtem Erkrankungsbeginn mit unzugeordnete Fälle ohne hinterlegten Erkrankungsbeginn verknüpfen. Aber was tun mit den verbleibenden 105 hartnäckigen Fällen?

Es stellt sich heraus, dass diese Fälle gegebenenfalls zwischen benachbarten Landkreisen verschoben wurden oder z.B. der Erkrankungsbeginn korrigiert wurde. Da eine manuelle Zuordnung sehr arbeitsintensiv wäre, gehen wir einen etwas anderen Weg. Wir definieren zunächst einen "Abstand" zwischen den obigen Kategorien. Je unplausibler wir finden, dass es eine Zuordnung zwischen Fällen geben könnte, desto größer soll der Abstand sein.

Ich habe mich entschlossen, die Einzelnen Aspekte der Kategorisierung durchzugehen und jeweils "Punkte" für Abweichungen verteilen:

Nun können wir für die verbleibenden Fälle (aufgeteilt in aktive Fälle, verstorbene und genesene Personen) jene Zuordnung suchen, die den Gesamtabstand minimiert.Für Bildungsprotzer: Wir minimieren die (1-Wasserstein-Distanz)[https://de.wikipedia.org/wiki/Wasserstein-Metrik] oder Earth Mover's Distance bzw. Lösen das Transportproblem. Glücklicherweise gibt es dafür eine sehr gute Python-Bibliothek. Damit haben wir auch die verbleibenden Fälle abgebildet.

Jetzt, da wir also Fälle zuordnen können, müssen wir noch das Übermittlungsdatum fortschreiben. Es bietet sich an, auf Kategorienebene zu arbeiten. Dann zählen wir für jede Kategorie die Übermittlungsdaten. Da wir nicht für alle Fälle Daten haben, gibt es für die älteren Fälle auch eine Spalte für unbekannt.Man könnte natürlich auch hier ein Imputationsverfahren benutzen, es kommt uns aber auf die Fälle im März und April nicht so an. Da wir nicht genau wissen, Fälle mit welchem Übermittlungsdatum gelöscht werden, kürzen wir anteilig alle Übermittlungsdaten, ebenso Verfahren wir bei wechsel der Kategorie im obigen Sinne.

Jetzt haben wir jede Menge spalten mit Übermittlungsdaten. Wäre es nicht einfacher, ein einzelnes Übermittlungsdatum zu haben? Das können wir, allerdings gibt es dann durch die Kürzungen für gelöschte Fälle und andere Zuordnungsschwierigkeiten gebrochene Fallzahlen.

Ich habe mich entschlossen, nur die RKI_COVID19 ab 30. 4. (wenn es Erkrankungsbeginne gibt) zu verwenden. Bestimmt könnte man auch weiter zurückgehen, allerdings würde ich das Refdatum (das vorher offenbar schonmal auftaucht) für die älteren Daten ignorieren.

Hier ist das Script für die Datenveredelung und eine Datei (der zweiten Sorte) gibt es auch.

Nowcasts für Bundesländer - Ergebnisse

Hier sind die Ergebnisse unseres Nowcasting. Wie erwartet, sind die Zahlen nahe jenen des RKI, aber nicht identisch. Bitte beachten Sie obige ausführliche Erläuterungen zu den Abweichungen. Die Hauptursache scheint das unterschiedliche Nowcasting-Verfahren zu sein (wir arbeiten daran). Außerdem ist es immer möglich, dass Fehler in der Programmierung zu falschen Resultaten führen. Nehmen Sie den Coda (unten) und prüfen Sie ihn, bevor Sie mir die Zahlen glauben.

Es gibt auch eine CSV-Datei mit Gemeldeten, Imputierten und Nowcast-Ergebnissen. Bitte beachten Sie die obigen Hinweise zur Nutzung dieser Daten. Ich versuche, diese aktuell zu halten. Ältere Versionen sind, soweit vorhanden, unter 2020-03-18 statt latest abgelegt. Die Daten beruhen auf den vorgestellten Verfahren und öffentlichen Daten des RKI.

Wir können das auch für Bayern ansehen. Dazu habe ich die Imputations- und Nowcasting-Verteilungen auf Deutschland-Basis benutzt. Auch andere Bundesländer sind unter ihrem Kürzel verfügbar.

Schließlich gibt es noch auf besonderen Wunsch eine Grafik für die Bundesländer im Norden:

Schließlich, für die volle Transparenz, gibt es auch noch den Code als Jupyter-Notebook und Python-Script. Für Rückmeldungen und Fehlerkorrekturen bin ich sehr dankbar. Wie man sieht, benötigt der Code die RKI_COVID19.csv für alte Datenstände (z.B. ab 29. April). Falls diese jemand öffentlich zur Verfügung stellen mag, bin ich sehr dankbar.