Epidemie

Pół szklanki mocnego kodu

Zanim dopadnie nas grypa

Czy lubicie długoterminowe prognozy pogody? Odbija się w nich głęboko zakorzeniona ludzka wiara, że odległą przyszłość można przewidzieć - na przekór efektowi motyla. No cóż, sam muszę przyznać, że cieszę się, gdy grudniowe prognozy przewidują piękną, słoneczną i mroźną aurę na zimowe ferie; ale jeśli zapowiadają ponury mokry standard, wtedy ratuję się myślą, że to przecież tylko prognoza...

Oczywiście można przewidywać nie tylko pogodę: na przykład zimą byłoby całkiem ciekawie przepowiedzieć rozwój corocznej epidemii grypy. I właśnie teraz, mimo braku wiedzy medycznej, uzbrojony jedynie w serię liczb, prosty model i kilka linii kodu (wybiorę moje ulubione środowisko obliczeniowe Octave) spróbuję przewidzieć, jak będzie przebiegać w tym sezonie.

Szczerze mówiąc, zdziwiłbym się, gdyby przyszłość posłuchała się moich prognoz: użyte przeze mnie narzędzia są raczej mało wyrafinowane. Ale... gdyby przewidywania się sprawdziły... czyż nie byłoby zabawnie?

Skorzystamy z publicznie dostępnych danych NIZP-PZH na temat raportowanej przez lekarzy w całej Polsce liczby zachorowań na grypę (jak je zdobyć samemu, pisaliśmy w Delcie 10/2019). Zaczniemy od naszkicowania wykresu zachorowań w ostatnich latach, zobacz rysunek obok. (Dane zbierane są cztery razy w miesiącu).

Wykres jest mało regularny, niemniej możemy zauważyć, że w pierwszym kwartale każdego roku ma silny pik, a w sierpniu - wyraźne minimum. Tabelka obok pokazuje daty, gdy zostało osiągnięte ekstremum, a dane do niej można wygenerować kodem:

W programie przyjęliśmy, że w zmiennej infected znajdują się liczby zachorowań, które wcześniej ściągnęliśmy z NIZP-PZH, a w days - daty odpowiadających pomiarów. Funkcja findpeaks z pakietu signal znajduje położenia loc lokalnych maksimów danych zawartych w wektorze infected. Dodatkowo wymagamy, by odległość między maksimami wynosiła co najmniej 24 pomiary, czyli 6 miesięcy (jasne, są przecież odległe z grubsza o rok). W trzeciej linii wydobywamy z danych tylko te odpowiadające znalezionym maksimom, które w kolejnej linii wypisujemy na ekranie.

Minimum zawsze przypada na drugi tydzień sierpnia. Maksimum nie zachowuje się aż tak regularnie, ale jak możemy obliczyć, średnia odległość między maksimami na przestrzeni ostatnich 6 lat wynosi około 368 dni - czyli z grubsza też jeden rok (co w końcu nie powinno nas aż tak dziwić...).

Rzuca się w oczy, że od kilku lat piki zgłaszanej liczby zachorowań systematycznie rosną (niemal liniowo). Nie będąc lekarzem, ale za to trochę wyobrażając sobie, jak działa sprawozdawczość, przypuszczałbym, że za tym wzrostem nie stoi realny wzrost liczby wszystkich zachorowań na grypę w Polsce - ale że zapewne jest on efektem coraz lepiej działającego systemu zbierania danych (coraz więcej lekarzy odsyła stosowne formularze?). Gdyby więc liniowy trend miał utrzymać się jeszcze w tym sezonie, łatwo będzie przewidzieć orientacyjną szczytową liczbę zgłoszeń:

Całą robotę - dopasowanie wykresu prostej do maksimów - wykonuje pierwsza linia powyższego kodu: funkcja polyfit wyznacza metodą najmniejszych kwadratów wielomian $P$ stopnia pierwszego, najmniej oddalony od punktów pomiarowych $|(xk,yk),$ gdzie $xk =$ maxdays(k), a $yk =$ maxinfected(k). W drugiej linii obliczamy wartość tego wielomianu w dniu następnego maksimum zachorowań (przyjęliśmy, że wypadnie dokładnie za rok od ostatniego).

Naturalnie w dłuższej perspektywie liniowy wzrost nie może się utrzymać, co zrozumiemy, przedłużając kolorowa prostą w lewo: "przewiduje" ona, że już w 2008 roku lekarze powinni byli zgłosić... ujemną liczbę zachorowań.

Modelowanie w wybranym sezonie

Powyższe miało charakter czysto rozrywkowy, więc czas na coś poważniejszego. Prawie każdy sezon chorobowy ma podobny przebieg (rysunek obok):

dosyć szybki wzrost zachorowań od połowy sierpnia,
z początkiem października następuje stabilizacja liczby zachorowań, a nawet - w pierwszym tygodniu stycznia - jej chwilowy spadek,
od początku stycznia - kilkutygodniowa faza szybkiego wzrostu liczby zachorowań,
osiągnięcie szczytowej liczby zachorowań (najczęściej w lutym),
szybkie zmniejszanie się liczby zachorowań gdzieś tak do końca kwietnia,
powolny spadek do kolejnego minimum.

Czy więc dałoby się - dysponując bieżącymi danymi ze strony internetowej NIZP-PZH - przewidzieć rozwój tegorocznej grypy? Do tego wykorzystamy następujący, bardzo prosty, model zachorowań. Niech $|St$ oznacza liczbę osób podatnych na wirusa (ale na razie zdrowych) w dniu $t,$ z kolei $It$ - liczbę osób zarażonych (i przy tym zarażających), a $R t$ - liczbę osób, które już nie mogą być zarażone ani zarażać, bo np. wyzdrowiały i nabrały odporności. W modelu przyjmujemy następującą zależność $|S,I,R$ w dniu następnym od tychże w dniu obecnym:

$pict$

Gdybyśmy znali współczynniki modelu: $a, b,$ to moglibyśmy go użyć do wyznaczenia przyszłych wartości $S, I,R$ tak, jak w poniższej funkcji:

To bardzo zgrubny model, ale mamy nadzieję, że okaże się wystarczający, jeśli ograniczymy zakres jego stosowalności. Nie możemy liczyć na to, że uda się nam w modelu dopuszczającym tylko jeden pik zachorowań odtworzyć wieloletni przebieg choroby. Ale wymodelowanie okresu największej zachorowalności w sezonie, powiedzmy: od początku stycznia do końca kwietnia, wydaje się celem jak najbardziej realnym.

Ponieważ $1/a$ odpowiada średniej długości okresu, gdy zainfekowany zaraża, to możemy z dobrym skutkiem przyjąć, że $1/a = 2,5$ dnia (potem bierzemy zwolnienie i kładziemy się do łóżka), więc $a = 0,4$ (na dzień).

Nie znając dobrej wartości współczynnika $b,$ dopasujemy ją tak, by wartości $|I$ wyznaczone przez model jak najlepiej przybliżały te, które już znamy z raportów NIZP-PZH.

W tym celu definiujemy funkcję

b ( F(b) = [I1(b),I2(b),...,In(b)],

która dla zadanego $|b$ wyznacza (za pomocą powyższego modelu) wartości $|It$ w punktach pomiarów w zadanym przez nas okresie. Następnie, korzystając z funkcji leasqr z pakietu optim, znajdziemy wartość $b,$ która zminimalizuje sumę kwadratów odchyleń wartości przewidywanych od wartości zmierzonych:

(I1(b) −infected(1))2 + ...+(In(b) − infected(n))2 min!

Najpierw ograniczamy zbiór danych infected, days do interesującego nas okresu (implementację funkcji przytnijdane można sobie z grubsza wyobrazić), a po załadowaniu pakietu optim funkcja leasqr wyznacza $| b$ minimalizujące sumę kwadratów odchyleń (czwarty argument, $F,$ jest potrzebny, bo ta funkcja wyznacza wartości $I1(b), I2(b),...$ potrzebne do obliczenia odchyleń) - omówienie jej implementacji, wykorzystującej KMcK, pominiemy.

Wymodelowany przebieg grypy dla sezonów 14/15, 15/16 i 16/17. Pionowa kreska wskazuje realne maksimum zachorowań. Linia kolorem: wynik działania modelu dla $|b$ wyznaczonego na podstawie wszystkich danych. Czarna, przerywana linia: dla $b$ wyznaczonego na podstawie tylko pierwszych czterech pomiarów

Jak widać na wykresach powyżej, dla tak wyznaczonego $| b$ dostaliśmy kolorową krzywą, całkiem dobrze wpasowującą się w punkty pomiarowe (sezon 2016/2017 jest - jak widzieliśmy wcześniej - pod wieloma względami wyjątkowy i model nie przewiduje prawidłowo wielkości piku).

Okazuje się, że dla większości sezonów dostajemy podobną wartość optymalnego $b.$ Aby więc przewidywać przyszłe epidemie, możemy po prostu skorzystać ze "średniej" wartości $| b ≈ 0,434,$ co zamyka sprawę (po prostu startujemy symulację z wartości znanych w konkretnym dniu). Z drugiej strony wybór uśrednionego $b$ powoduje, że nie możemy reagować na ewentualne atypowości bieżącego sezonu (np. niesprzyjającą aurę lub szczególnie zjadliwy szczep wirusa, co w każdym przypadku powoduje zmianę współczynnika $b$ ). Wówczas lepiej dopasować nowe $|b$ do bieżących danych. Przykładowy wynik dopasowania $b$ na podstawie jedynie czterech pomiarów ze stycznia (linia przerywana na wykresach obok) pokazuje, że możemy się wtedy sporo pomylić; jednak im więcej pomiarów wykorzystamy, tym lepsze będą przewidywania.

A zatem... zaglądając do mojej szklanej/krzemowej kuli...

przepowiadam, że:

Następny sezon grypy... już się zaczął: w połowie sierpnia!
Od października do końca grudnia raportowana liczba zachorowań nie przekroczy... 140 tysięcy w tygodniu! [Odczytany liniowy trend z pierwszego wykresu].
Szczyt zachorowań przypadnie... w drugim tygodniu lutego! [Po prostu zakładam, że stanie się to za rok od ostatniego]. Będzie wtedy zgłaszanych aż 400 tysięcy zachorowań tygodniowo! [Z liniowego trendu na pierwszym wykresie; jednak nie wierzę w to zbytnio i mam nadzieję, że naprawdę będzie ich znacznie mniej... Chociaż, z drugiej strony, grypa jest chorobą zdradziecką, więc po poprzednim roku ciszy teraz może zaatakować ze zdwojoną siłą].
Przebieg zachorowalności w I kwartale 2020 roku... będę na żywo monitorować i prognozować na mojej stronie internetowej już od teraz, co tydzień uwzględniając nowe dane! [Będę dopasowywać na bieżąco parametry modelu, minimalizując błąd średniokwadratowy, podobnie jak opisano powyżej].

No dobra, czas kończyć pisanie. Już w trakcie poczułem się - jakoś tak - niewyraźnie: boli mnie głowa, drapie w gardle i... tak jakby... łamie w kościach? Czyżby -

Aaa... apsik!...