Porozumění trhu práce je dobrou příležitostí k uplatnění statistického programovacího jazyka R.

Erko nám umožňuje několik věcí:

V rámci modelu budu uvažovat geometrickou posloupnost – růst mezd o stálé procento. Tento vztah odpovídá ekonomické teorii lépe, než závislost lineární (růst mezd o stálou částku).

Prvním krokem je načtení knihoven a akvizice surových dat. Základním vstupem pro analýzu bude vstupem standardní datová sada číslo 110079 – Zaměstnanci a průměrné hrubé měsíční mzdy podle odvětví. Dataset má kvartální periodu.

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")

Když jsme datovou sadu načetli, tak se na hrubo seznámíme s její strukturou a obsahem:

glimpse(raw_mzdy)
## Rows: 6,880
## Columns: 16
## $ idhod        <chr> "741383707", "741383708", "741383709", "741383713", "7413…
## $ hodnota      <dbl> 23546, 24057, 27242, 22691, 24135, 24635, 27830, 10640, 1…
## $ stapro_kod   <chr> "5958", "5958", "5958", "5958", "5958", "5958", "5958", "…
## $ mj_cis       <chr> "78", "78", "78", "78", "78", "78", "78", "78", "78", "78…
## $ mj_kod       <chr> "00200", "00200", "00200", "00200", "00200", "00200", "00…
## $ typosoby_kod <chr> "200", "200", "200", "200", "200", "200", "200", "200", "…
## $ odvetvi_cis  <chr> "5103", "5103", "5103", "5103", "5103", "5103", "5103", "…
## $ odvetvi_kod  <chr> "P", "P", "P", "P", "P", "P", "P", "Q", "Q", "H", "H", "H…
## $ rok          <int> 2012, 2012, 2012, 2013, 2013, 2013, 2013, 2000, 2000, 200…
## $ ctvrtletí    <chr> "2", "3", "4", "1", "2", "3", "4", "1", "2", "1", "2", "3…
## $ uzemi_cis    <chr> "97", "97", "97", "97", "97", "97", "97", "97", "97", "97…
## $ uzemi_kod    <chr> "19", "19", "19", "19", "19", "19", "19", "19", "19", "19…
## $ stapro_txt   <chr> "Průměrná hrubá mzda na zaměstnance", "Průměrná hrubá mzd…
## $ mj_txt       <chr> "Kč", "Kč", "Kč", "Kč", "Kč", "Kč", "Kč", "Kč", "Kč", "Kč…
## $ typosoby_txt <chr> "přepočtený", "přepočtený", "přepočtený", "přepočtený", "…
## $ odvetvi_txt  <chr> "Vzdělávání", "Vzdělávání", "Vzdělávání", "Vzdělávání", "…

Dataset obsahuje v normalizované formě více statistických veličin (mzdu a počet zaměstnanců) ve dvou dimenzích (surová data, a jejich přepočet na plné úvazky / FTEs = Full Time Employees). Pro eliminaci duplicit bude před vlastním modelováním nutno nastavit filtr.

Mzdy jsou evidované po odvětvích; pro zjednodušení se v úvodní vizualizaci zaměříme na tři, u kterých je větší pravděpodobnost využití technik jazyka R: ICT pracovníky, vědecké pracovníky a učitele a učitelky.

# odliji stranou čistší dataset
clean_mzdy <- raw_mzdy %>% 
  filter(stapro_kod == '5958'  # průměrná mzda / zahazuju počty zaměstnanců
         & typosoby_kod == '200'# přepočet na ekvivalent plného úvazku
         & odvetvi_kod %in% c('P', 'M', 'J')) %>% # tři vybrané sektory
  # konverze konvence z roku + čtvrtletí na prosté datum (první den kvartálu)
  mutate(datum = as.Date(as.yearqtr(paste0(rok, "Q", ctvrtletí))))

# základní overview graficky
ggplot(data = clean_mzdy, aes(x = datum, y = hodnota, color = odvetvi_txt)) +
  geom_point() +
  scale_y_continuous(limits = c(0, 70000), 
                     labels = scales::dollar_format(prefix = "", 
                                                    suffix = " Kč",
                                                    big.mark = " ")) +  
  labs(title = "Průměrná hrubá měsíční mzda v čase",
       color = "Sektor podle ČSÚ") +
  theme_minimal() +
  theme(axis.title = element_blank(),
        plot.title = element_text(hjust = 1/2, 
                                  size = 14),
        legend.position = c(0.795, 0.175),
        legend.direction = "vertical",
        legend.title = element_text(hjust = 1/2),
        legend.background = element_rect(fill = "white",
                                         color = NA))

Z grafu můžeme vypozorovat tři fáze cyklu mezd:

  • růst mezi rokem 2000 (počátek známé historie) a rokem 2010
  • spíše stagnaci mezi lety 2010 a 2015
  • opětovný růst od roku 2015 dále

Pro vlastní model zvolíme období růstu po roce 2015, kterému můžeme říkat Éra Andreje Babiše – a podle toho, jak moc panu premiérovi fandíme, budeme uvažovat, že to byl právě on, kdo:

  • konečně zařídil peníze pro naše lidi (je to pašák!)
  • nezodpovědně roztočil inflační spirálu (je to neřád!)

Na další práci s modelem už náš názor na pana premiéra vliv mít nebude.

Pro model vybereme jeden obor, a sice ICT pracovníky neboli ajťáky.

# připravíme podkladový dataset "ajťáci v časech Andreje Babiše" pro trénink modelu
era_babise <- clean_mzdy %>% 
  filter(datum >= as.Date("2015-01-01") 
         & odvetvi_kod == 'J') %>% # pouze ICT sektor
  arrange(datum) %>% # setřídím dle data / pro jistotu, kvůli sekvenci v příštím kroku
  mutate(sekvence = 1:n()) # pořadové číslo kvartálu v rámci Éry A. B.

Pro vlastní matematický model budeme uvažovat geometrickou posloupnost - jinými slovy budeme předpokládat, že mzdy ajťáků vzrostou z kvartálu na kvartál pokaždé o stejné procento.

Taková posloupnost se dá popsat vzorcem a * (1 + r) ^ sekvence, kde a je počáteční hodnota (mzda v IT na začátku éry Andreje Babiše), r je míra růstu z kvartálu na kvartál a sekvence je pořadové číslo kvartálu v rámci éry Andreje Babiše.

Pro nalezení konkrétních hodnot parametrů a a r použijeme techniku nejmenších čtverců, konkrétně funkci stats::nls(). Povinnými vstupy jsou podkladová data a vzorec očekávané závislosti v tildové notaci. Pro rychlejší konvergenci modelu můžeme doplnit iniciální hodnoty parametrů.

# natrénujeme matematický model přes stats::nls()
matematicky_model <- nls(hodnota ~ a * (1 + r)^sekvence, # vzorec, v tildové notaci
                         data = era_babise, # podkladová data
                         # vstupní odkad parametrů / odhad přes palec
                         start = list(a = 50000, 
                                      r = .01)) 

# shrnutí modelu - hodnoty parametrů, významnost &c.
summary(matematicky_model)
## 
## Formula: hodnota ~ a * (1 + r)^sekvence
## 
## Parameters:
##    Estimate Std. Error t value Pr(>|t|)    
## a 4.703e+04  5.795e+02   81.14  < 2e-16 ***
## r 1.231e-02  7.507e-04   16.40 1.53e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1571 on 24 degrees of freedom
## 
## Number of iterations to convergence: 3 
## Achieved convergence tolerance: 3.634e-08

Model je jednoduchý a teoreticky dobře podložený, takže nás nepřekvapí že pro oba hledané parametry vyšly odhady jako statisticky významné.

  • modelová mzda v IT sektoru na počátku éry nám vyšla 47 025 Kč; pro kontext: skutečnost Q4 2014 byla 48 744 Kč, což představuje rozdíl 3.5%.
  • modelové tempo kvartálního růstu hrubých mezd v ICT nám vyšlo 1.23%.

Získaný model můžeme snadno protáhnout do budoucna, a porovnat graficky se skutečností:

# podklad pro uplatnění matematického modelu
jiny_cas <- data.frame(sekvence = 0:30) %>% 
  mutate(datum = as.Date("2015-01-01") + months(3 * sekvence))

# predikce - uplatnění matematického modelu na nový vstup
jiny_cas$hodnota <- predict(matematicky_model,
                            newdata = jiny_cas)


# grafické podání zprávy o skutečnosti a modelu
ggplot() +
  geom_line(data = jiny_cas, aes(x = datum, 
                                 y = hodnota), 
            color = "red") +
  geom_point(data = era_babise, aes(x = datum, 
                                   y = hodnota),
             pch = 4, color = "gray50") +
  scale_y_continuous(limits = c(40, 70) * 1000, # menší rozsah osy, ať vyniknou epsilonky
                     labels = scales::dollar_format(prefix = "", 
                                                    suffix = " Kč",
                                                    big.mark = " ")) +
  labs(title = "Model a skutečnost průměrné mzdy v českém ICT sektoru v čase") +
  theme_minimal() +
  theme(axis.title = element_blank(),
        plot.title = element_text(hjust = 1/2, 
                                  size = 14))

Vidíme, že model aproximuje skutečnost poměrně přesně; největší rozdíly jsou v prvních kvartálech roku, kdy je výrazná sezónnost (v reálném světě jí táhnou v tomto čtvrtletí vyplácené roční odměny, což je vlastnost kterou v rámci modelu neuvažujeme).

Když jsme model připravili, tak je čas ho použít. Obecně se nabízí dvě možnosti uplatnění:

  • předpověď průměrné mzdy v budoucnu (po anglicku prediction)
  • porozumění obecné problematice růstu mezd (po anglicku inference)

První bod bude relevantní jen pro někoho – ideálně průměrná mzda je umělý konstrukt, za který nikdo reálně nepracuje – ale druhý bod je zajímavější.

Čeští ajťáci si mohou snadno porovnat realitu změny svého platu s trhem jako celkem. Pokud rostou pod trhem, tak svoji relativní pozici v poli všech ICT zaměstnanců ztrácí; pokud rostou nad trhem, tak se v pelotonu posunují dopředu.

České konvence pracovního trhu nepočítají s kvartální úpravou mezd, obvyklejší je revize roční; očekávanou výši ročního nárůstu ze získané hodnoty koeficientu r získáme snadno.

# očekávaný roční nárůst mezd z titulu inflace
(1 + coef(matematicky_model)[["r"]]) ^ 4 - 1
## [1] 0.05016941

Jinými slovy pro ideálně průměrnou IT pozici očekáváme každoroční navýšení o 5.02% čistě z titulu běhu času a posunu trhu jako celku - za to, že poběžíme na místě s Červenou královnou.

V některých kontextech není obvyklé automatické navyšování, a o lepší peníze si musí člověk říct. Což s sebou nese určité tření, a nejde to dělat příliš často – je tedy vhodnější méně časté navýšení o větší částku.

V takovém případě může být zajímavé vědět, kolik času zabere trhu kumulativní nárůst o určitou hodnotu – řekněme 10%. Odpověď na tuto otázku nám dá opět koeficient r našeho modelu, jen musíme použít logaritmy.

# čas v kvartálech pro nárůst tržní mzdy o 10%
log(1 + 10/100) / log(1 + coef(matematicky_model)[["r"]])
## [1] 7.788132

Jinými slovy ideálně průměrný ajťák může očekávat nárůst mzdy o 10% někdy mezi 7. a 8. kvartálem od počátku.

Příklad vývoje mezd v IT byl vědomě lehčího rázu; nicméně věřím, že jsem v jeho rámci předvedl jak lehkost, s jakou lze přistupovat k datům ČSÚ přímo z pohodlí vaší R Session, tak eleganci erkových modelovacích nástrojů a přesvědčivost grafických výstupů. QED.