![[Photo gallery]](/staff/tp/icons/ico_gallery.png){kind=icon}
W mediach podają albo to co było wczoraj, albo razem od początku świata. Zwłaszcza to łącznie od początku ma niewielką wartość, bo są kraje gdzie pandemia wygasa a są gdzie wręcz przeciwnie. Co z tego, że we Włoszech liczba ofiar to 35 tysięcy jak od miesięcy jest to 5--25 osób dziennie czyli tyle co PL, gdzie zmarło do tej pory kilkanaście razy mniej chorych. Podobnie w Chinach, gdzie cała afera się zaczęła, zmarło we wrześniu 15 osób. Media ponadto lubią duże liczby więc zwykle podają zmarłych, a nie zmarłych na 1mln ludności. Nawet w danych ECDC nie ma tego wskaźnika zresztą. Nie ma też wskaźnika zgony/przypadki, który mocno niefachowo i zgrubie pokazywałby skalę zagrożenia (im większa wartość tym więcej zarażonych umiera)
Taka filozofia stoi za skryptem załączonym niżej. Liczy on zmarłych/1mln ludności oraz współczynnik zgony/przypadki dla 14 ostatnich dni. Wartości przedstawia na wykresach: punktowym, histogramie oraz pudełkowym. W podziale na kontynenty i/lub kraje OECD/pozostałe. Ponieważ krajów na świecie jest ponad 200, to nie pokazuje danych dla krajów w których liczba zmarłych jest mniejsza od 100, co eliminuje małe kraje albo takie, w których liczba ofiar COVID19 też jest mała.
Wynik obok. Skrypt poniżej
![[Zdjęcie T. Przechlewskiego]](/personal/rysunki/home-icon.jpg){kind=icon}
#!/usr/bin/env Rscript
#
library("ggplot2")
library("dplyr")
#
## Space before/after \n otherwise \n is ignored
url <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
surl <- sprintf ("retrived %s from %s", tt, url)
## population/continent data:
c <- read.csv("ecdc_countries_names.csv", sep = ';', header=T, na.string="NA" )
## COVID data
d <- read.csv("covid19_C.csv", sep = ';', header=T, na.string="NA",
colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));
# today <- Sys.Date()
# rather newest date in data set:
today <- as.Date(last(d$date))
# 14 days ago
first.day <- today - 14
tt<- format(today, "%Y-%m-%d")
fd<- format(first.day, "%Y-%m-%d")
# select last row for every country
##t <- d %>% group_by(id) %>% top_n(1, date) %>% as.data.frame
## top_n jest obsolete w nowym dplyr
t <- d %>% group_by(id) %>% slice_tail(n=1) %>% as.data.frame
# cont. select id/totald columns, rename totald to alld
t <- t %>% select(id, totald) %>% rename(alld = totald) %>% as.data.frame
# join t to d:
d <- left_join(d, t, by = "id")
## extra column alld = number of deaths at `today'
##str(d)
# only countries with > 99 deaths:
d <- d %>% filter (alld > 99 ) %>% as.data.frame
## OECD membership data
o <- read.csv("OECD-list.csv", sep = ';', header=T, na.string="NA" )
## Remove rows with NA in ID column (just to be sure)
## d <- d %>% filter(!is.na(id))
d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)
cat ("### Cleaning newc/newd: assign NA to negatives...\n")
d$newc[ (d$newc < 0) ] <- NA
d$newd[ (d$newd < 0) ] <- NA
## Filter last 14 days only
d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > fd) %>% as.data.frame
## Last day
last.obs <- last(d$date)
## Format lable firstDay--lastDay (for title):
last14 <- sprintf ("%s--%s", fd, last.obs)
## Sum-up cases and deaths
t <- d %>% group_by(id) %>% summarise( tc = sum(newc, na.rm=TRUE), td = sum(newd, na.rm=TRUE)) %>% as.data.frame
## Add country populations
t <- left_join(t, c, by = "id")
## Add OECD membership
## if membership column is not NA = OECD member
t <- left_join(t, o, by = "id")
t$access[ (t$access > 0) ] <- 'oecd'
t$access[ (is.na(t$access)) ] <- 'non-oecd'
str(t)
## Deaths/Cases ratio %%
t$totr <- t$td/t$tc * 100
## Deaths/1mln
t$tot1m <- t$td/t$popData2019 * 1000000
## PL row
dPL <- t[t$id=="PL",]
##str(dPL)
## Extract ratios for PL
totr.PL <- dPL$totr
totr.PL
tot1m.PL <- dPL$tot1m
## Set scales
pScale <- seq(0,max(t$totr, na.rm=T), by=1)
mScale <- seq(0,max(t$tot1m, na.rm=T), by=10)
## Graphs
## 1. deaths/cases ratio (dot-plot)
## color=as.factor(access)): draw OECD/nonOECD with separate colors
p1 <- ggplot(t, aes(x = reorder(country, totr) )) +
### One group/one color:
##geom_point(aes(y = totr, color='navyblue'), size=1) +
### Groups with different colors:
geom_point(aes(y = totr, color=as.factor(access)), size=1) +
xlab("country") +
ylab("%") +
ggtitle(sprintf("Covid19 death cases ratio (%s)", last14),
subtitle="Only countries with 100 deaths and more | pale red line = PL level") +
theme(axis.text = element_text(size = 4)) +
theme(plot.title = element_text(hjust = 0.5)) +
##theme(legend.position="none") +
scale_color_discrete(name = "Membership") +
## Indicate PL-level of d/c ratio:
geom_hline(yintercept = totr.PL, color="red", alpha=.25, size=1) +
labs(caption=surl) +
scale_y_continuous(breaks=pScale) +
##coord_flip(ylim = c(0, 15))
coord_flip()
ggsave(plot=p1, file="totr_p14.png", width=10)
## 2. deaths/cases ratio (histogram)
p2 <- ggplot(t, aes(x = totr) ) +
geom_histogram(binwidth = 0.5, fill='navyblue', alpha=.5) +
###
ylab("N") +
xlab("%") +
ggtitle(sprintf("Covid19 deaths/cases ratio (%s)", last14),
subtitle="Only countries with 100 deaths and more | pale red line = PL level") +
scale_x_continuous(breaks=pScale) +
scale_y_continuous(breaks=seq(0, 40, by=2)) +
geom_vline(xintercept = totr.PL, color="red", alpha=.25, size=1) +
##coord_cartesian(ylim = c(0, 30), xlim=c(0, 8))
labs(caption=surl)
ggsave(plot=p2, file="totr_h14.png", width=10)
## 3. deaths/1m (dot-plot)
## color=as.factor(continent)): draw continents with separate colors
p3 <- ggplot(t, aes(x = reorder(country, tot1m) )) +
geom_point(aes(y = tot1m, color=as.factor(continent)), size =1 ) +
xlab("") +
ylab("deaths") +
ggtitle(sprintf("Covid19 deaths per 1 million (%s)", last14),
subtitle="Only countries with 100 deaths and more | red pale line = PL level") +
theme(axis.text = element_text(size = 6)) +
geom_hline(yintercept = tot1m.PL, color="red", alpha=.25, size=1) +
labs(caption=surl) +
scale_color_discrete(name = "Continent") +
scale_y_continuous(breaks=mScale) +
coord_flip()
ggsave(plot=p3, file="totr_m14.png", height=10)
## 3. deaths/1m (box-plots for continents)
p4 <- ggplot(t, aes(x=as.factor(continent), y=tot1m, fill=as.factor(continent))) +
geom_boxplot() +
ylab("deaths") +
xlab("continent") +
ggtitle(sprintf("Covid19 deaths per 1 million (%s)", last14),
subtitle="Only countries with 100 deaths and more") +
labs(caption=surl) +
scale_y_continuous(breaks=mScale) +
theme(axis.text = element_text(size = 6)) +
theme(legend.position="none")
ggsave(plot=p4, file="totr_c14.png", width=10)
## Mean/Median/IQR per continent
ts <- t %>% group_by(as.factor(continent)) %>%
summarise(Mean=mean(tot1m), Median=median(tot1m), Iqr=IQR(tot1m))
ts
Wartość zmarli/przypadki zależy ewidentnie od słynnej liczby testów. Taki Jemen wygląda, że wcale nie testuje więc ma mało przypadków. Z drugiej strony ci co dużo testują mają dużo przypadków, w tym słynne bezobjawowe więc też nie bardzo wiadomo co to oznacza w praktyce. Bo chyba się testuje po coś, a nie żeby testować dla samego testowania. Więc: dużo testów -- na mój niefachowy rozum -- to wczesne wykrycie, a mało to (w teorii przynajmniej) dużo nieujawnionych/nieizolowanych zarażonych, co powinno dawać w konsekwencji coraz większą/lawinowo większą liczbę przypadków, a co za tym idzie zgonów a tak nie jest.
Konkretyzując rysunek: Ameryka na czele liczby zgonów z medianą prawie 30/1mln, Europa poniżej 3,5/1mln a Azja coś koło 1,4/1mln (średnie odpowiednio 26,4/6,00/5,1.)
Przy okazji się okazało, że mam dplyra w wersji 0.8 co jest najnowszą wersją
w moim Debianie. No i w tej wersji 0.8 nie ma funkcji slice_head()
oraz slice_tail(), która jest mi teraz potrzebna.
Z instalowaniem pakietów R w Linuksie to różnie bywa,
ale znalazłem
na stronie Faster package installs on Linux with Package Manager beta release coś takiego:
RStudio Package Manager 1.0.12 introduces beta support for the binary format of R packages. In this binary format, the repository's packages are already compiled and can be `installed' much faster -- `installing' the package amounts to just downloading the binary!
In a configured environment, users can access these packages using `install.packages' without any changes to their workflow. For example, if the user is using R 3.5 on Ubuntu 16 (Xenial):
> install.packages("dplyr", repos = "https://demo.rstudiopm.com/cran/__linux__/xenial/latest")
Zadziałało, ale (bo jest zawsze jakieś ale) powtórzenie powyższego na raspberry pi nie zadziała. Wszystko się instaluje tyle że architektury się nie zgadzają. Więc i tak musiałem wykonać
> install.packages("dplyr")
## sprawdzenie wersji pakietów
> sessionInfo()
Instalacja trwała 4 minuty i zakończyła się sukcesem.
Jeszcze jeden skrypt do rysowania danych nt COVID19:
#!/usr/bin/env Rscript
# Przekazywanie argumentów z wiersza poleceń
# np.: Rscript --vanilla c19.R --iso PL -clean
library("optparse")
#
library("ggplot2")
library("dplyr")
library("scales")
library("ggpubr")
#
# parametr wygładzania (loess)
spanV <- 0.25
# UWAGA: przed/po \n musi być odstęp inaczej nie działa
surl <- "https://www.ecdc.europa.eu/en/publications-data/ \n download-todays-data-geographic-distribution-covid-19-cases-worldwide"
c0 <- 'PL'
option_list <- list(
make_option(c("-i", "--iso"), action="store", type="character", default=c0, help="country ISO2 code"),
make_option(c("-c", "--clean"), action="store_true", default=T, help="extra clean data?")
);
opt_parser <- OptionParser(option_list=option_list);
opt <- parse_args(opt_parser);
c0 <- opt$iso
dataClean <- opt$clean
# wczytanie danych
# date;id;country;newc;newd;totalc;totald
d <- read.csv("covid19_C.csv", sep = ';', header=T, na.string="NA",
colClasses = c('factor', 'factor', 'factor', 'character', 'character', 'numeric', 'numeric'));
str(d)
d$newc <- as.numeric(d$newc)
d$newd <- as.numeric(d$newd)
## Dane zawierają wartości ujemne co jest bez sensu
## z opcją --clean te ujemne wartości są zamieniane na NA
if ( dataClean ) {
cat ("### Cleaning newc/newd: assign NA to negatives...\n")
d$newc[ (d$newc < 0) ] <- NA
d$newd[ (d$newd < 0) ] <- NA
}
## Współczynniki zmarli/zakażeni
d$newr <- d$newd/d$newc * 100
d$totr <- d$totald/d$totalc * 100
## Wartości > 50% są zamieniane na NA (zwykle >50% wynika z błędnych danych)
if ( dataClean ) {
cat ("### Cleaning newc/newd: assign NA to newr/totr higher than 50...\n")
d$newr[ (d$newr > 50) ] <- NA
d$totr[ (d$totr > 50) ] <- NA
}
## Pomiń obserwacje wcześniejsze niż 15/02
d <- d %>% filter(as.Date(date, format="%Y-%m-%d") > "2020-02-15") %>% as.data.frame
d0 <- d %>% filter (id == c0) %>% as.data.frame
t0 <- d0 %>% group_by(id) %>% summarise(cc = sum(newc, na.rm=T), dd=sum(newd, na.rm=T))
lab0c <- toString(paste (sep=" = ", t0$id, t0$cc))
lab0d <- toString(paste (sep=" = ", t0$id, t0$dd))
## koniecznie dodać na.rm=T bo inaczej zwraca NA (jeżeli znajdzie NA)
maxCC <- max (d0$newc, na.rm=T)
maxDD <- max (d0$newd, na.rm=T)
maxRR <- max (d0$totr, na.rm=T)
last.obs <- last(d0$date)
first.date <- first(d0$date)
fstDay <- as.Date(first.date)
last.totr <- last(d0$totr)
max.newr <- max(d0$newr, na.rm=T)
## Przykład dodania 14 dni do daty
## srcDay <- as.Date(first.date) +14
## https://stackoverflow.com/questions/10322035/r-adding-days-to-a-date
srcDay <- as.Date(last.obs)
## Nazwa pliku wynikowego
## c19_ISO2_DATA.png, gdzie DATA jest datą ostatniej obserwacji
## np.: c19_SE_2020-09-16.png
c0o <- sprintf ("c19_%s_%s.png", c0, last.obs)
## Rysunek1: nowe przypadki
pc0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newc)) +
geom_point(aes(group = id, color = id), size=.8) +
geom_smooth(method="loess", se=F, span=spanV) +
theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
annotate("text", x = fstDay, y = 0.95 * maxCC,
label = sprintf("Total: %i cases", t0$cc), hjust = 0, vjust=0,
alpha=0.3, color='steelblue', size=6) +
xlab(label="") +
## Nie drukuj legendy
theme(legend.position="none") +
ggtitle(sprintf("%s: new confirmed cases (%s)", t0$id, last.obs))
## Rysunek2: nowe zgony
pd0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newd)) +
geom_point(aes(group = id, color = id), size=.8) +
geom_smooth(method="loess", se=F, span=spanV) +
theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
annotate("text", x = fstDay, y = 0.95 * maxDD, label = sprintf("Total: %i deaths", t0$dd), hjust = 0, vjust=0,
alpha=0.3, color='steelblue', size=6) +
scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
xlab(label="") +
theme(legend.position="none") +
ggtitle(sprintf ("%s: deaths (%s)", t0$id, last.obs))
## Rysunek3: nowe zgony/przypadki *100%
pr0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=newr)) +
geom_point(aes(group = id, color = id), size=.8) +
geom_smooth(method="loess", se=F, span=spanV) +
theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
annotate("text", x = fstDay, y = 0.95 * max.newr, label = sprintf("Maximum: %.2f %%", max.newr), hjust = 0, vjust=0,
alpha=0.3, color='steelblue', size=6) +
scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
xlab(label="") +
ylab(label="%") +
theme(legend.position="none") +
ggtitle(sprintf ("%s: deaths/cases %% (%s)", t0$id, last.obs) )
## Rysunek4: łączne zgony/przypadki *100%
prt0 <- ggplot(d0, aes(x= as.Date(date, format="%Y-%m-%d"), y=totr)) +
geom_point(aes(group = id, color = id), size=.8) +
geom_smooth(method="loess", se=F, span=spanV) +
theme(plot.subtitle=element_text(size=8, hjust=0, color="black")) +
annotate("text", x = fstDay, y = 0.95 * maxRR, label = sprintf("Average: %.2f %%", last.totr), hjust = 0, vjust=0,
alpha=0.3, color='steelblue', size=6) +
scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
xlab(label="") +
ylab(label="%") +
theme(legend.position="none") +
annotate("text", x = srcDay, y = 0, label = surl, hjust = 1, alpha=.3, size=3) +
ggtitle(sprintf ("%s: total deaths/cases %% (%s)", t0$id, last.obs) )
p00 <- ggarrange(pc0,pd0, pr0, prt0, ncol=2, nrow=2)
ggsave(plot=p00, c0o, width=15)
Użycie:
Rscript --vanilla c19.R --iso PL
albo:
for i in 'AU' 'BR' 'IN' 'US' 'ES' 'SE' 'PL' 'DE' 'GB'; Rscript --vanilla c19.R --iso $i done
Przykładowy wynik dla Hiszpanii
Jak widać jakość danych jest katastrofalna: pojawiające się liczne zera to w rzeczywistości brak danych. Zwykle sobota/niedziela zero a potem sruuuu 30 tysięcy bo za trzy dni razem. Wszyscy są zmęczeni w tej Hiszpanii pandemię i nawet nie chce im się danych podsyłać do ECDC?
Reuters donosi, że While Russia has confirmed the world's fourth largest tally of coronavirus cases, it has a relatively low death toll from the associated disease, COVID-19[...] But data released by the Rosstat State Statistics Service on Sept. 4 show there were 57,800 excess deaths between May and July, the peak of the outbreak. The figure was calculated by comparing fatalities over those three months in 2020 with the average number of May-July deaths between 2015 and 2019. The excess total is more than three times greater than the official May-July COVID-19 death toll of 15,955. Co mnie zmotywowało do poszukania danych źródłowych.
Okazało się że dość prosto jest je pobrać ze strony
Urzędu Statystycznego Federacji Rosyjskiej (https://eng.gks.ru/
a konkretnie https://showdata.gks.ru/finder/)
Ponieważ okazało się to aż tak proste, to pobrałem
nie tylko zgony ale także urodzenia. W okresie 2015--2020 (miesięcznie).
Dane są w formacie XSLX. Kolumny to miesiące wiersze poszczególne regiony i ogółem Federacja. Eksportuję do CSV używając LibreOffice i wybieram wiersz z danymi dla całej federacji. W rezultacie mam trzy wiersze: rok-miesiąc, urodzenia i zgony. Obracam moim super skryptem do transpozycji plików CSV:
month;urodzenia;zgony 2015-01-01;149269.99;174722.99 2015-02-01;144591.99;156737 2015-03-01;160974.99;175809.99 ...
Trzy wykresy liniowe przedstawiające dynamikę zgonów, urodzin i przyrostu naturalnego (czyli różnicy między urodzinami a zgonami)
library("dplyr")
library("ggplot2")
library("scales")
spanV <- 0.5
## przyrost naturalny
d <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
d$diff <- d$urodzenia - d$zgony
pz <- ggplot(d, aes(x= as.Date(month), y=zgony)) +
geom_line(color="steelblue", size=.6, alpha=.5) +
geom_point(color="steelblue", size=.8) +
## Trend
geom_smooth(method="loess", se=F, span=spanV, color="red3") +
## Skala na osi X-ów
scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
xlab(label="") +
ylab(label="deaths") +
geom_vline(xintercept = as.Date("2015-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2016-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2017-07-01"), alpha=.25, size=1)+
geom_vline(xintercept = as.Date("2018-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2019-07-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2020-07-01"), alpha=.25, size=1) +
labs(caption="https://showdata.gks.ru/finder/") +
ggtitle("Russian Federation: Deaths (ths)",
subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: july", spanV))
pu <- ggplot(d, aes(x= as.Date(month), y=urodzenia)) +
geom_line(color="steelblue", size=.6, alpha=.5) +
geom_point(color="steelblue", size=.8) +
geom_smooth(method="loess", se=F, span=spanV, color="red3") +
scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
xlab(label="") +
ylab(label="births") +
labs(caption="https://showdata.gks.ru/finder/") +
geom_vline(xintercept = as.Date("2015-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2016-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2017-08-01"), alpha=.25, size=1)+
geom_vline(xintercept = as.Date("2018-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2019-08-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2020-08-01"), alpha=.25, size=1) +
ggtitle("Russian Federation: Births (ths)",
subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: august", spanV))
pp <- ggplot(d, aes(x= as.Date(month), y=diff)) +
geom_line(color="steelblue", size=.6, alpha=.5) +
geom_point(color="steelblue", size=.8) +
geom_smooth(method="loess", se=F, span=spanV, color="red3") +
scale_x_date( labels = date_format("%y/%m"), breaks = "4 months") +
geom_vline(xintercept = as.Date("2015-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2016-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2017-05-01"), alpha=.25, size=1)+
geom_vline(xintercept = as.Date("2018-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2019-05-01"), alpha=.25, size=1) +
geom_vline(xintercept = as.Date("2020-05-01"), alpha=.25, size=1) +
xlab(label="") +
ylab(label="balance") +
labs(caption="https://showdata.gks.ru/finder/") +
ggtitle("Russian Federation: natural balance (ths)",
subtitle= sprintf("red line: loess with span=%.2f ; gray vertical: may", spanV))
Wykres liniowy dla miesięcy jednoimiennych (MJ):
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)
size0 <- 1.2
size1 <- 0.6
size2 <- 0.6
tt <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
tt$rok <- substr(tt$month, 1, 4)
tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month, 6, 7),
substr(tt$month, 9, 10) ), format="%m-%d")
## Obliczenie przyrostu naturalnego
tt$diff <- tt$urodzenia - tt$zgony
## Obliczenie
tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame
pf <- ggplot() +
ggtitle("Russian Federation: births 2015--2020") +
geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2020'), alpha=.5, size=size0) +
geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2020'), alpha=.5, size=size0) +
geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2019'), alpha=.5, size=size1) +
geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2019'), alpha=.5, size=size0) +
geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2018'), alpha=.5, size=size1) +
geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2018'), alpha=.5, size=size1) +
geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2017'), alpha=.5, size=size1) +
geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2017'), alpha=.5, size=size1) +
###
geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2016'), alpha=.5, size=size1) +
geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2016'), alpha=.5, size=size1) +
geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2015'), alpha=.5, size=size2) +
geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = urodzenia, colour = '2015'), alpha=.5, size=size1) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
labs(caption="https://showdata.gks.ru/finder/") +
ylab(label="ths") +
xlab(label="") +
labs(colour = "") +
theme(legend.position="top") +
theme(legend.text=element_text(size=10));
qf <- ggplot() +
ggtitle("Russian Federation: deaths 2015-2020") +
geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony, colour = '2020'), alpha=.5, size=size0) +
geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = zgony, colour = '2020'), alpha=.5, size=size0) +
geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony, colour = '2019'), alpha=.5, size=size1) +
geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = zgony, colour = '2019'), alpha=.5, size=size0) +
geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony, colour = '2018'), alpha=.5, size=size1) +
geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = zgony, colour = '2018'), alpha=.5, size=size1) +
geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony, colour = '2017'), alpha=.5, size=size1) +
geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = zgony, colour = '2017'), alpha=.5, size=size1) +
###
geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony, colour = '2016'), alpha=.5, size=size1) +
geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = zgony, colour = '2016'), alpha=.5, size=size1) +
geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony, colour = '2015'), alpha=.5, size=size2) +
geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = zgony, colour = '2015'), alpha=.5, size=size2) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
labs(caption="https://showdata.gks.ru/finder/") +
ylab(label="ths") +
xlab(label="") +
labs(colour = "") +
theme(legend.position="top") +
theme(legend.text=element_text(size=10));
qf <- ggplot() +
ggtitle("Russian Federation: natural balance 2015-2020") +
geom_line( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff, colour = '2020'), alpha=.5, size=size0) +
geom_point( data = tt.yy2020, mapping = aes(x=yyyymmdd, y = diff, colour = '2020'), alpha=.5, size=size0) +
geom_line( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff, colour = '2019'), alpha=.5, size=size1) +
geom_point( data = tt.yy2019, mapping = aes(x=yyyymmdd, y = diff, colour = '2019'), alpha=.5, size=size0) +
geom_line( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff, colour = '2018'), alpha=.5, size=size1) +
geom_point( data = tt.yy2018, mapping = aes(x=yyyymmdd, y = diff, colour = '2018'), alpha=.5, size=size1) +
geom_line( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff, colour = '2017'), alpha=.5, size=size1) +
geom_point( data = tt.yy2017, mapping = aes(x=yyyymmdd, y = diff, colour = '2017'), alpha=.5, size=size1) +
###
geom_line( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff, colour = '2016'), alpha=.5, size=size1) +
geom_point( data = tt.yy2016, mapping = aes(x=yyyymmdd, y = diff, colour = '2016'), alpha=.5, size=size1) +
geom_line( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff, colour = '2015'), alpha=.5, size=size2) +
geom_point( data = tt.yy2015, mapping = aes(x=yyyymmdd, y = diff, colour = '2015'), alpha=.5, size=size2) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
labs(caption="https://showdata.gks.ru/finder/") +
ylab(label="ths") +
xlab(label="") +
labs(colour = "") +
theme(legend.position="top") +
theme(legend.text=element_text(size=10));
Wreszcie wykres kombinowany liniowo-słupkowy dla miesięcy jednoimiennych (MJ).
Linie to 1) średnia liczby zgonów w latach 2015--2020,
2) średnia liczba zgonów w latach 2015--2020 plus/minus odchylenie
standardowe liczby zgonów (dla MJ); 3) liczba zgonów w roku 2020. Wykres słupkowy to
różnica między liczbą zgonów w roku 2020
a maksymalną liczbą zgonów w latach 2015-2019 ORAZ liczba zgonów z powodu COVID19
(z bazy ECDC). Ten skrypt korzysta z danych z pliku covid_ru.csv, który
powstał przez zagregowanie danych dziennych dostępnych ze strony ECDC (dla Federacji Rosyjskiej)
oraz policzenie średnich/odchyleń
dla jednoimiennych miesięcy z danych w pliku UZRT.csv
library(ggplot2)
library(dplyr)
library(tidyr)
library(scales)
options(dplyr.print_max = 1e9)
#
tt <- read.csv("covid_ru.csv", sep = ';', header=T, na.string="NA");
tt$diff <- tt$urodzenia - tt$zgony
cols <- c("mean19"="navyblue","mean+sd19"="#3591d1", "mean-sd19"="#3591d1",
"deaths20"='brown4', "diff19"="red", "c19d"='cyan')
pf <- ggplot(data=tt) +
ggtitle("Russian Federation: deaths 2015--2020",
subtitle='mean19/sd19: mean/sd 2015--19; diff19: 2020 - max(2015-019); deaths20: deaths 2020; c19d: C19 deaths') +
geom_line(mapping = aes(x=yyyymmdd, y = mean, colour = 'mean19'), alpha=.5, size=.8) +
geom_line(mapping = aes(x=yyyymmdd, y = mean + sd, colour = 'mean+sd19'), alpha=.5, size=.8) +
geom_line(mapping = aes(x=yyyymmdd, y = mean - sd, colour = 'mean-sd19'), alpha=.5, size=.8) +
geom_text(mapping = aes(x=yyyymmdd, label=difflabel, y=diff), vjust=-1.0, size=3 ) +
##geom_line(mapping = aes(x=yyyymmdd, y = c19d, colour = 'c19d'), alpha=.5, size=1.2) +
labs(caption = "Source: https://showdata.gks.ru/finder/; https://www.ecdc.europa.eu/en/covid-19-pandemic") +
##
geom_bar(mapping = aes(x=yyyymmdd, y = diff, fill = 'diff19' ), color='brown', stat="identity", alpha=.25) +
geom_bar(mapping = aes(x=yyyymmdd, y = c19d, fill = 'c19d' ), color ='navyblue', stat="identity", alpha=.25) +
###
geom_line(mapping = aes(x=yyyymmdd, y = deaths, colour = 'deaths20'), alpha=.5, size=1.2) +
scale_x_date( labels = date_format("%y/%m"), breaks = "2 months") +
ylab(label="ths") +
xlab(label="") +
##labs(colour = "Legend") +
scale_colour_manual(name="Lines: ", values=cols) +
scale_fill_manual(name="Bars: ", values=cols) +
theme(legend.position="right") +
theme(legend.text=element_text(size=12));
Skrypt pomocniczy liczący średnie 2015--2019 dla miesięcy jednoimiennych itp rzeczy
library(dplyr)
library(tidyr)
## Agreguje zgodny wg miesięcy jednoimiennych (MJ) dla lat 2015-19
## drukuje średnią/sd/max dla MJ
## Drukuje z2020 - max(z2015--2019)
tt <- read.csv("UZRT.csv", sep = ';', header=T, na.string="NA");
tt$rok <- substr(tt$month, 1, 4)
tt$yyyymmdd <- as.Date(sprintf ("%s-%s", substr(tt$month, 6, 7), substr(tt$month, 9, 10) ), format="%m-%d")
tt$diff <- tt$urodzenia - tt$zgony
tt.yy2020 <- tt %>% filter(rok==2020) %>% as.data.frame
tt.yy2019 <- tt %>% filter(rok==2019) %>% as.data.frame
tt.yy2018 <- tt %>% filter(rok==2018) %>% as.data.frame
tt.yy2017 <- tt %>% filter(rok==2017) %>% as.data.frame
tt.yy2016 <- tt %>% filter(rok==2016) %>% as.data.frame
tt.yy2015 <- tt %>% filter(rok==2015) %>% as.data.frame
## max zgony
tt.yy2020$max.so.far <- pmax (tt.yy2019$zgony, tt.yy2018$zgony,
tt.yy2017$zgony, tt.yy2016$zgony, tt.yy2015$zgony)
## średnia
tt.yy2020$mean <- (tt.yy2019$zgony + tt.yy2018$zgony
+ tt.yy2017$zgony + tt.yy2016$zgony + tt.yy2015$zgony)/5
## odchylenie std:
tt.yy2020$sqmean <- (tt.yy2019$zgony^2 + tt.yy2018$zgony^2 +
tt.yy2017$zgony^2 + tt.yy2016$zgony^2 + tt.yy2015$zgony^2)/5
tt.yy2020$sq <- sqrt(tt.yy2020$sqmean - tt.yy2020$mean^2)
tt.yy2020$diff20 <- tt.yy2020$zgony - tt.yy2020$max.so.far
sprintf ("%s;%.2f;%.2f;%.2f;%.2f",
tt.yy2020$yyyymmdd, tt.yy2020$diff20, tt.yy2020$sqmean, tt.yy2020$sq, tt.yy2020$max.so.far )
MZ to Ministerstwo Zdrowia. Specyfiką PL jest brak danych publicznych nt. Pandemii.. Na stronie GIS nie ma nic, a stacje wojewódzkie publikują jak chcą i co chcą. Ministerstwo zdrowia z kolei na swojej stronie podaje tylko dane na bieżący dzień, zaś na amerykańskim Twitterze publikuje komunikaty na temat. To co jest na stronie znika oczywiście następnego dnia zastępowane przez inne cyferki. Do tego są tam tylko dzienne dane zbiorcze o liczbie zarażonych i zmarłych, (w podziale na województwa). Nie ma na przykład wieku, płci itp... To co jest na Twitterze z kolei ma formę tekstowego komunikatu postaci: Mamy 502 nowe i potwierdzone przypadki zakażenia #koronawirus z województw: małopolskiego (79), śląskiego (77), mazowieckiego (75)... [...] Z przykrością informujemy o śmierci 6 osób zakażonych koronawirusem (wiek-płeć, miejsce zgonu): 73-K Kędzierzyn-Koźle, 75-M Łańcut, 92-K Lipie, 72-M, 87-M, 85-K Gdańsk.
Czyli podają zarażonych ogółem i w podziale na województwa oraz dane indywidualne zmarłych w postaci płci i wieku oraz miejsca zgonu (miasta żeby było inaczej niż w przypadku podawania zakażeń.)
No więc chcę wydłubać dane postaci 92-K z tweetów publikowanych
przez Ministerstwo Zdrowia. W tym celu za pomocą
tweepy pobieram cały streamline (+3200 tweetów zaczynających się jeszcze w 2019 roku),
ale dalej to już działam w Perlu, bo w Pythonie jakby mniej komfortowo się czuję.
Po pierwsze zamieniam format json na csv:
use JSON;
use Data::Dumper;
use Time::Piece;
use open ":encoding(utf8)";
use open IN => ":encoding(utf8)", OUT => ":utf8";
binmode(STDOUT, ":utf8");
## ID = tweeta ID; date = data;
## repid -- odpowiedź na tweeta o numerze ID
## text -- tekst tweeta
print "id;date;repid;text\n";
while (<>) {
chomp();
$tweet = $_;
my $json = decode_json( $tweet );
#print Dumper($json);
$tid = $json->{"id"};
$dat = $json->{"created_at"};
## Data jest w formacie rozwlekłym zamieniamy na YYYY-MM-DDTHH:MM:SS
## Fri Oct 04 14:48:25 +0000 2019
$dat = Time::Piece->strptime($dat,
"%a %b %d %H:%M:%S %z %Y")->strftime("%Y-%m-%dT%H:%M:%S");
$rep = $json->{"in_reply_to_status_id"};
$ttx = $json->{"full_text"}; $ttx =~ s/\n/ /g;
## Zamieniamy ; na , w tekście żeby użyć ; jako separatora
$ttx =~ s/;/,/g; ####
print "$tid;$dat;$rep;$ttx\n";
Komunikaty dłuższe niż limit Twittera są dzielone na kawałki, z których każdy jest odpowiedzią na poprzedni, np:
1298900644267544576;2020-08-27T08:30:20;1298900642522685440;53-M, 78-K i 84-K Kraków. Większość osób ... 1298900642522685440;2020-08-27T08:30:20;1298900640714948608;67-K Lublin (mieszkanka woj. podkarpackiego), 85-K Łańcut,... 1298900640714948608;2020-08-27T08:30:19;1298900639586680833;kujawsko-pomorskiego (24), świętokrzyskiego (18), opolskiego... 1298900639586680833;2020-08-27T08:30:19;;Mamy 887 nowych i potwierdzonych przypadków zakażenia #koronawirus z województw: ...
Czyli tweet 1298900639586680833 zaczyna, 1298900640714948608 jest odpowiedzią na 1298900639586680833,
a 1298900642522685440 odpowiedzią na 1298900640714948608 itd. Tak to sobie wymyślili...
W sumie chyba niepotrzebnie, ale w pierwszym kroku agreguję podzielone komunikaty
w ten sposób, że wszystkie odpowiedzi są dołączane
do pierwszego tweeta (tego z pustym polem in_reply_to_status_id):
## nextRef jest rekurencyjna zwraca numer-tweeta,
## który jest początkiem wątku
sub nextRef {
my $i = shift;
if ( $RR{"$i"} > 0 ) {
return ( nextRef( "$RR{$i}" ) );
} else { return "$i" }
}
### ### ###
while (<>) { chomp();
($id, $d, $r, $t) = split /;/, $_;
$TT{$id} = $t;
$RR{$id} = $r;
$DD{$id} = $d;
}
### ### ###
for $id ( sort keys %TT ) {
$lastId = nextRef("$id");
$LL{"$id"} = $lastId;
$LLIds{"$lastId"} = "$lastId";
}
### ### ###
for $id (sort keys %TT) {
## print "### $DD{$id};$id;$LL{$id};$TT{$id}\n"; }
$TTX{$LL{"$id"}} .= " " . $TT{"$id"};
$DDX{$LL{"$id"}} .= ";" . $DD{"$id"};
}
### ### ###
for $i (sort keys %TTX) {
$dates = $DDX{"$i"};
$dates =~ s/^;//; ## pierwszy ; jest nadmiarowy
@tmpDat = split /;/, $dates;
$dat_time_ = $tmpDat[0];
($dat_, $time_) = split /T/, $dat_time_;
$ffN = $#tmpDat + 1;
$collapsedTweet = $TTX{$i};
print "$i;$dat_;$time_;$ffN;$collapsedTweet\n";
}
Zapuszczenie powyższego powoduje konsolidację wątków, tj. np. powyższe 4 tweety z 2020-08-27 połączyły się w jeden:
1298900639586680833;2020-08-27;08:30:19;4; Mamy 887 nowych i potwierdzonych przypadków zakażenia #koronawirus z województw: małopolskiego (233), śląskiego (118), mazowieckiego (107), [...] Liczba zakażonych koronawirusem: 64 689 /2 010 (wszystkie pozytywne przypadki/w tym osoby zmarłe).
Teraz wydłubujemy tylko tweety z frazą nowych i potwierdzonych albo nowe i potwierdzone:
## nowych i potwierdzonych albo nowe i potwierdzone ## (MZ_09.csv to CSV ze `skonsolidowanymi' wątkami) cat MZ_09.csv | grep 'nowych i pot\|nowe i pot' > MZ_09_C19.csv wc -l MZ_09_C19.csv 189 MZ_09_C19.csv
Wydłubanie fraz wiek-płeć ze skonsolidowanych wątków jest teraz proste:
perl -e '
while (<>) {
($i, $d, $t, $c, $t) = split /;/, $_;
while ($t =~ m/([0-9]+-[MK])/g ) {
($w, $p) = split /\-/, $1;
print "$d;$w;$p\n";
}
}' MZ_09_C19.csv > C19D.csv
wc -l C19PL_down.csv
1738
Plik wykazał 1738 osób. Pierwsza komunikat jest z 16 kwietnia. Ostatni z 31. sierpnia. Pierwsze zarejestrowane zgony w PL odnotowano 12 marca (albo 13 nieważne). W okresie 12 marca --15 kwietnia zmarło 286 osób. Dodając 286 do 1738 wychodzi 2024. Wg MZ w okresie 12.03--31.08 zmarło 2039. Czyli manko wielkości 15 zgonów (około 0,5%). No trudno, nie chce mi się dociekać kto i kiedy pogubił tych 15...
Równie prostym skryptem zamieniam dane indywidualne na tygodniowe
#!/usr/bin/perl -w
use Date::Calc qw(Week_Number);
while (<>) {
chomp();
##if ($_ =~ /age/) { next } ##
my ($d, $w, $p ) = split /;/, $_;
my ($y_, $m_, $d_) = split /\-/, $d;
my $week = Week_Number($y_, $m_, $d_);
$DW{$week} += $w;
$DN{$week}++;
$DD{$d} = $week;
$YY{$week} = 0;
$YY_last_day{$week} = $d;
## wiek wg płci
$PW{$p} += $w;
$PN{$p}++;
}
for $d (keys %DD) {
$YY{"$DD{$d}"}++; ## ile dni w tygodniu
}
print STDERR "Wg płci/wieku (ogółem)\n";
for $p (keys %PW) {
$s = $PW{$p}/$PN{$p};
printf STDERR "%s %.2f %i\n", $p, $s, $PN{$p};
$total += $PN{$p};
}
print STDERR "Razem: $total\n";
print "week;deaths;age;days;date\n";
for $d (sort keys %DW) {
if ($YY{$d} > 2 ) {## co najmniej 3 dni
$s = $DW{$d}/$DN{$d};
printf "%s;%i;%.2f;%i;%s\n", $d, $DN{$d}, $s, $YY{$d}, $YY_last_day{$d};
}
}
Co daje ostatecznie (week -- numer tygodnia w roku;
meanage -- średni wiek zmarłych;
deaths -- liczba zmarłych w tygodniu;
days -- dni w tygodniu;
date -- ostatni dzień tygodnia):
week;deaths;meanage;days;date 16;55;77.07;4;2020-04-19 17;172;75.09;7;2020-04-26 18;144;77.29;7;2020-05-03 19;123;76.46;7;2020-05-10 20;126;76.40;7;2020-05-17 21;71;76.37;7;2020-05-24 22;68;78.12;7;2020-05-31 23;93;75.73;7;2020-06-07 24;91;75.93;7;2020-06-14 25;109;77.24;7;2020-06-21 26;83;75.06;7;2020-06-28 27;77;74.09;7;2020-07-05 28;55;76.91;7;2020-07-12 29;54;77.33;7;2020-07-19 30;48;76.52;7;2020-07-26 31;60;74.88;7;2020-08-02 32;76;77.17;7;2020-08-09 33;71;73.11;7;2020-08-16 34;77;75.61;7;2020-08-23 35;79;74.33;7;2020-08-30
To już można na wykresie przedstawić:-)
library("dplyr")
library("ggplot2")
library("scales")
##
spanV <- 0.25
d <- read.csv("C19D_weekly.csv", sep = ';', header=T, na.string="NA")
first <- first(d$date)
last <- last(d$date)
period <- sprintf ("%s--%s", first, last)
d$deaths.dailymean <- d$deaths/d$days
cases <- sum(d$deaths);
max.cases <- max(d$deaths)
note <- sprintf ("N: %i (source: twitter.com/MZ_GOV_PL)", cases)
pf <- ggplot(d, aes(x= as.Date(date), y=meanage)) +
geom_bar(position="dodge", stat="identity", fill="steelblue") +
scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50,55,60,65,70,75,80)) +
xlab(label="week") +
ylab(label="mean age") +
ggtitle(sprintf ("Mean age of COVID19 fatalities/Poland/%s", period), subtitle=note )
note <- sprintf ("N: %i (source: twitter.com/MZ_GOV_PL)", cases)
pg <- ggplot(d, aes(x= as.Date(date), y=deaths.dailymean)) +
geom_bar(position="dodge", stat="identity", fill="steelblue") +
scale_x_date( labels = date_format("%m/%d"), breaks = "2 weeks") +
scale_y_continuous(breaks=c(0,5,10,15,20,25,30,35,40,45,50)) +
xlab(label="week") +
ylab(label="daily mean") +
ggtitle(sprintf("Daily mean number of COVID19 fatalities/Poland/%s", period), subtitle=note )
ggsave(plot=pf, file="c19dMA_w.png")
ggsave(plot=pg, file="c19dN_w.png")
Wynik tutaj:
Średnia Wg płci/wieku/ogółem. Ogółem -- 75,9 lat (1738 zgonów) w tym: mężczyźni (914 zgonów) -- 73.9 lat oraz kobiety (824) -- 78.4 lat. Stan pandemii na 31.08 przypominam. Apogeum 2 fali... Średnia wieku w PL (2019 rok) 74,1 lat (M) oraz 81,8 lat (K).