Využívanie niektorých knižníc

rm(list=ls())
library(lmtest)   #  podpora regresie
library(outliers) # analyza odlahlych hodnot (outliers)
library(gptstudio)
library(kableExtra)
library(knitr)
library(dplyr)
library(broom)
library(corrplot)

Príprava údajov - import z csv súboru

Súbor Life_Expectancy_Data obsahuje databázu determinantov očakávanej dĺžky života. Import údajov urobíme nasledovne

# import the dataset and create a data.frame udaje
#
udaje_svet <- read.csv("udaje_svet/Life-Expectancy-Data-Updated.csv",header=TRUE,sep=",",dec=".",check.names = TRUE)
head(udaje_svet)
##      Country                        Region Year Infant_deaths Under_five_deaths
## 1    Turkiye                   Middle East 2015          11.1              13.0
## 2      Spain                European Union 2015           2.7               3.3
## 3      India                          Asia 2007          51.5              67.9
## 4     Guyana                 South America 2006          32.8              40.5
## 5     Israel                   Middle East 2012           3.4               4.3
## 6 Costa Rica Central America and Caribbean 2006           9.8              11.2
##   Adult_mortality Alcohol_consumption Hepatitis_B Measles  BMI Polio Diphtheria
## 1        105.8240                1.32          97      65 27.8    97         97
## 2         57.9025               10.35          97      94 26.0    97         97
## 3        201.0765                1.57          60      35 21.2    67         64
## 4        222.1965                5.68          93      74 25.3    92         93
## 5         57.9510                2.89          97      89 27.0    94         94
## 6         95.2200                4.19          88      86 26.4    89         89
##   Incidents_HIV GDP_per_capita Population_mln Thinness_ten_nineteen_years
## 1          0.08          11006          78.53                         4.9
## 2          0.09          25742          46.44                         0.6
## 3          0.13           1076        1183.21                        27.1
## 4          0.79           4146           0.75                         5.7
## 5          0.08          33995           7.91                         1.2
## 6          0.16           9110           4.35                         2.0
##   Thinness_five_nine_years Schooling Economy_status_Developed
## 1                      4.8       7.8                        0
## 2                      0.5       9.7                        1
## 3                     28.0       5.0                        0
## 4                      5.5       7.9                        0
## 5                      1.1      12.8                        1
## 6                      1.9       7.9                        0
##   Economy_status_Developing Life_expectancy
## 1                         1            76.5
## 2                         0            82.8
## 3                         1            65.4
## 4                         1            67.0
## 5                         0            81.7
## 6                         1            78.2
#

O čom je kvantilová regresia teória - príklad

Cieľom tohto Notebooku je predstaviť kvantilovú regresiu.

Táto metóda je užitočná, keď nechceme opísať iba priemerný vplyv vysvetľujúcej premennej na závislú premennú, ale aj to, ako sa tento vplyv líši pre pozorovania nachádzajúce sa v rôznych častiach podmieneného rozdelenia závislej premennej. Napríklad namiesto otázky:

Ako \(x\) ovplyvňuje priemer hodnoty \(y\)?

sa môžeme opýtať:

Ako \(x\) ovplyvňuje 10. percentil, medián alebo 90. percentil hodnoty \(y\)?

Toto je obzvlášť užitočné, keď:

1. Základná myšlienka kvantilovej regresie

V klasickej lineárnej regresii odhadovanej metódou OLS modelujeme podmienený priemer:

\[ E(y_i \mid x_i) = \beta_0 + \beta_1 x_i. \]

Naproti tomu v kvantilovej regresii modelujeme podmienený kvantil:

\[ Q_{\tau}(y_i \mid x_i) = \beta_0(\tau) + \beta_1(\tau) x_i, \]

kde:

Namiesto odhadu jednej regresnej priamky pre podmienený priemer môžeme odhadnúť niekoľko regresných priamok pre rôzne časti rozdelenia.

2. Prečo je to užitočné?

Predpokladajme, že príjem ovplyvňuje výdavky domácností. OLS nám hovorí, ako príjem mení priemerné výdavky. Kvantilová regresia však môže ukázať, či má príjem:

Týmto spôsobom nám kvantilová regresia umožňuje študovať heterogenitu efektov.

3. Princíp odhadu

Metóda najmenších štvorcov minimalizuje súčet druhých mocnín rezíduí:

\[ \min_{\beta} \sum_{i=1}^n (y_i - x_i'\beta)^2. \]

Kvantilová regresia používa inú funkciu strát, nazývanú kontrolná funkcia strát (alebo asymetrická absolútna strata):

\[ \min_{\beta} \sum_{i=1}^n \rho_{\tau}(y_i - x_i'\beta), \]

kde

\[ \rho_{\tau}(u) = u\big(\tau - I(u < 0)\big). \]

To znamená:

Dôležitým praktickým dôsledkom je, že kvantilová regresia je často odolnejšia voči odľahlým hodnotám v závislej premennej ako OLS.

4. OLS versus kvantilová regresia

OLS

OLS odpovedá na otázku:

O koľko sa zmení podmienený priemer premennej \(y\), keď sa premenná \(x\) zvýši o jednu jednotku?

Kvantilová regresia

Kvantilová regresia odpovedá na otázku:

O koľko sa zmení podmienený \(\tau\)-kvantil premennej \(y\), keď sa premenná \(x\) zvýši o jednu jednotku?

Koeficient sklonu v kvantilovej regresii sa preto interpretuje podobne ako metóda najmenších štvorcov (OLS), ale vzťahuje sa na konkrétny kvantil, a nie na priemer.

5. Príklad so simulovanými údajmi

V nasledujúcom príklade generujeme údaje, v ktorých sa rozptyl premennej \(y\) zvyšuje s premennou \(x\). To znamená, že podmienené rozdelenie sa rozširuje pre vyššie hodnoty premennej \(x\). Toto je situácia, v ktorej je kvantilová regresia veľmi informatívna.

# install.packages("quantreg")   # uncomment if needed
library(quantreg)
library(ggplot2)
library(dplyr)
set.seed(123)

n <- 300
x <- runif(n, 0, 10)

# Heteroskedastic data: variability grows with x
u <- rnorm(n, mean = 0, sd = 0.8 + 0.5 * x)
y <- 2 + 1.5 * x + u

data_qr <- data.frame(x = x, y = y)
head(data_qr)
##          x         y
## 1 2.875775  8.076534
## 2 7.883051 17.471011
## 3 4.089769  9.079732
## 4 8.830174  9.986489
## 5 9.404673 15.449741
## 6 0.455565  2.395162

Scatterplot údajov

ggplot(data_qr, aes(x = x, y = y)) +
  geom_point(alpha = 0.6) +
  labs(
    title = "Simulované údaje pre kvantilovú regresiu",
    x = "Vysvetľujúca premenná x",
    y = "Vysvetľovaná premenná y"
  ) +
  theme_minimal()

6. Odhad OLS a kvantilovej regresie

We first estimate the standard OLS model, and then quantile regressions for:

m_ols <- lm(y ~ x, data = data_qr)
m_q25 <- rq(y ~ x, tau = 0.25, data = data_qr)
m_q50 <- rq(y ~ x, tau = 0.50, data = data_qr)
m_q75 <- rq(y ~ x, tau = 0.75, data = data_qr)

výsledky OLS

summary(m_ols)
## 
## Call:
## lm(formula = y ~ x, data = data_qr)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -11.1828  -1.9437  -0.0017   1.7881  13.5687 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.05825    0.42229   4.874 1.78e-06 ***
## x            1.49394    0.07372  20.264  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.577 on 298 degrees of freedom
## Multiple R-squared:  0.5795, Adjusted R-squared:  0.5781 
## F-statistic: 410.6 on 1 and 298 DF,  p-value: < 2.2e-16

výsledky kvantilovej regresie

{r q25-summary2 summary(m_q25)

summary(m_q50)
## 
## Call: rq(formula = y ~ x, tau = 0.5, data = data_qr)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 2.32488      2.07794  2.59917 
## x           1.40400      1.35303  1.50171
summary(m_q75)
## 
## Call: rq(formula = y ~ x, tau = 0.75, data = data_qr)
## 
## tau: [1] 0.75
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 2.56628      2.31909  3.16484 
## x           1.82015      1.63000  1.93044

7. Porovnanie odhadnutých koeficientov

coef_table <- rbind(
  OLS = coef(m_ols),
  Q25 = coef(m_q25),
  Q50 = coef(m_q50),
  Q90 = coef(m_q75)
)

coef_table
##     (Intercept)        x
## OLS    2.058254 1.493939
## Q25    1.656366 1.154944
## Q50    2.324878 1.403998
## Q90    2.566281 1.820151

Interpretácia

Ak sa koeficienty sklonu líšia v rámci kvantilov, potom vplyv \(x\) nie je rovnomerný v rámci podmieneného rozdelenia \(y\).

Napríklad:

  • ak je sklon väčší v 75. percentile ako v 25. percentile,
  • potom vyššie hodnoty \(x\) sú spojené so silnejším nárastom v hornej časti rozdelenia ako v dolnej časti.

To naznačuje, že vplyv \(x\) je heterogénny.

8. Grafické porovnanie regresných priamok

ggplot(data_qr, aes(x = x, y = y)) +
  geom_point(alpha = 0.5) +
  geom_abline(
    intercept = coef(m_ols)[1],
    slope = coef(m_ols)[2],
    linewidth = 1
  ) +
  geom_abline(
    intercept = coef(m_q25)[1],
    slope = coef(m_q25)[2],
    linetype = "dashed",
    color = "blue",
    linewidth = 1
  ) +
  geom_abline(
    intercept = coef(m_q50)[1],
    slope = coef(m_q50)[2],
    linetype = "dotdash",
    color = "green",
    linewidth = 1
  ) +
  geom_abline(
    intercept = coef(m_q75)[1],
    slope = coef(m_q75)[2],
    linetype = "longdash",
    color = "red",
    linewidth = 1
  ) +
  labs(
    title = "Vyrovnané priamky OLS a kvanilových regresií",
    subtitle = "čierna = OLS, Modrá, zelená, červená - Q25,Q50, Q75",
    x = "x",
    y = "y"
  ) +
  theme_minimal()

9. Odhadovanie viacerých kvantilov naraz

Veľmi často chceme odhadnúť celú sadu kvantilov. Napríklad od 25. do 75. percentilu.

taus <- seq(0.1, 0.9, by = 0.1)
m_many <- rq(y ~ x, tau = taus, data = data_qr)
summary(m_many)
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.1
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 1.05638      0.83518  1.60419 
## x           0.86295      0.63656  0.97183 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.2
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 1.32109      0.64882  1.91774 
## x           1.14242      0.92742  1.22773 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.3
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 1.65101      1.11134  1.85084 
## x           1.24082      1.16062  1.38763 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.4
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 1.91165      1.40716  2.46532 
## x           1.38161      1.20784  1.44957 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 2.32488      2.07794  2.59917 
## x           1.40400      1.35303  1.50171 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.6
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 2.47811      2.02466  2.79854 
## x           1.51904      1.41448  1.64766 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.7
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 2.56724      2.33806  3.04354 
## x           1.66664      1.58339  1.83000 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.8
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 2.64593      2.31301  3.30625 
## x           1.92603      1.69915  2.07245 
## 
## Call: rq(formula = y ~ x, tau = taus, data = data_qr)
## 
## tau: [1] 0.9
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) 2.98539      2.54812  3.91652 
## x           2.15588      1.90586  2.34556

Graf koeficientov sklonu naprieč kvantilmi

coef_many <- data.frame(
  tau = taus,
  intercept = coef(m_many)[1, ],
  slope = coef(m_many)[2, ]
)

coef_many
##          tau intercept     slope
## tau= 0.1 0.1  1.056384 0.8629467
## tau= 0.2 0.2  1.321093 1.1424229
## tau= 0.3 0.3  1.651011 1.2408208
## tau= 0.4 0.4  1.911652 1.3816066
## tau= 0.5 0.5  2.324878 1.4039984
## tau= 0.6 0.6  2.478112 1.5190388
## tau= 0.7 0.7  2.567240 1.6666354
## tau= 0.8 0.8  2.645930 1.9260282
## tau= 0.9 0.9  2.985390 2.1558783
ggplot(coef_many, aes(x = tau, y = slope)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  geom_hline(yintercept = coef(m_ols)[2], linetype = "dashed") +
  labs(
    title = "Sklon koeficientu naprieč kvantilmi",
    subtitle = "čiarkovaná horizontálna čiara = OLS sklon",
    x = "Quantile (tau)",
    y = "Odhadnuté sklony regresných koeficientov"
  ) +
  theme_minimal()

Tento graf je veľmi užitočný pre výučbu, pretože ukazuje, či sa vplyv \(x\) mení v rámci rozdelenia \(y\).

10. Interpretácia koeficientov

Predpokladajme, že pre \(\tau = 0,90\) je odhadovaný sklon \(1,9\). Potom je interpretácia:

Zvýšenie \(x\) o jednu jednotku je spojené so zvýšením približne o 2,2 jednotky v podmienenom 75. percentile \(y\).

Toto nie je to isté ako tvrdenie, že 75% pozorovaní sa zvýši o 2,2 jednotky. Toto tvrdenie sa vzťahuje na podmienený kvantil závislej premennej.

11. Výhody kvantilovej regresie

Hlavné výhody:

  1. Poskytuje bohatší pohľad na vzťah medzi premennými.

  2. Je užitočná pri heteroskedasticite.

  3. Je informatívna, keď sa vplyv líši v rámci rozdelenia.

  4. Je menej citlivá na odľahlé hodnoty v závislej premennej ako OLS.

  5. Nevyžaduje normalitu rezíduí rovnakým spôsobom ako mnohé klasické postupy založené na OLS.

12. Obmedzenia

Treba spomenúť aj niektoré obmedzenia:

  1. Interpretácia je náročnejšia ako v OLS.

  2. Výsledky v extrémnych kvantiloch môžu byť v malých vzorkách menej stabilné.

13. Praktický záver

Kvantilná regresia je užitočným rozšírením klasickej regresnej analýzy. Zatiaľ čo OLS sa zameriava na podmienený priemer, kvantilová regresia nám umožňuje študovať celé podmienené rozdelenie závislej premennej.

Je obzvlášť cenná, keď:

14. Life-expectancy príklad

Príprava údajov - import z csv súboru

Súbor Life_Expectancy_Data obsahuje databázu determinantov očakávanej dĺžky života. Import údajov urobíme nasledovne

# import the dataset and create a data.frame udaje_svet
#
udaje_svet <- read.csv("udaje_svet/Life-Expectancy-Data-Updated.csv",header=TRUE,sep=",",dec=".",check.names = TRUE)
head(udaje_svet)
##      Country                        Region Year Infant_deaths Under_five_deaths
## 1    Turkiye                   Middle East 2015          11.1              13.0
## 2      Spain                European Union 2015           2.7               3.3
## 3      India                          Asia 2007          51.5              67.9
## 4     Guyana                 South America 2006          32.8              40.5
## 5     Israel                   Middle East 2012           3.4               4.3
## 6 Costa Rica Central America and Caribbean 2006           9.8              11.2
##   Adult_mortality Alcohol_consumption Hepatitis_B Measles  BMI Polio Diphtheria
## 1        105.8240                1.32          97      65 27.8    97         97
## 2         57.9025               10.35          97      94 26.0    97         97
## 3        201.0765                1.57          60      35 21.2    67         64
## 4        222.1965                5.68          93      74 25.3    92         93
## 5         57.9510                2.89          97      89 27.0    94         94
## 6         95.2200                4.19          88      86 26.4    89         89
##   Incidents_HIV GDP_per_capita Population_mln Thinness_ten_nineteen_years
## 1          0.08          11006          78.53                         4.9
## 2          0.09          25742          46.44                         0.6
## 3          0.13           1076        1183.21                        27.1
## 4          0.79           4146           0.75                         5.7
## 5          0.08          33995           7.91                         1.2
## 6          0.16           9110           4.35                         2.0
##   Thinness_five_nine_years Schooling Economy_status_Developed
## 1                      4.8       7.8                        0
## 2                      0.5       9.7                        1
## 3                     28.0       5.0                        0
## 4                      5.5       7.9                        0
## 5                      1.1      12.8                        1
## 6                      1.9       7.9                        0
##   Economy_status_Developing Life_expectancy
## 1                         1            76.5
## 2                         0            82.8
## 3                         1            65.4
## 4                         1            67.0
## 5                         0            81.7
## 6                         1            78.2
#
attach(udaje_svet)
udaje_2010 <- udaje_svet[udaje_svet$Year == 2010,]
m_ols <- lm(Life_expectancy ~ Alcohol_consumption+Schooling, data = udaje_svet)
m_q25 <- rq(Life_expectancy ~ Alcohol_consumption+Schooling, tau = 0.25, data = udaje_svet)
m_q50 <- rq(Life_expectancy ~ Alcohol_consumption+Schooling, tau = 0.50, data = udaje_svet)
m_q75 <- rq(Life_expectancy ~ Alcohol_consumption+Schooling, tau = 0.75, data = udaje_svet)
summary(m_ols)
## 
## Call:
## lm(formula = Life_expectancy ~ Alcohol_consumption + Schooling, 
##     data = udaje_svet)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.286  -3.658   1.221   4.735  13.945 
## 
## Coefficients:
##                     Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         52.06419    0.31319 166.241  < 2e-16 ***
## Alcohol_consumption -0.19727    0.03798  -5.195  2.2e-07 ***
## Schooling            2.32476    0.04768  48.759  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.375 on 2861 degrees of freedom
## Multiple R-squared:  0.5409, Adjusted R-squared:  0.5405 
## F-statistic:  1685 on 2 and 2861 DF,  p-value: < 2.2e-16

výsledky kvantilovej regresie

summary(m_q25)
## 
## Call: rq(formula = Life_expectancy ~ Alcohol_consumption + Schooling, 
##     tau = 0.25, data = udaje_svet)
## 
## tau: [1] 0.25
## 
## Coefficients:
##                     Value    Std. Error t value  Pr(>|t|)
## (Intercept)         46.93625  0.65772   71.36219  0.00000
## Alcohol_consumption -0.16586  0.06361   -2.60726  0.00917
## Schooling            2.47174  0.08935   27.66427  0.00000
summary(m_q50)
## 
## Call: rq(formula = Life_expectancy ~ Alcohol_consumption + Schooling, 
##     tau = 0.5, data = udaje_svet)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                     Value     Std. Error t value   Pr(>|t|) 
## (Intercept)          53.48395   0.40799  131.09256   0.00000
## Alcohol_consumption  -0.22454   0.04734   -4.74342   0.00000
## Schooling             2.31527   0.05642   41.03632   0.00000
summary(m_q75)
## 
## Call: rq(formula = Life_expectancy ~ Alcohol_consumption + Schooling, 
##     tau = 0.75, data = udaje_svet)
## 
## tau: [1] 0.75
## 
## Coefficients:
##                     Value     Std. Error t value   Pr(>|t|) 
## (Intercept)          60.49015   0.30462  198.57277   0.00000
## Alcohol_consumption  -0.05500   0.02376   -2.31506   0.02068
## Schooling             1.76661   0.03696   47.79383   0.00000
m_ols <- lm(Life_expectancy ~ Alcohol_consumption+Schooling, data = udaje_2010)
m_q25 <- rq(Life_expectancy ~ Alcohol_consumption+Schooling, tau = 0.25, data = udaje_2010)
m_q50 <- rq(Life_expectancy ~ Alcohol_consumption+Schooling, tau = 0.50, data = udaje_2010)
m_q75 <- rq(Life_expectancy ~ Alcohol_consumption+Schooling, tau = 0.75, data = udaje_2010)
coef_table <- rbind(
  OLS = coef(m_ols),
  Q25 = coef(m_q25),
  Q50 = coef(m_q50),
  Q75 = coef(m_q75)
)

coef_table
##     (Intercept) Alcohol_consumption Schooling
## OLS    53.35024          -0.2323919  2.219635
## Q25    48.04923          -0.2777595  2.447294
## Q50    54.40012          -0.1462916  2.190997
## Q75    61.42832          -0.0813932  1.703866
ggplot(udaje_svet, aes(x = Alcohol_consumption, y = Life_expectancy)) +
  geom_point(alpha = 0.5) +
  geom_abline(
    intercept = coef(m_ols)[1],
    slope = coef(m_ols)[2],
    linewidth = 1
  ) +
  geom_abline(
    intercept = coef(m_q25)[1],
    slope = coef(m_q25)[2],
    linetype = "dashed",
    color = "blue",
    linewidth = 1
  ) +
  geom_abline(
    intercept = coef(m_q50)[1],
    slope = coef(m_q50)[2],
    linetype = "dotdash",
    color = "green",
    linewidth = 1
  ) +
  geom_abline(
    intercept = coef(m_q75)[1],
    slope = coef(m_q75)[2],
    linetype = "longdash",
    color = "red",
    linewidth = 1
  ) +
  labs(
    title = "Vyrovnané priamky OLS a kvanilových regresií",
    subtitle = "čierna = OLS, Modrá, zelená, červená - Q25,Q50, Q75",
  x = "x",
  y ="y"
  ) +
  theme_minimal()

15. Záverečná poznámka

V aplikovanej ekonometrii sa kvantilová regresia často používa pri analýze:

Jej hlavnou silnou stránkou je, že umožňuje ekonometrovi povedať nielen čo sa deje v priemere, ale aj čo sa deje v rôznych častiach rozdelenia.