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)
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
#
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ď:
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.
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.
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.
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 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.
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
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()
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)
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
{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
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
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:
To naznačuje, že vplyv \(x\) je heterogénny.
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()
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
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\).
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.
Hlavné výhody:
Poskytuje bohatší pohľad na vzťah medzi premennými.
Je užitočná pri heteroskedasticite.
Je informatívna, keď sa vplyv líši v rámci rozdelenia.
Je menej citlivá na odľahlé hodnoty v závislej premennej ako OLS.
Nevyžaduje normalitu rezíduí rovnakým spôsobom ako mnohé klasické postupy založené na OLS.
Treba spomenúť aj niektoré obmedzenia:
Interpretácia je náročnejšia ako v OLS.
Výsledky v extrémnych kvantiloch môžu byť v malých vzorkách menej stabilné.
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ď:
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
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()
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.