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/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_q10 <- rq(y ~ x, tau = 0.10, data = data_qr)
m_q50 <- rq(y ~ x, tau = 0.50, data = data_qr)
m_q90 <- rq(y ~ x, tau = 0.90, 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
summary(m_q10)
##
## Call: rq(formula = y ~ x, tau = 0.1, 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
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_q90)
##
## Call: rq(formula = y ~ x, tau = 0.9, 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_table <- rbind(
OLS = coef(m_ols),
Q10 = coef(m_q10),
Q50 = coef(m_q50),
Q90 = coef(m_q90)
)
coef_table
## (Intercept) x
## OLS 2.058254 1.4939386
## Q10 1.056384 0.8629467
## Q50 2.324878 1.4039984
## Q90 2.985390 2.1558783
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_q10)[1],
slope = coef(m_q10)[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_q90)[1],
slope = coef(m_q90)[2],
linetype = "longdash",
color = "red",
linewidth = 1
) +
labs(
title = "Vyrovnané priamky OLS a kvanilových regresií",
subtitle = "čierna = OLS, Modrá, zelená, červená - Q10,Q50, Q90",
x = "x",
y = "y"
) +
theme_minimal()
Veľmi často chceme odhadnúť celú sadu kvantilov. Napríklad od 10. do 90. 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 1,9 jednotky v podmienenom 90. percentile \(y\).
Toto nie je to isté ako tvrdenie, že 90 % pozorovaní sa zvýši o 1,9 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
#
udaje_svet <- read.csv("udaje/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_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)
summary(m_ols)
##
## Call:
## lm(formula = Life_expectancy ~ Alcohol_consumption + Schooling,
## data = udaje_2010)
##
## Residuals:
## Min 1Q Median 3Q Max
## -20.041 -3.505 1.081 4.166 12.041
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 53.3502 1.2244 43.572 <2e-16 ***
## Alcohol_consumption -0.2324 0.1542 -1.507 0.134
## Schooling 2.2196 0.1866 11.892 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.999 on 176 degrees of freedom
## Multiple R-squared: 0.5437, Adjusted R-squared: 0.5385
## F-statistic: 104.9 on 2 and 176 DF, p-value: < 2.2e-16
summary(m_q25)
##
## Call: rq(formula = Life_expectancy ~ Alcohol_consumption + Schooling,
## tau = 0.25, data = udaje_2010)
##
## tau: [1] 0.25
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 48.04923 43.92813 51.49500
## Alcohol_consumption -0.27776 -0.67920 0.23554
## Schooling 2.44729 1.79355 2.83861
summary(m_q50)
##
## Call: rq(formula = Life_expectancy ~ Alcohol_consumption + Schooling,
## tau = 0.5, data = udaje_2010)
##
## tau: [1] 0.5
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 54.40012 53.19757 55.68274
## Alcohol_consumption -0.14629 -0.47345 0.05897
## Schooling 2.19100 2.04482 2.46750
summary(m_q75)
##
## Call: rq(formula = Life_expectancy ~ Alcohol_consumption + Schooling,
## tau = 0.75, data = udaje_2010)
##
## tau: [1] 0.75
##
## Coefficients:
## coefficients lower bd upper bd
## (Intercept) 61.42832 57.15053 63.96004
## Alcohol_consumption -0.08139 -0.52826 0.19289
## Schooling 1.70387 1.30309 2.33771
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
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.