Covid19 marimekko chart

Dane

library(ggplot2)
library(dplyr)

owid <- 'https://covid.ourworldindata.org/data/owid-covid-data.csv'
dat <- "2021/02/20"
default_violet <- '#C77CFF'
default_cyan <- '#00BFC4'

Dane dotyczące zgonów z powodu COVID19 pobrano z https://covid.ourworldindata.org/data/owid-covid-data.csv (szczegóły za chwilę). Plik fao_names.csv zawiera nazwy krajów także pobrane z bazy FAO, do których dodani (za pomocą skryptu Perlowego) następnie dodano do nich kody państw (https://pl.wikipedia.org/wiki/ISO_3166-1_alfa-2)

n <- read.csv("fao_names.csv", sep = ';',  header=T, na.string="NA");

Także z bazy FAO pobrano liczbę ludności. Obie ramki połączono. FAO ma problem polityczny z Chinami, w szczególności nie istnieje tam Tajwan więc niezbędne są pewne poprawki:

popbase <- read.csv("FaoPopBase2018.csv", sep = ';',  header=T, na.string="NA");
n <- left_join(n, popbase, by="country") %>% 
  filter (country != 'China') %>%
  mutate(country = recode(country,
                          `Bolivia (Plurinational State of)`='Bolivia',
                          `Micronesia (Federated States of)`='Micronesia',
                          `United Kingdom of Great Britain and Northern Ireland`='UK',
                          `Venezuela (Bolivarian Republic of)`='Venezuela',
                          `China Hong Kong SAR`='Hong Kong',
                          `China Macao SAR`='Macao',
                          `China Taiwan Province of`='Taiwan',
                          `China mainland`='China',
                          `Iran (Islamic Republic of)`='Iran',
                          `Democratic People's Republic of Korea`='North Korea',
                          `United Republic of Tanzania`='Tanzania'
  ))

Teraz można wczytać dane dotyczące COVID (najbardziej aktualne dotyczą 2021-02-20) Dane te też są zanieczyszone zawierają mianowicie agregaty typu Europa czy Unia Europejska. Należy te agregaty usunąć, żeby się potem nie sumowało trzy razy to samo:

d <- read.csv("covid_20210220.csv", sep = ';',  header=T, na.string="NA");
## należy usunąćbo mącą agregaty
aggregates <- c('Africa', 'Asia', 'Europe', 'European Union', 
                'International', 'North America', 'Oceania', 
                'South America', 'World')
d <- left_join(d, n, by='iso') %>% filter (! country.x %in% aggregates )

Liczymy współczynniki zgonów na 1mln mieszkańców (dr). Ponieważ kolumna pop zawiera liczbę ludności w tysiącach to żeby przejść na miliony należy oczywiście podzielić wartości totald wartości przez tysiąc a nie przez milion:

## ponieważ w bazie FAO ludność jest w tys
million <- 1000

d <- d %>% 
  select(totalc, totald, country.y, iso, pop) %>%
  mutate(popm = pop / million) %>%
  ## Oblicz współczynniki na 1mln
  mutate(cr = totalc/popm, dr = totald/popm) %>%
  ## Tylko kraje wykazujące zmarłych
  filter(totald > 0) %>% as.data.frame
totalWC <- sum(d$totald, na.rm = T)
totalPop <- sum (d$pop, na.rm=T)

Liczba krajów: 180. Łączna liczba zgonów na świecie 2460804 (2021/02/20)

Ponieważ wykres dla wszystkich krajów byłby nieczytelny ograniczymy się do krajów które reprezentują 90% wszystkich zgonów (w porządku malejącym):

## 90% przypadków na świecie
d90 <- d %>%  
  ## oblicz udziały
  mutate(totaldf = totald/totalWC *100 ) %>%
  ## uporządkuj malejąco
  arrange(desc(totaldf)) %>%
  ## policz szereg kumulowany
  mutate(totaldC = cumsum(totaldf)) %>% 
  filter (totaldC < 90)
c90n <- nrow(d90)
c90df <- sum(d90$totaldf)

Takich krajów jest 32 (reprezentują dokładnie 89.54% zgonów na świecie i mieszka tam łącznie 3700342.218 osób czyli 49.14%).

Wykres marimekko

Marimekko to wykres pokazujący jednocześnie zjawisko w ujęciu jednostkowym i globalnym. Taki histogram, tyle że szerokości słupków są proporcjonalne do jednej cechy (a nie stałe) a wysokości są proporcjonalne do wielkości drugiej cechy. Jedną cechą jest zwykle jakiś iloraz a drugą mianownik tego ilorazu. W przypadku wizualizacji zgonów: wysokości słupków to liczba zgonów na milion, podczas gdy szerokości to liczba mieszkańców w kraju. Pole słupka jest zatem proporcjonalne do liczby zgonów.

Żeby wykreślić takie marimekko trzeba dokonać elementarnych przekształceń:

d90 <- d90 %>% arrange (dr) %>% 
  mutate( ww = cumsum(popm)) %>%
  mutate( wm= ww - popm) %>%
  mutate( wt = wm + (ww - wm)/2) %>%
  mutate (isoh = as.character(iso), isol = as.character(iso), isoll=as.character(iso))

Zmienne isoh/isol będą potrzebne do wstawienia etykiet w dwu wierszach. Otóż gdyby etykiety wstawić w jednym to by na siebie zachodziły. Żeby tego nie było niektóre przesuwamy do drugiego wiersza (por wykres); niestety trzeba to robić ręcznie:

d90$isoh[ (d90$iso %in% c('CL', 'EC', 'HU', 'SE') ) ] <- ""
d90$isol[ (! d90$iso %in% c('CL', 'EC', 'HU', 'SE') ) ] <- ""

p9  <- ggplot(d90, aes(ymin = 0)) +
  ylab(label="Covid deaths per 1mln") +
  xlab(label="population (mln)") +
  ggtitle(sprintf("Deaths from COVID19"), 
          subtitle=sprintf("%i krajów z największą liczbą zgonów (łącznie %.1f%% zgonów na świecie) \nsource %s (%s)", 
                           c90n, c90df, owid, dat)) +
  geom_rect(aes(xmin = wm, xmax = ww, ymax = dr, fill = iso)) +
  ##
  geom_text(aes(x = wt, y = 0, label = isoh), vjust=+0.5, hjust=+1.25, size=2.0, angle = 90) +
  geom_text(aes(x = wt, y = 0, label = isol), vjust=+0.5, hjust=2.95, size=2.0, angle = 90 ) + 
  ##
  ylim(-100, NA) +
  theme(legend.position = "none")
p9

Wykres wygląda efektownie, ale cierpi na wszystkie wady wykresów powierzchniowych mianowicie małe różnice są niewidoczne. Owszem dokładnie widać które kraje mają jakie wielkości zgonów na milion ale jakie są różnice w wielkościach zgonów to już nie jest tak dobrze widoczne, np jak się ma łączna liczba zgonów w UK i w Meksyku…

Zamiast marimekko

Więc oczywiście zawsze można zamiast marimekko wyrysować obie cechy oddzielnie:

library("ggpubr")
## zgony na 1mln
p0 <- d90 %>% 
  ggplot(aes(x= reorder(as.factor(country.y), dr) , y=-dr )) +
  geom_bar(alpha=.3, fill=default_cyan, stat="identity") +
  xlab(label="") +
  ylab(label="deaths") +
  ggtitle(sprintf("Deaths from COVID19 (per 1mln)"), 
          subtitle=sprintf("source %s (%s)", owid, dat)) +
  theme(axis.text = element_text(size = 8), 
        axis.ticks.y = element_blank(),
        axis.text.y=element_blank()) +
  theme(plot.title = element_text(hjust = 0.5)) +
  ##scale_y_continuous(breaks=seq(0,30, by=2)) +
  coord_flip()
## udziały w całości
p1 <- d90 %>% 
  ggplot(aes(x= reorder(as.factor(country.y), dr) , y=totaldf )) +
  geom_bar(alpha=.3, fill=default_cyan, stat="identity") +
  xlab(label="") +
  ylab(label="%") +
  ggtitle(sprintf("Deaths from COVID19 (%% of world total)"), 
          subtitle=sprintf("source %s (%s)", owid, dat)) +
  theme(axis.text = element_text(size = 8)) +
  theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(breaks=seq(0,30, by=2)) +
  coord_flip()

## https://rpkgs.datanovia.com/ggpubr/reference/ggarrange.html
ggarrange(p0, p1, ncol=2, widths = c(2.5, 3) )

Albo nawet na jednym ale to wymaga pewnej ekwilibrystyki. Np możemy przedstawić wielkość zgonów na 1mln jako wielkość względną, konkretnie względem wielkości maksymalnej. Wtedy wszystkie słupki odpowiadające zgonom/1mln będą się zawierały w przedziale 0–100. Wielkości udziałów w zgonach też są w tym przedziale:

### Normalize to 100
### zi = (xi – min(x)) / (max(x) – min(x)) * 100
drMin <-  min(d90$dr)
drMax <-  max(d90$dr)
d90$drn10 <- (d90$dr - drMin)/ (drMax - drMin) * 10

p2 <- d90 %>% 
  ggplot(aes(x= reorder(as.factor(country.y), totaldf) , y=totaldf )) +
  geom_bar(size=2.5, alpha=.3, fill=default_cyan, stat="identity") +
  geom_point(aes(y=drn10), size=2.5, alpha=.6, color=default_violet) +
  xlab(label="") +
  ylab(label="% of world total") +
  ggtitle(sprintf("Deaths from COVID19 (%% of world total)"), subtitle=sprintf("source %s (%s)", owid, dat)) +
  theme(axis.text = element_text(size = 8)) +
  ##theme(plot.title = element_text(hjust = 0.5)) +
  scale_y_continuous(breaks=seq(0,30, by=2), 
                     sec.axis = sec_axis(~ . * 10, name='death ratio % (BE=100)')) +
  coord_flip()

p2

Dorysowujemy drugą oś (sec.axis) o skali równej pierwszej przemnożonej przez 10 (wynik widać na wykresie). Druga cecha jest wizualizowana kropkami, żeby się słupki nie nakładały na siebie…

Zakończenie

Szczegółowe dane dla 32 krajów są następujące (%zgonów to udział danego kraju w zgonach ogółem–nie sumuje się do 100 bo, krajów jest tylko 32)

library(knitr)
d90 <- d90 %>% arrange(totald) %>% select(iso, country.y, totald, dr, pop, totaldf)
kable(d90,  
      col.names = c('kod', 'kraj', 'zgony', 'zgony/1mln', 'ludność', '%zgonów'))

kod	kraj	zgony	zgony/1mln	ludność	%zgonów
PH	Philippines	12068	113.15370	106651.394	0.4904088
PK	Pakistan	12601	59.37474	212228.286	0.5120684
SE	Sweden	12649	1268.49771	9971.638	0.5140190
IQ	Iraq	13245	344.62033	38433.600	0.5382387
HU	Hungary	14252	1468.14334	9707.499	0.5791603
NL	Netherlands	15323	898.20605	17059.560	0.6226827
EC	Ecuador	15513	908.02359	17084.358	0.6304037
PT	Portugal	15897	1549.99033	10256.193	0.6460084
CZ	Czechia	19097	1790.50988	10665.677	0.7760472
RO	Romania	19795	1014.81002	19506.114	0.8044119
CL	Chile	19974	1066.46534	18729.160	0.8116859
CA	Canada	21631	583.44587	37074.562	0.8790217
BE	Belgium	21887	1906.17146	11482.178	0.8894248
UA	Ukraine	26404	596.75240	44246.156	1.0729827
TR	Turkey	27983	339.84661	82340.088	1.1371487
ID	Indonesia	34316	128.20238	267670.543	1.3945036
PL	Poland	42077	1109.57894	37921.592	1.7098883
PE	Peru	44690	1397.03138	31989.260	1.8160731
ZA	South Africa	48940	846.82242	57792.518	1.9887809
AR	Argentina	51122	1152.40475	44361.150	2.0774511
CO	Colombia	58685	1181.71087	49661.048	2.3847897
IR	Iran	59409	726.26972	81800.188	2.4142110
ES	Spain	67101	1437.07202	46692.858	2.7267917
DE	Germany	67883	816.64331	83124.418	2.7585700
RU	Russian Federation	81517	559.35457	145734.038	3.3126165
FR	France	83546	1285.51074	64990.511	3.3950693
IT	Italy	95486	1574.96729	60627.291	3.8802765
GB	UK	120593	1796.09734	67141.684	4.9005528
IN	India	156302	115.55309	1352642.280	6.3516639
MX	Mexico	179797	1424.80289	126190.788	7.3064332
BR	Brazil	245977	1174.28651	209469.323	9.9957981
US	United States of America	497648	1521.41144	327096.265	20.2229840

Ten dokument w formacie Rmd oraz dane są dostępne tutaj: https://github.com/hrpunio/Bookdown/tree/main/Misc

Covid19 marimekko chart

TP

23 02 2021

Dane

Wykres marimekko

Zamiast marimekko

Zakończenie