W tym dokumencie wizualizujemy zawartość rejestru zgonów z Powodu COVID19 w PL udostępnionego przez MZ pod adresem https://basiw.mz.gov.pl/index.html. Porównujemy ten rejestr z danymi z bazy danych nt. pandemii COVID z John Hopkins University.

Baza JHU

Baza danych zgonów COVID z John Hopkins University (liczba zgonów dziennie); filter wyciąga wiersze z danymi dla Polski; select wyciąga tylko kolumny date oraz newd (dzienna liczba zgonów z powodu COVID):

jhu_0 <- read.csv("covid19_JHU.csv", sep = ';',  header=T, na.string="NA") %>%
  filter ( country == 'Poland') %>%
  select (date, newd)

Pierwszy i ostatni dzień

first.jhu.day <- first(jhu_0$date)
last.jhu.day <- last(jhu_0$date)

last.jhu.day
## [1] 2021-12-24
## 724 Levels: 2020-01-01 2020-01-02 2020-01-03 2020-01-04 ... 2021-12-24

Baza BASiW

Baza BASiW (https://basiw.mz.gov.pl/index.html) zgonów z powodu COVID za rok 2020; Baza zawiera dane indywidualne, ale w bardzo dziwacznym formacie, mianowicie dla każdej kombinacji dnia, województwa, wieku i płci zawiera liczbę zgonów z powodu covid oraz z powodu covid i innych chorób. Ten dziwaczny format zamieniono na dane indywidualne (odpowiednio powielając wiersze). Dodatkowym utrudnieniem jest to, że kolumna wiek zawiera symbol 90+ (który nie jest liczbą)

Z tego drugiego powodu czytam wiek jako typ znakowy; następnie zamieniamy 95+ arbitralnie na 95. Następnie zamieniamy typ znakowy na liczbę (as.integer). Przekodowujemy zmienną zgony_19 (T zgon z powodu covid19; N zgon z powodu covid + inne choroby)

basiw2020 <- read.csv("basiw_2020.csv", sep = ';',  header=T,
                      colClasses=c('wojewodztwo' = 'character', 'wiek'='character'),
                      na.string="NA") %>%
  mutate ( wiek = case_when(wiek == '90+' ~ '95', wiek != '90+' ~ wiek )) %>%
select(data, woj=wojewodztwo, zgony_c19,  plec, vek=wiek) %>%
mutate (
  vek = as.integer(vek),
  woj = as.factor(woj),
  zgony_c19 = case_when(
    zgony_c19 == 0 ~ "N",
    zgony_c19 == 1 ~ "T" )) %>%
mutate ( zgony_c19 = as.factor(zgony_c19))

Baza BASiW rok 2021 podobnie jak w roku 2020 zawiera dane pseudo-indywidualne, ale trochę inne; Podobnie jak poprzednio najpierw zamieniane skryptem na dane idywidualne (w sensie jeden wiersz = jeden zgon)

Nie ma nazwy województwa (jak w roku 2020), ale jest nr_teryt. Płeć zawiera dwie wartości K, M + nieznana.

W związku z tym: teryt_woj importujemy jako character (żeby 02 nie zmieniło się w 2); Płeć też importujemy jako character, co pozwoli usunąć nienana i dzięki temu będzie to samo co w bazie w roku 2020.

basiw2021_0 <- read.csv("basiw_2021.csv", sep = ';',  header=T, 
                      colClasses=c('teryt_woj' = 'character', 'plec' = 'character',
                                   'wiek'='integer'),
                      na.string="NA") %>%
  mutate ( plec = 
             case_when(plec == 'nieznana' ~ NA_character_, 
                       plec != 'nieznana' ~ plec )) %>%
  select(data=data_rap_zgonu, woj=teryt_woj, 
         zgony_c19=czy_wspolistniejace, plec, vek=wiek) %>%
  mutate (plec = as.factor(plec))

Uwaga na mutate zmieniające wartości zmiennej płeć. Wszystkie wartości po znaku ~ w poleceniu case_when muszą być tego samego typu; dlatego wstawiamy NA_character_ a nie (beztypowe) NA.

Wczytujemy nazwy województw + numery teryt ze stosownego pliku csv:

woj <- read.csv("teryt_wojewodztwa_BB.csv", sep = ';',  header=T, 
                colClasses=c('teryt' = 'character'), na.string="NA")

Łączymy basiw2021 z ramką nazw województw:

basiw2021 <- left_join(basiw2021_0, woj, by=c('woj'='teryt')) %>%
  select(data, woj=nazwa, zgony_c19, plec, vek)

Wreszcie wszystkie kolumny w ramkach basiw2021 oraz basiw2020 są tego samego typu; łączymy ramki (w pionie). Nawiasem mówiąc, że należało wykonać aż tyle transformacji żeby otrzymać połączony i ujednolicony plik dla lat 2020 i 2021 jest miarą burdelu i niekompetencji w MZ. Mi by było wstyd publikować dane aż tak na odpierdol, no ale im nie jest:

basiw <- rbind(basiw2020, basiw2021)

## sprawdzamy co wyszło:
first.basiw.day <- first(basiw$data)
last.basiw.day <- last(basiw$data)
last.basiw.day
## [1] 2021-12-28
## 626 Levels: 2020-03-12 2020-03-13 2020-03-14 2020-03-16 ... 2021-12-28

Ustalmy która baza (basiw czy jhu) jest któtsza:

last.day <- last.jhu.day
if (as.Date(last.basiw.day) <= as.Date(last.jhu.day)) {
  last.day <- last.basiw.day
} 

last.day 
## [1] 2021-12-24
## 724 Levels: 2020-01-01 2020-01-02 2020-01-03 2020-01-04 ... 2021-12-24

last.day to starsza data z dwóch ostatnich;

Sprawdzamy ile mamy:

## zabici łącznie
total.deaths <- nrow(basiw)

Obcinamy bazy baziw i jhu do ostatniego wspólnego dnia' (czyli do dnialast.day`) a następnie agregujemy wyniki:

## BASIW database ###
## dane dzienne
basiw_d <- basiw %>% 
  filter (as.Date(data) <= as.Date(last.day) ) %>%
  group_by (data) %>%
  summarise ( n = n() )  %>%
  ungroup() 

## dane tygodniowe
basiw_w <- basiw %>%
  mutate(
    ww = substr(date2ISOweek(data), 1,8)) %>%
  group_by (ww) %>%
  summarise ( n = n() )  %>%
  ungroup() %>%
  mutate(dat = ISOweek2date(sprintf("%s-1", ww))) %>%
  ungroup() 

total.deaths.w <- sum(basiw_w$n)

## JHU database ###
## ta baza już zawiera dane dzienne

jhu <- jhu_0 %>%
  filter (as.Date(date) <= as.Date(last.day) )

## dane tygodniowe

jhu_w <- jhu %>%
  mutate(
     ww = substr(date2ISOweek(date), 1,8))  %>%
  group_by (ww) %>%
  summarise ( n.jhu = sum(newd) )  %>%
  ungroup() %>%
  mutate(dat = ISOweek2date(sprintf("%s-1", ww))) %>%
  ungroup()

Wykresy

Zgony dziennie baza basiw

p1 <- ggplot(basiw_d, aes(x= as.Date(data, format="%Y-%m-%d"), y=n)) + 
  geom_smooth(method="loess", se=F, span=spanV) +
  geom_line(size=.4, alpha=.5) +
  geom_point(size=1, alpha=.5) +
  xlab(label="") +
  ylab(label="age") +
  scale_x_date( labels = date_format("%y/%W"), breaks = mainBreaks) +
  theme_nikw()

p1

Zgony tygodniowe baza basiw

p2 <- ggplot(basiw_w, aes(x= as.Date(dat, format="%Y-%m-%d"), y=n)) + 
  geom_smooth(method="loess", se=F, span=spanV) +
  geom_line(size=.4, alpha=.5) +
  geom_point(size=1, alpha=.5) +
  xlab(label="") +
  ylab(label="age") +
  scale_x_date( labels = date_format("%y/%W"), breaks = mainBreaks) +
  theme_nikw()

p2

Porównanie basiw vs jhu (dane tygodniowe)

basiw_vs_jhu <- left_join(basiw_w, jhu_w, by='dat')

p3 <- basiw_vs_jhu %>% 
  pivot_longer(cols=c('n', 'n.jhu'), names_to = 'base', values_to = 'v') %>%
  ggplot (aes(x= as.Date(dat, format="%Y-%m-%d"), y=v, color=base)) + 
  ##geom_smooth(method="loess", se=F, span=spanV) +
  geom_line(size=.4, alpha=.25) +
  geom_point(size=1, alpha=.25) +
  xlab(label="") +
  ylab(label="age") +
  scale_x_date( labels = date_format("%y/%W"), breaks = mainBreaks) +
  theme_nikw()

p3
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_point).

ggsave(plot=p3, file='basiw_vs_jhu.png', width = 9)
## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).
p4 <- basiw_vs_jhu %>% 
  mutate  (d = n - n.jhu) %>%
  ggplot (aes(x= as.Date(dat, format="%Y-%m-%d"), y=d)) + 
  ##geom_smooth(method="loess", se=F, span=spanV) +
  geom_line(size=.4, alpha=.25) +
  geom_point(size=1, alpha=.25) +
  xlab(label="") +
  ylab(label="age") +
  scale_x_date( labels = date_format("%y/%W"), breaks = mainBreaks) +
  theme_nikw()

p4
## Warning: Removed 1 row(s) containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_point).

Generalnie się zgadza, co jest chyba jedynym pozytywem z całego tego ćwiczenia.