Raport4
Maja i Zuzia
2024-11-19
Regresja kwantylowa
Waszym zadaniem dzisiaj jest zamodelowanie - porównanie KMNK oraz regresji kwantylowej (różno-poziomowej) dla zmiennej “earnings” - wynagrodzenia. Dobierz i przetestuj predyktory, kwantyle dla modeli. Wykonaj testy różnic współczynnikow dla finalnych modeli.
library(CVXR)
library(AER)
library(stargazer)
library(WRTDStidal)
library(tidyverse)
library(kableExtra)
library(quantreg)
library(PogromcyDanych)
library(dplyr)
library(ggplot2)
library(sjPlot)
library(lmtest)
library(ggridges)
library(tseries)data("CPSSW9298")
dane92 <- CPSSW9298 %>%
filter(year=='1992')
dane98 <- CPSSW9298 %>%
filter(year=='1998')Dane Stock i Watson (2007) przedstawiają kilka podzbiorów utworzonych z marcowych badań Current Population Surveys (CPS) z danymi na temat związku zarobków i wykształcenia dla dwóch okresów: 1992 i 1998.
1992
Zależność wielkości wynagrodzenia od wieku ankietera
ggplot(dane92, aes(age, earnings))+
geom_point()+
geom_smooth(method=lm, col='coral1')+
geom_quantile(quantiles=0.5, col='steelblue')+
scale_x_continuous(breaks = seq(min(dane92$age), max(dane92$age), by = 1))+
labs(title='Wykres rozrzutu: zarobki ~ wiek')+
theme_light()
Czerwonym kolorem zaznaczona została regresja OLS a niebieskim regresja medianowa. Dają one zbliżone wyniki, natomiast nie wykazują prawidłowo zależności między wynagrodzeniem a wiekiem ankietera.
Dla szacowań OLS częsciowo wynika to z braku stałej wariancji dla reszt. Dla testu Breuscha-Pagana badającego homoskedastyczność reszt p-value jest poniżej poziomu istotności (0.05), zatem odrzucona została hipoteza zerowa na rzecz alternatywnej, mówiącej o wystąpieniu heteroskedastyczności.
W przypadku regresji medianowej, jest ona bardziej odporna na wartości odstające, natomiast sama, nie wyjaśnia dobrze zależności w rozkładzie.
lm_model <- lm(earnings~age, data=dane92)
mr_model <- rq(earnings~age, data=dane92, tau=.5)
bptest(lm_model)##
## studentized Breusch-Pagan test
##
## data: lm_model
## BP = 89.516, df = 1, p-value < 2.2e-16AIC <- AIC(lm_model,mr_model)| df | AIC | |
|---|---|---|
| lm_model | 3 | 47384.93 |
| mr_model | 2 | 47324.48 |
Regresja kwantylowa nie wymaga założenia o homoskedastyczności. Dla danych o zarobkach i wieku, zastosowano kwantyle 0.2, 0.4, 0.6 oraz 0.8. Dzięki temu, możemy zbadać np. jaki kształtują się zarobki w dolnym ogonie rozkładu wieku, czyli dla osób w wieku 25 i 26 lat.
ggplot(dane92, aes(age, earnings))+
geom_point()+
geom_smooth(method=lm, col='coral1')+
geom_quantile(quantiles=seq(.2,0.8, by=.2), col='grey',size=2,alpha=.3)+
geom_quantile(quantiles=0.5, col='steelblue')+
scale_x_continuous(breaks = seq(min(dane92$age), max(dane92$age), by = 1))+
labs(title='Model kwantylowy dla wieku')+
theme_light()
Rozkład zarobków według płci
ggplot(dane92, aes(x = earnings, fill = gender)) +
geom_density(alpha = 0.6, position = "identity", bins = 30) +
labs(title = "Rozkład zarobków według płci")+
theme_light() 
Zarówno dla mężczyzn jak i dla kobiet, rozkład zarobków jest odchylony w lewą stronę, czyli jest prawoskośny (brak normalności rozkładu).
lm <- lm(earnings~gender, data=dane92)
qm25 <- rq(earnings~gender, data=dane92, tau=.25)
qm50 <- rq(earnings~gender, data=dane92, tau=.5)
qm75 <- rq(earnings~gender, data=dane92, tau=.75)theme_set(theme_bw())
plot_models(lm, qm25,qm50,qm75,
show.values=TRUE, m.labels=c('OLS','QR 25%','QR 50%','QR 75%'),
legend.title='Model',title='Różnica w zarobkach kobiet i mężczyzn' )
Analizując powyższy wykres, dzięki regresji kwantylowej możemy zobaczyć o ile mniej zarabiają średnio kobiety w różnych przedziałach dochodów. Wraz ze wzrostem dochodów, różnica pomiędzy zarobkami kobiet, a mężczyzn rośnie, zatem luka płacowa powiększa się.
Regresja kwantylowa dla kilku zmiennych objaśniających
lm <- lm(log(earnings) ~ degree+gender+age, data = dane92)
q20 <- rq(log(earnings) ~ degree+gender+age, data = dane92, tau = 0.20)
q40 <- rq(log(earnings) ~ degree+gender+age, data = dane92, tau = 0.40)
q60 <- rq(log(earnings) ~ degree+gender+age, data = dane92, tau = 0.60)
q80 <- rq(log(earnings) ~ degree+gender+age, data = dane92, tau = 0.80)
# Tabela z porównaniem wyników trzech modeli:
stargazer(q20, q40, q60,q80, lm, title = "Wyniki regresji kwantylowych", type = "text", column.labels=c('QR 20%','QR 20%','QR 60%','QR 80%','OLS'))##
## Wyniki regresji kwantylowych
## =====================================================================================
## Dependent variable:
## -----------------------------------------------------------------
## log(earnings)
## quantile OLS
## regression
## QR 20% QR 20% QR 60% QR 80% OLS
## (1) (2) (3) (4) (5)
## -------------------------------------------------------------------------------------
## degreebachelor 0.402*** 0.403*** 0.376*** 0.364*** 0.375***
## (0.017) (0.013) (0.012) (0.012) (0.010)
##
## genderfemale -0.141*** -0.174*** -0.187*** -0.201*** -0.167***
## (0.016) (0.013) (0.012) (0.012) (0.010)
##
## age 0.020*** 0.027*** 0.031*** 0.032*** 0.026***
## (0.003) (0.002) (0.002) (0.002) (0.002)
##
## Constant 1.300*** 1.392*** 1.484*** 1.696*** 1.486***
## (0.086) (0.070) (0.062) (0.065) (0.053)
##
## -------------------------------------------------------------------------------------
## Observations 7,590 7,590 7,590 7,590 7,590
## R2 0.187
## Adjusted R2 0.187
## Residual Std. Error 0.430 (df = 7586)
## F Statistic 582.296*** (df = 3; 7586)
## =====================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01*Logarym przy zmiennej objaśnane wynika ze skośności danych.
Wszystkie zmienne objaśniające w modelu kwantylowym wykazują wysoką istotność, zatem wnoszą do modelu istotne informacje. Porównajmy zatem modele kwantylowe z modelem oszacowanym metodą najmniejszych kwadratów.
Porównaie reszt z regresji KMNK i kwantylowej:
## model kwantylowy
model1 <- rq(log(earnings) ~ degree+gender+age,tau = 0.20, data = dane92)
reszty1 <- resid(model1)
## bezwarunkowy (pusty) model kwantylowy
model2 <- rq(log(earnings) ~ 1, tau = 0.20,data=dane92)
reszty2 <- resid(model2)
r_rq<-round(goodfit(reszty1, reszty2, 0.20),4)
## r2 modelu KMNK dla porównania
model_lm <- lm(log(earnings) ~ degree+gender+age,data = dane92)
r_lm<-round(summary(model_lm)$r.squared,4)| Metric | Value |
|---|---|
| Pseudo R dla regresji kwantylowej | 0.0875 |
| R kwadrat dla KMNK | 0.1872 |
R kwardat dla regresji kwantylowej jest niski, ponieważ dane są skokowe (wiek).
Testy statystyczne dla regresji liniowej (OLS):
Liniowość - Test RESET Ramseya
Nie ma podstaw do odrzucenia H0: model liniowy, zgodny ze specyfikacjąresettest(lm)## ## RESET test ## ## data: lm ## RESET = 1.0356, df1 = 2, df2 = 7584, p-value = 0.3551Homoskedastyczność - Test Breuscha Pagana
Odrzucamy H0 na rzecz H1: Zmienność składnika losowego nie jest zależny od x-ów (heteroskedastyczność)bptest(lm)## ## studentized Breusch-Pagan test ## ## data: lm ## BP = 28.897, df = 3, p-value = 2.354e-06Normalność składnika losowego - Test Jarque’a-Berego
Odrzucamy H0 na rzecz H1: składnik losowy nie ma rozkładu normalnyjarque.bera.test(lm$residuals)## ## Jarque Bera Test ## ## data: lm$residuals ## X-squared = 228, df = 2, p-value < 2.2e-16Autokorelacja składnika losowego - Test Ljung-Boxa
Odrzucamy H0 na rzecz H1: występuje autokorelacja składnika losowegoBox.test(model1$residuals, type="Ljung-Box")## ## Box-Ljung test ## ## data: model1$residuals ## X-squared = 109.68, df = 1, p-value < 2.2e-16
Kryterium informacyjne Akaike:
modele <- list(lm,q20,q40,q60,q80)
modele_n <- c('OLS','QR 20%','QR 20%','QR 60%','QR 80%')
aic <- sapply(modele, AIC)
AIC <- data.frame(Model = modele_n,AIC = aic)| Model | AIC |
|---|---|
| OLS | 8750.943 |
| QR 20% | 11540.259 |
| QR 20% | 9645.561 |
| QR 60% | 9135.309 |
| QR 80% | 9891.324 |
Regresja kwantylowa dla 0.60 kwantyla okazuje się najniższa, co wskazuje na najlepze dopasowanie do danych.
1998
Rozkład zarobków
hist(dane98$earnings, main = "Rozkład zarobków w 1998 roku", xlab = "Zarobki", ylab = "Częstotliwość", col = "lightblue")
Na podstawie histogramu można zauważyć, że rozkład zarobków w roku 1998 był asymetryczny prawostronnie i większość osób zarabiała stosunkowo niskie kwoty, a tylko nieliczne osiągały znacznie wyższe dochody. Dominują niskie zarobki, przy czym największa częstotliwość przypada na przedziały w okolicach 5–15. Jednocześnie obecność długiego ogona po prawej stronie sugeruje istnienie niewielkiej liczby osób z bardzo wysokimi zarobkami, znacznie przekraczającymi średnią, co może być oznaką nierówności dochodowych.
Zależność wielkości wynagrodzenia od wieku ankietera
ggplot(dane98, aes(age, earnings))+
geom_point()+
geom_smooth(method=lm, col='violetred')+
geom_quantile(quantiles=c( 0.25, 0.5, 0.75, 0.90, 0.95, 0.99), col='cornflowerblue')+
scale_x_continuous(breaks = seq(min(dane92$age), max(dane98$age), by = 1))+
labs(title='Wykres rozrzutu: zarobki ~ wiek')+
theme_light()
Wykres rozrzutu pokazuje łagodny wzrost zarobków wraz z wiekiem, jednak związek ten nie jest silny, co potwierdza duży rozrzut danych w każdej grupie wiekowej. Linie kwantyli wskazują na znaczną rozpiętość zarobków – zarówno osoby młodsze (25–28 lat), jak i starsze (33–34 lata) mogą osiągać zarówno niskie, jak i bardzo wysokie dochody. Wyższe kwantyle (np. 0.95) pokazują, że wiek nie jest kluczowym czynnikiem dla najwyższych zarobków. Wiek ma ograniczony wpływ na zarobki, a ich zróżnicowanie może być bardziej zależne od innych czynników, takich jak wykształcenie, doświadczenie czy branża.
lm_model2 <- lm(earnings~age, data=dane98)
mr_model2 <- rq(earnings~age, data=dane98, tau=0.5)
bptest(lm_model)##
## studentized Breusch-Pagan test
##
## data: lm_model
## BP = 89.516, df = 1, p-value < 2.2e-16AIC <- AIC(lm_model,mr_model)| df | AIC | |
|---|---|---|
| lm_model | 3 | 47384.93 |
| mr_model | 2 | 47324.48 |
Sytacja przedstawionych modeli wygląda podobnie jak dla roku 1992, w obu przypadkach AIC jest bardzo zbliżone, z nieznacznie wyższą wartością dla modelu liniowego.
Rozkład zarobków według płci
ggplot(dane98, aes(x = earnings, fill = gender)) +
geom_density(alpha = 0.6, position = "identity", bins = 30) +
labs(title = "Rozkład zarobków według płci")+
theme_light() 
Rozkład zarobków według płci nie zmienił się znacznie od tego, co można było zauważyć w roku 1992. Zarówno w przypadku kobiet jak i mężczyn występuje brak normalności rozkładu i widoczna jest prawoskośność. W 1998 roku nadal w oczy rzucaja się nierówności płacowe kobiet i mężczyzn. Kobiety częściej zarabiają niskie kwoty i niewiele z nich osiąga bardzo wysokie zarobki - na płace w przedziale 40-50 mogą liczyć głównie mężczyźni.
lm2 <- lm(earnings~gender, data=dane98)
qm25 <- rq(earnings~gender, data=dane98, tau=.25)
qm50 <- rq(earnings~gender, data=dane98, tau=.5)
qm75 <- rq(earnings~gender, data=dane98, tau=.75)theme_set(theme_bw())
plot_models(lm2, qm25,qm50,qm75,
show.values=TRUE, m.labels=c('OLS','QR 25%','QR 50%','QR 75%'),
legend.title='Model',title='Różnica w zarobkach kobiet i mężczyzn' )
Na wykresie przedstawiono różnice w zarobkach kobiet i mężczyzn w różnych modelach regresji (OLS oraz kwantylowych dla kwantyli 25%, 50%, 75%). Wszystkie modele wskazują, że kobiety zarabiają mniej od mężczyzn, ale skala tej różnicy maleje w wyższych kwantylach zarobków. W najniższych zarobkach (QR 25%) różnica wynosi -2,99 jednostki, a w najwyższych (QR 75%) -1,44 jednostki. Oznacza to, że różnice w zarobkach są większe w dolnych partiach rozkładu zarobków, a mniejsze w górnych. Istotność statystyczna wyników sugeruje, że różnice są silnie znaczące we wszystkich przypadkach.
Regresja kwantylowa dla kilku zmiennych objaśniających
q25 <- rq(earnings~ age+gender+degree, data = dane98, tau = 0.25)
q50 <- rq(earnings ~ age+gender+degree, data = dane98, tau = 0.50)
q75 <- rq(earnings~ age+gender+degree, data = dane98, tau = 0.75)
# Tabela z porównaniem wyników trzech modeli:
stargazer(q25, q50, q75, title = "Wyniki regresji kwantylowych", type = "text",column.labels=c('QR 25%','QR 50%','QR 75%'))##
## Wyniki regresji kwantylowych
## ============================================
## Dependent variable:
## -----------------------------
## earnings
## QR 25% QR 50% QR 75%
## (1) (2) (3)
## --------------------------------------------
## age 0.184*** 0.302*** 0.481***
## (0.027) (0.030) (0.043)
##
## genderfemale -1.643*** -2.404*** -3.365***
## (0.150) (0.170) (0.242)
##
## degreebachelor 3.803*** 4.986*** 6.250***
## (0.166) (0.182) (0.280)
##
## Constant 3.124*** 2.830*** 1.442
## (0.792) (0.894) (1.293)
##
## --------------------------------------------
## Observations 5,911 5,911 5,911
## ============================================
## Note: *p<0.1; **p<0.05; ***p<0.01Po zastosowaniu najczęściej używanego podziału kwantylowego okazuje się, że stała dla ostatniego kwantyla jest nieistotna statystycznie. Sugerując się analizą danych za rok 1992, może okazać się, że lepsze modele dałoby się uzyskać przy zastosowaniu innego podziału kwantylowego oraz logarytmowanie zmiennej objaśnianej.
lm <- lm(log(earnings) ~ degree+gender+age, data = dane98)
q20 <- rq(log(earnings) ~ degree+gender+age, data = dane98, tau = 0.20)
q40 <- rq(log(earnings) ~ degree+gender+age, data = dane98, tau = 0.40)
q60 <- rq(log(earnings) ~ degree+gender+age, data = dane98, tau = 0.60)
q80 <- rq(log(earnings) ~ degree+gender+age, data = dane98, tau = 0.80)
# Tabela z porównaniem wyników modeli:
stargazer(q20, q40, q60,q80, lm, title = "Wyniki regresji kwantylowych", type = "text", column.labels=c('QR 20%','QR 20%','QR 60%','QR 80%','OLS'))##
## Wyniki regresji kwantylowych
## =====================================================================================
## Dependent variable:
## -----------------------------------------------------------------
## log(earnings)
## quantile OLS
## regression
## QR 20% QR 20% QR 60% QR 80% OLS
## (1) (2) (3) (4) (5)
## -------------------------------------------------------------------------------------
## degreebachelor 0.413*** 0.394*** 0.379*** 0.377*** 0.383***
## (0.016) (0.015) (0.013) (0.016) (0.012)
##
## genderfemale -0.169*** -0.179*** -0.197*** -0.211*** -0.180***
## (0.016) (0.015) (0.013) (0.015) (0.012)
##
## age 0.018*** 0.023*** 0.026*** 0.026*** 0.021***
## (0.003) (0.003) (0.002) (0.003) (0.002)
##
## Constant 1.504*** 1.663*** 1.783*** 2.043*** 1.788***
## (0.083) (0.079) (0.070) (0.083) (0.062)
##
## -------------------------------------------------------------------------------------
## Observations 5,911 5,911 5,911 5,911 5,911
## R2 0.183
## Adjusted R2 0.182
## Residual Std. Error 0.446 (df = 5907)
## F Statistic 440.464*** (df = 3; 5907)
## =====================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Po zastosowaniu wspomnianych zmian, wszystkie zmienne w powstałych modelach okazały się być istotne statystycznie.
Porównanie modeli kwantylowych oraz KMNK
## model kwantylowy
model3 <- rq(log(earnings) ~ degree+gender+age,tau = 0.20, data = dane98)
reszty3 <- resid(model3)
## bezwarunkowy (pusty) model kwantylowy
model4 <- rq(log(earnings) ~ 1, tau = 0.20,data=dane98)
reszty4 <- resid(model4)
r_rq2<-round(goodfit(reszty3, reszty4, 0.20),4)
## r2 modelu KMNK dla porównania
model_lm2 <- lm(log(earnings) ~ degree+gender+age,data = dane98)
r_lm2<-round(summary(model_lm2)$r.squared,4)| Metric | Value |
|---|---|
| Pseudo R dla regresji kwantylowej | 0.1011 |
| R kwadrat KMNK | 0.1828 |
Wyliczone R kwardat ponownie okazało się niskie, ze względu na charakter danych. W przypadku roku 1998,Psedu R kwadrat jest wyższe niż w roku 1992, jednak nadal niższe niż dla metody KMNK.
Kryterium informacyjne Akaike:
modele <- list(lm,q20,q40,q60,q80)
modele_n <- c('OLS','QR 20%','QR 20%','QR 60%','QR 80%')
aic <- sapply(modele, AIC)
AIC <- data.frame(Model = modele_n,AIC = aic)AIC %>%
kable("html", caption = "Porównanie poziomów AIC") %>%
kable_styling(full_width = FALSE)| Model | AIC |
|---|---|
| OLS | 7228.160 |
| QR 20% | 9221.416 |
| QR 20% | 7832.731 |
| QR 60% | 7517.501 |
| QR 80% | 8301.913 |
Model OLS ma najmniejszy AIC, co sugeruje, że pasuje do estymacji naszych danych.
Podsumowanie
Podsumowując, w wielu przypadkach regresja kwantylowa może dostarczyć bardziej wiarygodnych i szczegółowych wniosków niż regresja liniowa, szczególnie gdy dane są heteroskedastyczne, nie mają rozkładu normalnego lub zawierają wartości odstające. W analizowanym przypadku regresji wielorakiej, modele dają zbliżone wyniki, natomiast OLS nie spełnia większości założeń.
Dodatkowo, regresja kwantylowa dostarcza bardziej kompletnego obrazu zależności między zmiennymi i pozwala uzyskać głębsze zrozumienie badanego zjawiska, dokładnie jak w przypadku analizy zależności zarobków do płci.