Když jsem před lety sepisoval populárně poučný příspěvek o odhadování exponenciály metodou nejmenších čtverců – dělá se to přes stats::nls() – tak jsem použil vědomě lehčí formulace problému coby mezd ajťáků v časech Andreje Babiše. Tím jsem svojí exponenciálu zasadil do času a prostoru.
Rok se sešel s rokem, volební období s obdobím, a z obrazovek nás opět bombardují politici s tím že tentokráte určitě bude líp, stačí když dáme hlas jim (a ne těm druhým). Přišlo mi to jako zajímavá příležitost k tomu, abych se vrátil do staré řeky a oprášil techniky práce s geometrickou řadou v kontextu statistického programovacího jazyka R.
Cílem cvičení bude:
- připomenout techniku, jak pomocí balíčku
{czso}přistupovat k datům z Veřejné databáze Českého statistického úřadu přímo do R vaší session - načtený dataset zpracovat technikami světa
{tidyverse}, konkrétně balíčky{dplyr}pro datovou manipulaci a{ggplot2}pro statickou vizualizaci - o datasetu mezd a jejich vývoji v čase podat zprávu graficky
- nad datasetem mezd sestavit matematický model, a uplatnit jej
Protože mzdy a jejich vývoj v čase je složitý problém, který není snadné plně podchytit v jeho komplexnosti, zaměřím se na dva sektory ve kterých se pohybuju a se kterými mám žitou zkušenost:
- jako konzultant působím v sektoru
J - Informační a komunikační činnosti(“ajťáci”) - jako student a vyučující na VŠE patřím do sektoru
P - Vzdělávání(“učitelky”)
Základem je načtení dat:
library(czso) # protože staťák...
library(tidyverse) # kvůli dplyr & ggplot2
library(zoo) # pro konverzi datumů / kvartály v letech
# načtu dataset s id 110079 = Zaměstnanci a průměrné hrubé měsíční mzdy podle odvětví
raw_mzdy <- czso::czso_get_table("110079")
# odliji stranou čistší dataset
clean_mzdy <- raw_mzdy %>%
filter(stapro_kod == '5958' # průměrná mzda / zahazuju počty zaměstnanců
& rok >= 2016
& typosoby_kod == '200'# přepočet na ekvivalent plného úvazku
& odvetvi_kod %in% c('P', 'J')) %>% # vybrané sektory
# konverze konvence z roku + čtvrtletí na prosté datum (první den kvartálu)
mutate(datum = as.Date(as.yearqtr(paste0(rok, "Q", ctvrtleti)))) %>%
arrange(datum) %>% # setřídím dle data / pro jistotu, kvůli sekvenci v příštím kroku
group_by(odvetvi_kod) %>%
mutate(sekvence = rank(datum),
posledni = datum == max(datum)) %>% # poslední záznam - pro popisku grafů
ungroup()Nad načtenými daty připravím vizualizaci technikami ggplot2-u; je zajímavé pozorovat jak mzdy ajťáků i učitelek mají stejně jasný sezónní trend (1× ročně bonusy) ale se špičkami v různých čtvrtletích.
Pracovně se domnívám, že učitelky mají odměny spíše svázané s rozpočtem (vyplácené v Q4, aby se utratilo), kdežto ajťáci spíše s hospodářským výsledkem (vyplácené v Q1, po závěrce roku).
# základní overview graficky
ggplot(data = clean_mzdy, aes(x = datum, y = hodnota,
fill = odvetvi_txt)) +
geom_point(pch = 21, stroke = NA) +
scale_y_continuous(limits = c(0, NA),
labels = scales::dollar_format(prefix = "",
suffix = " Kč",
big.mark = " ")) +
scale_x_date(breaks = seq(as.Date("2015-01-01"),
as.Date("2025-01-01"),
by = "2 years"),
date_minor_breaks = "1 year",
date_labels = "%Y") +
labs(title = "Průměrná hrubá měsíční mzda v čase",
fill = "Sektor podle ČSÚ") +
theme_minimal() +
theme(legend.position = c(0.795, 0.175),
legend.direction = "vertical",
axis.title = element_blank())
Kromě výrazné roční periody je v grafech (hlavně na učitelkách) intuitivně cítit změna trendu kolem konce roku 2021. Tuto intuici můžeme ověřit odděleným natrénováním matematického modelu nad érami Andreje Babiše a Petra Fialy.
Předpoklad za modelem je, že se mzda se chová jako exponenciála (každý kvartál vzroste o stejné relativní procento) a ne jako přímka (každý kvartál vzroste o stejnou absolutní částku).
# omezený dataset "éra Andreje Babiše" pro trénink modelu
era_babise <- clean_mzdy %>%
filter(datum >= as.Date("2017-12-13") # počátek první vlády AB
& datum < as.Date("2021-12-17")) # konec druhé vlády AB
# omezený dataset "éra Petra Fialy" pro trénink modelu
era_fialy <- clean_mzdy %>%
filter(datum >= as.Date("2021-12-17")) # jmenování PF prezidentem MZ
# 4x matematický model přes stats::nls()
model_babis_p <- nls(hodnota ~ a * (1 + r)^sekvence,
data = subset(era_babise, odvetvi_kod == "P"),
start = list(a = 50000, r = .01))
model_babis_j <- nls(hodnota ~ a * (1 + r)^sekvence,
data = subset(era_babise, odvetvi_kod == "J"),
start = list(a = 50000, r = .01))
model_fiala_p <- nls(hodnota ~ a * (1 + r)^sekvence,
data = subset(era_fialy, odvetvi_kod == "P"),
start = list(a = 50000, r = .01))
model_fiala_j <- nls(hodnota ~ a * (1 + r)^sekvence,
data = subset(era_fialy, odvetvi_kod == "J"),
start = list(a = 50000, r = .01)) Když přepočteme kvartální růst mezd na roční ekvivalenty, tak se dostaneme na následující čísla:
- ajťákům v časech Andreje Babiše rostly mzdy o
4.74%ročně - ajťákům v časech Petra Fialy rostly mzdy o
8.91%ročně - učitelkám v časech Andreje Babiše rostly mzdy o
11.2%ročně - učitelkám v časech Petra Fialy rostly mzdy o
4.52%ročně
Je velká otázka, jak moc můžeme modelovaný nárůst (a jeho změnu mezi vládami) přisuzovat za zásluhu konkrétnímu premiérovi. Těžko tvrdit, že Babiš může za Covid, případně Fiala za Ukrajinu. Ale říct, že za časů jednoho bylo jednak, a za časů druhého druhak, jde snadno - datumy začátku a konce vlád jsou ložené.
A také si můžeme společně zaspekulovat: jak by vypadal svět, pokud by se trend nezměnil? Jak by vypadaly dnes mzdy ajťáků a učitelek, pokud by pokračovalo tempo předchozí vlády?
Pomocníkem v našem spekulování bude metoda stats::predict.nls(), jejímž uplatněním na natrénovaný matematický model získáme nová data – a s jejich pomocí obohatíme původní obrázek o původní a nový trend.
Pro snazší orientaci barvím éru Petra Fialy fialově, a éru Andreje Babiše korporátní zelenou.
# pomocný dataset pro oddělený výpočet trendů v érách
trendy <- clean_mzdy %>%
select(datum, sekvence, posledni) %>%
unique() %>%
mutate(era = case_when(datum <= as.Date("2017-12-13") ~ "sobotka",
datum <= as.Date("2021-12-17") ~ "babis",
T ~ "fiala"))
# trendy ajťáků v datasetu
trendy$ab_p <- predict(model_babis_p, newdata = trendy)
trendy$pf_p <- predict(model_fiala_p, newdata = trendy)
# trendy učitelů v datasetu
trendy$ab_j <- predict(model_babis_j, newdata = trendy)
trendy$pf_j <- predict(model_fiala_j, newdata = trendy)
# výsledný obrázek - původní ggplot, doplněný o trendové čáry
ggplot(data = clean_mzdy, aes(x = datum, y = hodnota)) +
geom_point(aes(fill = odvetvi_txt),
stroke = NA,
pch = 21) +
scale_y_continuous(limits = c(0, NA),
labels = scales::dollar_format(prefix = "",
suffix = " Kč",
big.mark = " ")) +
scale_x_date(breaks = seq(as.Date("2015-01-01"),
as.Date("2025-01-01"),
by = "2 years"),
date_minor_breaks = "1 year",
date_labels = "%Y") +
geom_line(data = trendy, aes(x = datum, y = ab_p),
color = "grey", lty = "dashed") +
geom_line(data = subset(trendy, era == "babis"),
aes(x = datum, y = ab_p, color = "AB")) +
geom_line(data = subset(trendy, era == "fiala"),
aes(x = datum, y = pf_p, color = "PF")) +
geom_line(data = trendy, aes(x = datum, y = ab_j),
color = "grey", lty = "dashed") +
geom_line(data = subset(trendy, era == "babis"),
aes(x = datum, y = ab_j, color = "AB")) +
geom_line(data = subset(trendy, era == "fiala"),
aes(x = datum, y = pf_j, color = "PF")) +
geom_text(data = subset(clean_mzdy, posledni),
aes(x = datum, y = hodnota, label = scales::comma(hodnota)),
nudge_x = 250) +
geom_text(data = subset(trendy, posledni),
aes(x = datum, y = ab_p, label = scales::comma(ab_p)),
color = "gray70",
nudge_x = 250) +
geom_text(data = subset(trendy, posledni),
aes(x = datum, y = ab_j, label = scales::comma(ab_j)),
color = "gray70",
nudge_x = 250) +
scale_color_manual("Éra vlády",
values = c("AB" = "#5c9234",
"PF" = "darkorchid")) +
labs(title = "Trendy mezd v érách Babiše a Fialy",
fill = "Sektor podle ČSÚ") +
theme_minimal() +
guides(color = guide_legend(position = "bottom",
direction = "horizontal")) +
theme(legend.position = c(0.795, 0.175),
legend.direction = "vertical",
axis.title = element_blank())
Z grafu je patrné, že změna trendu měla na ajťáky a učitelky opačný vliv:
- ajťáci se dnes proti trendu z časů Andreje Babiše pohybují přibližně +15 tisíc Kč / měsíc
- učitelky se dnes proti trendu z časů Andreje Babiše pohybují přibližně -20 tisíc Kč / měsíc
A protože jsme v Česku, kde se nehraje až tolik na to abych se měl absolutně dobře (co je absolutně dobře? jaká je jednotka dobrosti?) ale jestli se mám relativně líp jak soused, tak se podíváme na relativní poměr mezd ajťáků a učitelek:
rel_mzdy <- clean_mzdy %>%
# pivot z "dlouhého" formátu na "široký"
pivot_wider(id_cols = datum,
names_from = odvetvi_kod,
values_from = hodnota) %>%
mutate(pomer = P / J) # podíl učitelek ku ajťákům
ggplot(data = rel_mzdy, aes(x = datum, y = pomer)) +
annotate("rect",
xmin = as.Date("2017-12-13"), # počátek první vlády AB
xmax = as.Date("2021-12-17"), # konec druhé vlády AB
ymin = 0,
ymax = Inf,
fill = "#5c9234",
alpha = 1/5) +
annotate("rect",
xmin = as.Date("2021-12-17"), # jmenování PF premiérem
xmax = max(clean_mzdy$datum), # konec datová řady
ymin = 0,
ymax = Inf,
fill = "darkorchid",
alpha = 1/5) +
geom_smooth(se = F, color = "red", alpha = 2/3) +
geom_point(pch = 4, size = 3/4, color = "gray25") +
geom_point(data = rel_mzdy[c(which.max(rel_mzdy$pomer),
which.min(rel_mzdy$pomer)),],
color = "red") +
geom_text(data = rel_mzdy[which.max(rel_mzdy$pomer),],
aes(x = datum, y = pomer,
label = paste0(round(100 * pomer, 2), "%")),
nudge_y = .04) +
geom_text(data = rel_mzdy[which.min(rel_mzdy$pomer),],
aes(x = datum, y = pomer,
label = paste0(round(100 * pomer, 2), "%")),
nudge_y = -.04) +
scale_y_continuous(limits = c(0, NA),
labels = scales::percent) +
scale_x_date(breaks = seq(as.Date("2015-01-01"),
as.Date("2025-01-01"),
by = "2 years"),
date_minor_breaks = "1 year",
date_labels = "%Y") +
theme_minimal() +
theme(axis.title = element_blank()) +
labs(title = "Relativní srovnání mezd sektorů Vzdělávání a ICT",
subtitle = "v érách premiérů Babiše a Fialy ")
Graf je zrnitý (pamatujete, že učitelky mívají bonusy v jiném kvartálu než ajťáci?) ale dá se technikami erka vyhladit. A dvě informace z něj vyskakují zřetelně:
- historicky nejvyšší poměr mezd průměrné učitelky k průměrnému ajťákovi byl 2021-Q4 (tedy poslední kvartál éry Andreje Babiše)
- historicky nejmenší poměr mezd průměrné učitelky k průměrnému ajťákovi byl 2025-Q1 (tedy právě teď)
Díky matematickému modelování a technikám statistického programovacího jazyka R tak vidím, že změna poměrů z časů vlády Andreje Babiše k Petru Fialovi měla – alespoň co se mzdovéhovo vývoje týče – kromě vítězů také poražené. Učitelky, které se ptají “kde je moje dvacka?!” přitom lidsky chápu; stejně tak jako ajťáky kteří mlčí a šoupou nohama, že oni nic – to trh.