Śródroczny model zapadalności na grypę

Cel opracowania

Na moment pisania niniejszego opracowania od ponad szesciu miesiecy trwa epidemia koronawirusa. U jej początków, dzienna liczba zachorowań kształtowala się na poziomie ok 200-400 przypadków dziennie. W miarę zbliżania się do miesięcy jesiennych, liczba ta zaczęła rosnąć do ponad 3000 tysiecy przypadków dziennie. Większość państw, w tym Polska, podjęło szereg drastycznych środków, by ograniczyć rozmiary epidemii koronawirusa.

Jednak nawet oficjalne dane o zdrowiu publicznym, a także badania naukowe nad wirusami grypy skłaniają do zadania dwóch podstawowych pytań:

• czy zjawisko zachorowań na SARS-CoV-2 można w ogóle zaklasyfikowac jako epidemię?
• jesli nie, to czy obecne zachorowania nie są stałym zjawiskiem, nie diagnozowanym wcześniej.

Aby spróbować odpowiedzieć na te pytania (z pominięciem kwestii, jaka jest definicja epidemii wg WHO - i czy w ogóle istenieje takowa) warto:

• zbadać typową krzywą przebiegu zachorowań na grypę w ciągu roku,
• na tej podstawie zbudować model śródrocznych zachorowań na grypę oraz
• porównać dane o obecnych zachorowaniach na koronawirusa z danymi historycznymi i z modelem.

Dla zainteresowanych prześledzeniem poniższych rozważań, udostępniam kod w języku R służący do pobrania danych, ich przetworzenia i zaprezentowania.

Źródła danych

Badanie typowej krzywej przebiegu zachorowań na grypę w ciągu roku zostanie przeprowadzone na podstawie danych Państwowego Zakładu Higieny Narodowego Instytutu Zdrowia Publicznego. Dane te udostępnione są w formie dwutygodniowych raportów w formacie PDF na poprzedniej wersji strony tej instytucji. W 2020 roku NIZP rozpoczął udostępnianie danych poprzez serwisy EpiBaza i ProfiBaza (w trakcie rozwoju), lecz dane te są zagregowane i nie obejmują okresów śródrocznych. Stąd aby ustalić śródroczny przebieg zapadalności na wybrane choroby zakaźne w Polsce, trzeba sięgnąć do poprzedniej wersji strony. Dane do niniejszego opracowania zostały pobrane ze strony: http://wwwold.pzh.gov.pl/oldpage/epimeld/index_p.html

Zestaw danych zostanie przygotowany w dwóch krokach:

1. dane zostaną pobrane w formie plików PDF
2. z pobranych archiwów zostaną wyodrębnione dane dotyczące poszczególnych chorób.

# Wymagane biblioteki
library(rvest)
library(dplyr)

# Pobranie całego archiwum meldunków o zachorowaniach na choroby zakaźne, zakażeniach i zatruciach w Polsce (dwutygodniowe, kwartalne, półroczne, roczne)
url <- "http://wwwold.pzh.gov.pl/oldpage/epimeld/index_p.html"
pzholdwp <- read_html(url)
links <- pzholdwp %>% html_nodes("a") %>% html_attr("href")
links <- grep("index_mp.html$", links, value = TRUE)
links <- paste0("http://wwwold.pzh.gov.pl/oldpage/epimeld/",links)
yrs <- as.data.frame(stringr::str_extract(links, "\\d{4}"))
l2p <- data.frame()
for (i in 1:length(yrs[,1])){
url2 <- links[i]
pzhYR <- read_html(url2)
links2pdf <- pzhYR %>% html_nodes("a") %>% html_attr("href")
links2pdf <- grep("pdf$", links2pdf, value = TRUE)
links2pdf <- as.data.frame(paste0("http://wwwold.pzh.gov.pl/oldpage/epimeld/",yrs[i,],"/",links2pdf))
l2p <- rbind(l2p, links2pdf)
}
dir.create("pdf")
pth <- paste0(getwd(),"/pdf/")
for (i in 1:length(l2p[,1])){
fname <- stringr::str_sub(l2p[i,], -15, -1)
fname <- stringr::str_remove(fname, "[//]")
httr::GET(l2p[i,], httr::authenticate(":", ":", type="ntlm"),
          httr::write_disk(paste0(pth,fname), overwrite = FALSE),
          httr::progress()
          )
}

Na moment sporządzania tego raportu, na stronach PZH znajduje się 637 meldunków o zachorowaniach na choroby zakaźne, zakażeniach i zatruciach w Polsce (dwutygodniowe, kwartalne, półroczne, roczne). Zestaw danych będzie oparty o raporty dwutygodniowe, by dokładniej zamodelować trend śródroczny.

# "Scrapping" plików pdf
library(pdftools)
library(stringr)
# library(lubridate)
library(plotly)

lista_plikow <- list.files(path = "pdf/")
grypa <- data.frame()

for (i in 1:length(lista_plikow)){
text_rep_temp <- pdf_text(pdf = paste0("pdf/",lista_plikow[i]))
text_rep_temp <- c(unlist(strsplit(text_rep_temp, "\r\n")))
# write(text_rep_temp, file = "test.txt")
wanted_lines <- c(grep("Grypa",text_rep_temp)[1],
                  grep("zgłoszonych w okresie od |o zgłoszonych zachorowaniach za okres od",text_rep_temp)[1]
                  )
wanted <- text_rep_temp[get("wanted_lines")]
data_raportu <- str_extract(wanted[2], "\\d{2}.\\d{2}.\\d{4}")
data_raportu <- lubridate::dmy(data_raportu)
wanted[1]
ilosc_przyp <- as.numeric(str_extract(wanted[1], "\\s\\d+"))
wiersz <- data.frame(data_raportu, ilosc_przyp) 
grypa <- rbind(grypa, wiersz)
}

grypa <- grypa[complete.cases(grypa),]


#################################
# Dane po 2008 - zmiana formatu raportu (data jest we właściwościach pliku, dane o zachorowanian są prezentowane narastająco)

grypa08 <- data.frame()
for (i in 1:length(lista_plikow)){
text_rep_temp08 <- pdf_text(pdf = paste0("pdf/",lista_plikow[i]))
text_rep_temp08 <- c(unlist(strsplit(text_rep_temp08, "\r\n")))
wanted_lines08 <- c(grep("Grypa i podejrzenia",text_rep_temp08)[1]
                  )
wanted08_data_raportu <- pdf_info(pdf = paste0("pdf/",lista_plikow[i]))$keys$Title
wanted08 <- text_rep_temp08[get("wanted_lines08")]
data_raportu08 <- str_extract(wanted08_data_raportu, "\\d{2}.\\d{2}\\s\\d{4}|\\d{2}.\\d{2}.\\d{4}")
data_raportu08 <- lubridate::dmy(data_raportu08)
ilosc_przyp08 <- as.numeric(gsub("\\s","",(str_extract(wanted08[1], "(\\s\\d+){1,3}"))))
wiersz08 <- data.frame(data_raportu08, ilosc_przyp08) 
grypa08 <- rbind(grypa08, wiersz08)
}
grypa08 <- grypa08[complete.cases(grypa08),]

# obliczenie ilości nowych przypadków (uwaga na dane z końca roku)
nowe_pryp_08 <- diff(grypa08$ilosc_przyp08) 
nowe_pryp_08 <- c(grypa08[1,2],nowe_pryp_08)
grypa08 <- cbind(grypa08, nowe_pryp_08)

#poprawa danych z końca roku i pierwszego raportu z początku nowego roku
temp_ds <- subset(grypa08, nowe_pryp_08 <0)
temp_ds$nowe_pryp_08 <- temp_ds$ilosc_przyp08
grypa08a <- subset(grypa08, nowe_pryp_08 > 0)
grypa08a <- rbind(grypa08a, temp_ds)
grypa08a <- grypa08a[order(grypa08a$data_raportu08),]
grypa08a <- grypa08a[,c(1,3)]


# Połączenie danych do 2008 i po, w całość
cn <- colnames(grypa)
colnames(grypa08a) <- cn
grypa_calosc <- rbind(grypa, grypa08a)
grypa_calosc <- grypa_calosc[order(grypa_calosc$data_raportu),]

Po zebranu danych z plików pdf pobranych ze stron PZH, został utworzony zestaw danych zawierający ilość przypadków zachorowań na grypę w interwałach dwutygodniowych. Dane zawierają 580 obserwacji za okres od 1996-01-01 do 2020-09-30.

Analiza rozpoznawcza

Pierwsze spojrzenie na dane wyodrębnione z plików pdf:

Sezonowość dająca się zaobserwować jest oczywista. Niemniej odnośnie samego zestawu danych, trzeba zauważyć:

1. dane ze stycznia 1996 roku zawierają najprawdopodoniej narastające dane z 1995 roku (wskazuje na to ilość zaraportowanych przypadków: 1 775 089),
2. brak jest danych z pierwszej połowy 2009 roku,
3. wzorzec raportowania danych w drugiej połowie lat 90tych XX w. i w pierwszych latach XXI w. wskazywał skokowy wzrost zachorowań w styczniu. Natomiast od roku 2010 widać, że ilość zaraportowanych przypadków rośnie stopniowo.

Z powyższego wynika, że do opracowania śródrocznego modelu zachorowań na grypę należy użyć danych po roku 2009, gdyż są one dokładniejsze.

Trendy długoterminowe prezentuje wykres rocznej ilości przypadków:

Wydaje się, że do dalszej analizy śródrocznego przebiegu zachorowań warto pobrać w pierwszym podejściu dane od roku 2013, w którym widać istotnie większą ilość (ok 3 mln) niż w latach poprzednich, kiedy roczna ilość przypadków kształtowała sie na pomiędzy 300 tys. a 1.5 mln.

Łącząc dwie powyższe obserwacje (dokładniejszy model raportowania po roku 2009 oraz większą liczbę rocznych zachorowań po roku 2013) dalsza analiza będzie oparta na danych z lat 2010-2019.

Śródroczny model zachorowań na grypę

Normalizacja danych do zakresu 0 - 1 będzie pierwszym krokiem budowy modelu. Normalizacja danych pozwoli przeanalizować przebieg zachorowań w ciągu roku z różnych lat, niezależnie od łącznej liczby przypadków. Dodatkową operacją na danych będzie przenumerowanie miesięcy w ten sposów, aby połączyć połówki lat w jeden sezon. Lipiec jest pierwszym miesiącem sezonu, czerwiec ostatnim. Normalizacja zostanie przeprowadzona na danych w obrębie jednego sezonu.

# Normalizacja danych

library(dplyr)
normalizacja <- function(x) {return ((x - min(x)) / (max(x) - min(x)))}
grypa_calosc_rrmm <- grypa_calosc %>% group_by(rok, msc) %>% summarise(suma_przypadkow = sum(ilosc_przyp))
grypa_calosc_rrmm$nr_miesiaca_sezon <- ifelse(grypa_calosc_rrmm$msc<=6, grypa_calosc_rrmm$msc+6, grypa_calosc_rrmm$msc-6)
grypa_calosc_rrmm <- grypa_calosc_rrmm %>% mutate(etykieta_sezonu = 
                                                    ifelse(nr_miesiaca_sezon <= 6, 
                                                           paste0(rok, "/",rok+1),
                                                           paste0(rok-1, "/",rok)
                                                           ))
grypa_calosc_rrmm$msc_nazwa <- month(grypa_calosc_rrmm$msc, label = TRUE)

norm_data <- grypa_calosc_rrmm %>% group_by(etykieta_sezonu) %>% summarise(suma_przypadkow_norm = normalizacja(suma_przypadkow))

grypa_calosc_rrmm$sum_przy_norm <- norm_data$suma_przypadkow_norm

Wybrane dane z lat 2008-2019 pokazują następujący przebieg zachorowalności na grypę w ciągu roku:

Jak widać, śródroczny przebieg grypy w poszczególnych sezonach - jak wspomniano na wstępnie - przypomina prawoskośną krzywą dzwonowatą rozkładu normalnego. Jednak poszczególne lata dość znacząco się od siebie różniły - szczególnie w miesiącach wrzesień - grudzień:

Miesiąc	Współczynnik zmienności
wrz	0.38
paź	0.29
lis	0.28
gru	0.27

Zaprezentowany współczynnik zmienności to stosunek odchylenia standardowego do średniej. Zachorowania we wrześniu w poszczególnch latach różniły się od swojej średniej 38%. Niemniej ogólny przebieg sezonu grypowego w tych latach wskazuje na występowanie wyraźnego wzorca.

Drugim krokiem będzie dobranie odpowiedniej funkcji i jej parametrów, która jest najlepiej dopasowana do znormalizowanych danych. Na zestawie danych określonym powyżej, został zastosowany algorym Quantile Random Forest (metoda numeryczna), dzięki której po normalizacji wyznaczono typowy przebieg sezonu grypowego. ¹ Ilustruje to poniższy wykres:

Z powyższego wykresu widać, że model dość dobrze odzwierciedla przebieg sezonu. W początkowej fazie wartości są niedoszacowane, jednak od miesiąca 6 do 10 znajdują się mniej więcej w środku obserwowanych wartości. Potwierdzają to także podstawowe parametry modelu (należy je czytać w kontekście przytoczonego powyżej współczynnika zmienności):

mtry	RMSE	Rsquared	MAE	RMSESD	RsquaredSD	MAESD
2	0.2598494	0.4511912	0.1803596	0.024928	0.0736318	0.0155539

Wartość Rsquared na poziomie 0.45 oznacza, że model wyjaśnia badane zjawisko w 45%. Wydawałoby się, że jest to słaby wynik. Jednak biorąc pod uwagę współczynnik zmienności kształtujący się pomiędzy 30% a 40% w poszczególnych miesiącach, wartość 45% dla modelu można uznać za satysfakcjonującą.

Prognoza QRF na sezon 2020 / 2021

Ponieważ cel niniejszego opracowania to budowa Śródrocznego modelu zachorowalności na grypę, nie zaś modelu wyjaśniającego czynniki decydujące o zachrowalności, prognoza na sezon 2020/2021 będzie wykonana uproszczonymi metodami. Wyjaśnienie przyczyn zjawiska zachorowalności to osobne zagadnienie. Prezentowane tutaj podejście jest podejściem czysto numerycznym, zatem dane zostaną potraktowane jako zwykłe szeregi czasowe.

Aby móc wykorzystać zbudowany model do prognozy zachrowań na grypę w sezonie 2020/2021 należy najpierw założyć, lub przewidzieć maksymalny poziom zachorowań w danym sezonie (najgorszy miesiąc w sezonie). Jest to potrzebne, by powrócić od danych znormalizowanych (z przedziału 0-1) do danych w skali pierwotnej (rzeczywistej). Dalej są zaprezentowane trzy scenariusze, by wyznaczyć maksymalny poziom zachorowań w sezonie. Podobnie jak w przypadku budowy śródrocznego modelu, prognoza będzie oparta na danych z lat 2010-2019. Do progozowania został użyty pakiet forecast a w jego ramach model ETS (Exponential smoothing state space model)². Poniższa tabela prezentuje prognozę maksymalnej rocznej (sezonowej) ilości przypadków wraz z przedziałami ufności:

	Point Forecast	Lo 80	Hi 80
11	847383	521187.5	1173578

Prognozowana ilość zachorowań w szczycie sezonu 2020/2021 to 847,3 tys. osób. Z prawdopodobieństwem 80% można oczekiwać, że liczba zachorowań w szczycie, będzie mieścić się w przedziale: 521,2 do 1 173,6 tys. osób. Posługując się tą prognozą oraz modelem wyznaczonym w poprzednim punkcie, mozna wyznaczyć prawdopodobny przebieg sezonu grypowego:

Zgodnie z numerycznym modelem QRF, sezon grypowy 2020/2021 będzie miał swój szczyt w lutym 2021. Łaczna ilość zachorowań według wariantu pośredniego wyniesie 3 588 814.

Szczegółową prognozę w miesiącach prezentuje poniższa tabela:

Miesiąc	Dzienne zachorowania, wariant minimalny	Miesięcznie, wariant minimalny	Dzienne zachorowania, wariant pośredni	Miesięcznie, wariant pośredni	Dzienne zachorowania, wariant maksymalny	Miesięcznie, wariant maksymalny
lip	1362	40854	1401	42043	1441	43233
sie	1303	39096	1303	39096	1303	39096
wrz	1726	51767	2011	60340	2297	68914
paź	4988	149640	7481	224437	9974	299235
lis	7026	210780	10898	326946	14770	443112
gru	7508	225228	11706	351170	15904	477112
sty	11733	352002	18791	563723	25848	775443
lut	17373	521187	28246	847383	39119	1173578
mar	14238	427146	22990	689712	31743	952276
kwi	6315	189457	9707	291196	13098	392934
maj	2104	63118	2646	79372	3188	95626
cze	1985	59554	2447	73396	2908	87238

Wnioski dotyczące SARS-CoV-2

Odnosząc powyższe rozważania odnośnie zachorowalności na grypę, na epidemię koronawirusa, należy przytoczyć najważniejsze badania w tym zakresie. Badania te pozwolą oszacować, jaka ilość zachorowań może w rzeczywistości dotyczyć koronawirusa. Warto poznać:

1. “Virus–virus interactions impact the population dynamics of influenza and the common cold”, PNAS December 26, 2019 116 (52) 27142-27150; first published December 16, 2019, Sema Nickbakhsh, Colette Mair, Louise Matthews, Richard Reeve, Paul C. D. Johnson, Fiona Thorburn, Beatrix von Wissmann, Arlene Reynolds, James McMenamin, Rory N. Gunson, and Pablo R. Murcia PNAS
2. “Epidemiology of Seasonal Coronaviruses: Establishing the Context for the Emergence of Coronavirus Disease 2019”, The Journal of Infectious Diseases, 01 Jun 2020, 222(1):17-25, EUROPE PMC

Najważniejszy wniosek płynący z tych badań to:

udział koronawirusów w corocznych zachorowaniach na grypę wynosi od 7 do 15 procent

Dzięki tym badaniom oraz zbudowanej prognozie, można w prosty sposób wyznaczyć prawodopodobną liczbę zachorowań na koronawirusa zarówno w latach ubiegłych, jak i w sezonie 2020/2021. Przy wszystkich zastrzeżeniach, że obecny szczep koronawirusa jest inny od poprzednich, należy zauważyć, że nawet “sezonowe” koronawirusy różnią się od siebie. Biorąc pod uwagę nowoczesną wiedzę molekularną o wirusach, można z dość dużą pewnością stwierdzić, że koronawirus z początków 2020 roku zdążył w ciągu kilku miesięcy zmutować i już różni się od początkowej swojej wersji.

Sezon	Suma zachorowań	Zachorowania na koronawirusa (15%)	Zachorowania na koronawirusa (7%)
2010/2011	1055300	158295	738710
2011/2012	1040955	156143	728668
2012/2013	2959318	443898	2071523
2013/2014	2724400	408660	1907080
2014/2015	3786125	567919	2650288
2015/2016	4050970	607646	2835679
2016/2017	4803745	720562	3362622
2017/2018	5431585	814738	3802109
2018/2019	4635847	695377	3245093
2019/2020	4103763	615564	2872634

Jak wynika z tej prostej kalkulacji, w sezonie 2019/2020 liczba chorych na koronawirusa wynosła pomiędzy 287,2 tys. a 615,5 tys. co daje średnią dzienną zachorowalność w ciągu sezonu między ok 800 a 1670 przypadków dziennie.

Model zachrowalności na grypę użyty do prognozy ilości przypadków koronawirusa w sezonie 2020/2021 wskazuje na następujący teoretyczny przebieg zachorowań (wg 15% udziału w infekcjach):

Według modelu pośredniego, we wrześniu powinno być wykryte 51,8 tys. przypadków grypy. w tym 8,8 tys. przypadków koronawirusa. Jak wynika z danych PZH we wrześniu wykryto 217,6 tys. przypadków grypy, w tym 26,1 tys. przypadków koronawirusa:

### Budowa zestawu danych dla sezonu 2020/2021 

grypa20 <- data.frame()
for (i in 1:length(lista_plikow)){
text_rep_temp20 <- pdf_text(pdf = paste0("pdf/",lista_plikow[i]))
text_rep_temp20 <- c(unlist(strsplit(text_rep_temp20, "\r\n")))
wanted_lines20 <- c(
                  grep("Zakażenia SARS-CoV-2",text_rep_temp20)[1]
                  )
wanted_lines20
wanted20_data_raportu <- pdf_info(pdf = paste0("pdf/",lista_plikow[i]))$keys$Title
wanted20 <- text_rep_temp20[get("wanted_lines20")]
data_raportu20 <- str_extract(wanted20_data_raportu, "\\d{2}.\\d{2}\\s\\d{4}|\\d{2}.\\d{2}.\\d{4}")
data_raportu20 <- lubridate::dmy(data_raportu20)
ilosc_przyp20 <- as.numeric(gsub("\\s","",(str_extract(wanted20[1], "(\\s\\d+){1,3}"))))
wiersz20 <- data.frame(data_raportu20, ilosc_przyp20) 
grypa20 <- rbind(grypa20, wiersz20)
}
grypa20 <- grypa20[complete.cases(grypa20),]

# obliczenie ilości nowych przypadków (covida, dane od lipca)
nowe_pryp_20 <- diff(grypa20$ilosc_przyp20) 
nowe_pryp_20 <- c(grypa20[1,2],nowe_pryp_20)
grypa20 <- cbind(grypa20, nowe_pryp_20)
grypa20 <- subset(grypa20, grypa20$data_raportu20 > "2020-06-30")
grypa20$msc <- month(grypa20$data_raportu20)
grypa20 <- grypa20 %>% group_by(msc) %>% summarise(suma_covid = sum(nowe_pryp_20))
gry2021 <- subset(grypa_calosc_rrmm, grypa_calosc_rrmm$etykieta_sezonu=="2020/2021"&grypa_calosc_rrmm$msc>6)
gry2021 <- gry2021[,c(1:6)]
gry2021 <- cbind(gry2021, grypa20$suma_covid)
gry2021$suma_przypadkow_z_cov <- gry2021$suma_przypadkow + gry2021$...7

Miesiąc	Suma przypadków	Grypa: pozostałe szczepy	Grypa: koronawirus
lip	51331	39096	12235
sie	92087	70564	21523
wrz	217602	191476	26126

Udział koronawirusa w łącznej liczbie przypadków wynosi 12% i jest zgodny z obserwacjami i badaniami przytoczonymi powyżej: “Virus–virus interactions impact the population…” oraz “Epidemiology of Seasonal Coronaviruses: …”.

Rozwój zachorowań na grypę w sezonie 2020/2021 w porównaniu do modelu i ostatnich 5 lat wygląda na wykresie następująco:

Na powyższym wykresie wyraźnie widać dwie rzeczy:

1. model teoretyczny zedycowanie niedoszacowuje rozwój zachorowań w początkowych miesiącach sezonu,
2. sezon grypowy 2020/2021 jest zdecydowanie łagodniejszy niż sezony 2017/2018, 2018/2019, 2019/2020.

Na koniec niniejszego opracowania przytoczę jeszcze ostatni wykres, podobny jak poprzedni, tyle, że zawierający dane odnośnie zachorowań na grypę spowodowanych koronawirusem:

Dane historyczne odnośnie zachorowań są obliczone przy założeniu 15% udziału tego szczepu wirusa w ogólnych zachorowaniach. Sezon 2020/2021 jest nie tylko łagodniejszy (udział koronawirusa wynosi 12%) ale i ogólna liczba zachorowań jest znacząco mniejsza niż w latach ubiegłych.

---

---

1. Journal of Machine Learning Research 7 (2006) 983–99, Nicolai Meinshausen “Quantile Regression Forests”.↩

2. International Journal of Forecasting, Volume 18, Issue 3, July–September 2002, Pages 439-454, Rob J Hyndman, Anne B Koehler, Ralph D Snyder, Simone Grosea↩