Nieklasyczne metody statystyki

Regresja kwantylowa

Dlaczego kwantylowa?

Dlaczego potrzebujemy regresji kwantylowej (QR)?

W szczególności, QR:

• jest odporna na punkty odstające i wpływowe

• nie zakłada stałej wariancji (znanej jako homoskedastyczność) dla zmiennej odpowiedzi lub reszt

• nie zakłada normalności ale główną zaletą QR w porównaniu z regresją liniową (LR) jest to, że QR bada różne wartości zmiennej odpowiedzi, a nie tylko średnią, i dostarcza w związku z tym pełniejszego obrazu związków między zmiennymi!

Wprowadzenie

Regresja kwantylowa (ang. quantile regression) została zaproponowana przez Koenkera i Bassetta (1978). Szczególny przypadek regresji kwantylowej dla kwantyla rzędu 0,5 (czyli mediany) jest równoważny estymatorowi LAD (ang. Least Absolute Deviation) – minimalizuje sumę bezwzględnych błędów.
Wprowadzenie różnych kwantyli regresji daje pełniejszy opis rozkładów warunkowych zwłaszcza w przypadku rozkładów asymetrycznych lub uciętych.

Regresja kwantylowa jest kolejną wariacją na temat najmniejszych kwadratów . Stratą jest współczynnik \(l_1\) funkcji:

\[ \phi(u) = \tau\max(u,0) - (1-\tau)\max(-u,0) = \frac{1}{2}|u| + \left(\tau - \frac{1}{2}\right)u, \]

gdzie \(\tau \in (0,1)\) oznacza konkretny kwantyl. Problemem jak poprzednio jest minimalizacja całkowitej straty resztowej. Model ten jest powszechnie stosowany w ekologii, ochronie zdrowia i innych dziedzinach, gdzie sama średnia nie wystarcza do uchwycenia złożonych zależności między zmiennymi.

Wymagania

Wymagana jest jedna liczbowa zmienna zależna. Zmienna przewidywana musi być zmienną ilościową. Predyktory mogą być zmiennymi ilościowymi lub sztucznymi zmiennymi w przypadku predyktorów jakościowych. Aby można było uruchomić analizę, wymagany jest wyraz wolny lub co najmniej jeden predyktor.

Regresja kwantylowa nie czyni założeń dotyczących rozkładu zmiennej przewidywanej i jest odporna na wpływ obserwacji odstających.

Analiza kwantylowa jest pokrewna regresji metodą najmniejszych kwadratów.

Przykład 1.

Wykorzystamy przykład z pakietu quantreg.

Jaki jest związek między całkowitym dochodem gospodarstwa domowego a odsetkiem dochodów wydatkowanych na żywność? Prawo Engela w ekonomii głosi, że w miarę wzrostu dochodów, część dochodów wydatkowanych na żywność spada, nawet jeśli wydatki na żywność bezwzględnie rosną. Stosując regresję kwantylową do tych danych, można określić, jakie wydatki na żywność ponosi 90% rodzin (dla 100 rodzin z danym dochodem), gdy nie interesują nas średnie wydatki na żywność.

Dane, które wykorzystamy - to zbiór “engel” - dane dotyczące wydatków na żywność. Jest to zbiór danych regresyjnych składający się z 235 obserwacji dotyczących dochodów i wydatków na żywność dla belgijskich gospodarstw domowych klasy robotniczej.

Powyższy wykres przedstawia dopasowanie regresji kwantylowej dla \(\tau = (0.1, 0.25, 0.5, 0.75, 0.90, 0.95)\). Dopasowanie KMNK to gruba czarna linia.

Poniżej znajduje się tabela z oszacowanymi współczynnikami.

knitr::kable(fits, format = "html", caption = "Oszacowania z KMNK oraz `quantreg`") %>%
    kable_styling("striped") %>%
    column_spec(1:8, background = "#ececec")

Oszacowania z KMNK oraz quantreg

	OLS	\(\tau_{0.10}\)	\(\tau_{0.25}\)	\(\tau_{0.50}\)	\(\tau_{0.75}\)	\(\tau_{0.90}\)	\(\tau_{0.95}\)
(Intercept)	147.4753885	110.1415742	95.4835396	81.4822474	62.3965855	67.3508721	64.1039632
income	0.4851784	0.4017658	0.4741032	0.5601806	0.6440141	0.6862995	0.7090685

Ok, możemy to zrobić bardziej przejrzyście i sformatować w ładnej tabeli wyników:

## 
## Wyniki regresji kwantylowych
## ==========================================
##                   Dependent variable:     
##              -----------------------------
##                         foodexp           
##                 (1)       (2)       (3)   
## ------------------------------------------
## income       0.474***  0.560***  0.644*** 
##               (0.029)   (0.028)   (0.023) 
##                                           
## Constant     95.484*** 81.482*** 62.397***
##              (21.392)  (19.251)  (16.305) 
##                                           
## ------------------------------------------
## Observations    235       235       235   
## ==========================================
## Note:          *p<0.1; **p<0.05; ***p<0.01

Finalnie, zaprezentujmy wyłącznie te 3 modele na wykresie:

Przykład 2.

Tutaj przeprowadzimy testy użycia pakietu quantreg, wykorzystując wbudowany zbiór danych “mtcars”. Zmienna “mpg” oznacza spalanie samochodów (mile/galon).

Zamodulejmy zależność regresyjną dla tej zmiennej od kilku predyktorów.

Najpierw oszacujmy regresję KMNK:

kmnk <- lm(mpg ~ disp + hp + factor(am) + factor(vs), data = mtcars)
summary(kmnk)

## 
## Call:
## lm(formula = mpg ~ disp + hp + factor(am) + factor(vs), data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.7981 -1.9532  0.0111  1.5665  5.6321 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 24.832119   2.890418   8.591 3.32e-09 ***
## disp        -0.008304   0.010087  -0.823  0.41757    
## hp          -0.037623   0.013846  -2.717  0.01135 *  
## factor(am)1  4.419257   1.493243   2.960  0.00634 ** 
## factor(vs)1  2.052472   1.627096   1.261  0.21794    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.812 on 27 degrees of freedom
## Multiple R-squared:  0.8104, Adjusted R-squared:  0.7823 
## F-statistic: 28.85 on 4 and 27 DF,  p-value: 2.13e-09

Teraz oszacujmy warunkowe regresje kwantylowe na różnych kwantylach, błąd standardowy uzyskany przez bootstrap.

Zauważ, że istnieje gradient we współczynnikach kwantylowych hp, jak również disp. Znak disp odwraca się, również współczynnik na czynniku am jest różny w zależności od kwantyli:

kwantyle <- c(0.25, 0.50, 0.75)
reg_kwantylowa <- rq(mpg ~ disp + hp + factor(am),tau = kwantyle,data = mtcars)
summary(reg_kwantylowa, se = "boot")

## 
## Call: rq(formula = mpg ~ disp + hp + factor(am), tau = kwantyle, data = mtcars)
## 
## tau: [1] 0.25
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 25.34665  1.59892   15.85231  0.00000
## disp        -0.02441  0.00766   -3.18815  0.00351
## hp          -0.01672  0.01596   -1.04767  0.30375
## factor(am)1  1.39719  1.32335    1.05580  0.30008
## 
## Call: rq(formula = mpg ~ disp + hp + factor(am), tau = kwantyle, data = mtcars)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 27.49722  1.83557   14.98024  0.00000
## disp        -0.02253  0.01598   -1.40963  0.16966
## hp          -0.02713  0.02389   -1.13586  0.26564
## factor(am)1  3.37328  2.07836    1.62305  0.11578
## 
## Call: rq(formula = mpg ~ disp + hp + factor(am), tau = kwantyle, data = mtcars)
## 
## tau: [1] 0.75
## 
## Coefficients:
##             Value    Std. Error t value  Pr(>|t|)
## (Intercept) 28.06384  1.80366   15.55938  0.00000
## disp         0.00445  0.01399    0.31813  0.75275
## hp          -0.06662  0.01762   -3.78175  0.00075
## factor(am)1  7.91402  2.30915    3.42724  0.00190

Testy współczynników

Użyjemy funkcji rq.anova z pakietu regresji kwantylowej, aby przeprowadzić test WALDA. Pamiętaj, że test WALDA mówi, że biorąc pod uwagę nieograniczone oszacowania modelu, przetestujemy hipotezę zerową mówiącą, że współczynniki spełniają pewne liniowe ograniczenia.

Aby ją przetestować, użyjemy obiektu zwróconego z uruchomienia rq z różnymi liczbami kwantyli i ustawimy opcję joint na true lub false. Gdy joint jest true: “równość współczynników kierunkowych powinna być wykonana jako wspólne testy na wszystkich parametrach nachylenia”, gdy joint jest false: “należy zgłaszać oddzielne testy na każdym z parametrów nachylenia”.

Zauważ, że testy kwantylowe są testami “linii równoległej”. Oznacza to, że powinniśmy wyjąć różne x-wyrazy_wolne dla każdego kwantyla, ponieważ reprezentują one poziomy rozkładów warunkowych. Jeśli jednak współczynniki kwantyli dla współczynnikow są takie same, to nie ma efektów specyficznych dla kwantyli, wystarczą efekty średnie.

Badanie statystycznej różnicy między 25. i 50. kwantylem warunkowym:

Biorąc pod uwagę powyższe oszacowania kwantyli, różnica między kwantylami 0,25 i 0,50 istnieje, ale czy są one wystarczająco duże, aby być statystycznie różne? Jaka jest wartość p? Przeglądając poniższe wyniki, nie są one statystycznie różne!

Po pierwsze, joint = TRUE. To nie jest testowanie, czy współczynnik na disp jest taki sam jak współczynnik na hp. To jest wspólne testowanie, czy współczynniki dla różnych kwantyli disp i różnych kwantyli hp są takie same dla każdej zmiennej.

kwantyle <- c(0.25, 0.50)
reg_kwantylowa <- rq(mpg ~ disp + hp + factor(am),tau = kwantyle, data = mtcars)
anova(reg_kwantylowa, test = "Wald", joint=TRUE)

## Quantile Regression Analysis of Deviance Table
## 
## Model: mpg ~ disp + hp + factor(am)
## Joint Test of Equality of Slopes: tau in {  0.25 0.5  }
## 
##   Df Resid Df F value Pr(>F)
## 1  3       61  0.8421 0.4761

Po drugie, joint = False:

anova(reg_kwantylowa, test = "Wald", joint=FALSE)

## Quantile Regression Analysis of Deviance Table
## 
## Model: mpg ~ disp + hp + factor(am)
## Tests of Equality of Distinct Slopes: tau in {  0.25 0.5  }
## 
##             Df Resid Df F value Pr(>F)
## disp         1       63  0.0305 0.8619
## hp           1       63  0.5461 0.4627
## factor(am)1  1       63  1.3500 0.2497

Badanie statystycznej różnicy między 25, 50 i 75 kwantylem warunkowym:

Pierwszy kwartyl i mediana nie wydają się być statystycznie różne, teraz dołączymy trzeci kwartyl. Jak widać wcześniej, kwartyle wspólnie wykazują gradient. Teraz możemy zobaczyć, że disp, hp i am są oddzielnie statystycznie różne.

Po pierwsze, joint = TRUE:

kwantyle <- c(0.25, 0.50, 0.75)

reg_kwantylowa <- rq(mpg ~ disp + hp + factor(am),tau = kwantyle, data = mtcars)

anova(reg_kwantylowa, test = "Wald", joint=TRUE)

## Quantile Regression Analysis of Deviance Table
## 
## Model: mpg ~ disp + hp + factor(am)
## Joint Test of Equality of Slopes: tau in {  0.25 0.5 0.75  }
## 
##   Df Resid Df F value   Pr(>F)   
## 1  6       90  3.3173 0.005367 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Po drugie, joint = False:

anova(reg_kwantylowa, test = "Wald", joint=FALSE)

## Quantile Regression Analysis of Deviance Table
## 
## Model: mpg ~ disp + hp + factor(am)
## Tests of Equality of Distinct Slopes: tau in {  0.25 0.5 0.75  }
## 
##             Df Resid Df F value   Pr(>F)   
## disp         2       94  5.4903 0.005558 **
## hp           2       94  6.7221 0.001868 **
## factor(am)1  2       94  7.2758 0.001154 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dobroć dopasowania

Możemy obliczyć współczynniki dobroci dopasowania regresji kwantylowej z wykorzystaniem reszt i reszt bezwarunkowych:

goodfit(resid, resid_nl, tau)

Miara dobroci dopasowania dla regresji kwantylowej jest szacowana jako 1 minus stosunek sumy odchyleń bezwzględnych w modelach w pełni sparametryzowanych do sumy odchyleń bezwzględnych w zerowym (bezwarunkowym) modelu kwantylowym.

Wartości te są przydatne do porównań między modelami kwantylowymi, ale nie są porównywalne ze standardowymi współczynnikami determinacji. Te ostatnie oparte są na wariancji odchyleń kwadratowych, natomiast wartości dobroci dopasowania dla regresji kwantylowej oparte są na odchyleniach bezwzględnych. Wartości dobroci dopasowania zawsze będą mniejsze niż wartości R².

## model kwantylowy
model1 <- rq(mpg ~ disp + hp + factor(am),tau = 0.5, data = mtcars)
reszty1 <- resid(model1)

## bezwarunkowy (pusty) model kwantylowy
model2 <- rq(mpg ~ 1, tau = 0.5,data=mtcars)
reszty2 <- resid(model2)

goodfit(reszty1, reszty2, 0.5)

## [1] 0.5403311

## r2 modelu KMNK dla porównania
model_lm <- lm(mpg ~ disp + hp + factor(am), data = mtcars)

summary(model_lm)$r.squared

## [1] 0.7992061

Zadanie

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.

W przypadku problemów - obejrzyj video tutorial (włącz polskie napisy) oraz wejdź na jego stronę ze źródłami. Możesz również wykorzystać w/w przykłady.

library(PogromcyDanych)
data("CPSSW9298")

dane_1992 <- filter(CPSSW9298, CPSSW9298$year==1992)
dane_1998 <- filter(CPSSW9298, CPSSW9298$year==1998)

Model KMNK dla roku 1992 oraz 1998

Pierwszy model regresji liniowej (KMNK) został stworzony dla roku 1992, gdzie zmienną zależną jest earnings, a predyktorami są: degree (wykształcenie), gender (płeć), oraz age (wiek).Analogiczny model został stworzony dla danych z 1998 roku.

model_92 <- lm(earnings ~ degree + gender + age, data = dane_1992)
summary(model_92)

## 
## Call:
## lm(formula = earnings ~ degree + gender + age, data = dane_1992)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -13.172  -3.324  -0.606   2.580  32.651 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     0.62877    0.61537   1.022    0.307    
## degreebachelor  4.38740    0.11758  37.314   <2e-16 ***
## genderfemale   -2.00554    0.11581 -17.318   <2e-16 ***
## age             0.34221    0.02041  16.771   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.981 on 7586 degrees of freedom
## Multiple R-squared:  0.1956, Adjusted R-squared:  0.1953 
## F-statistic: 614.9 on 3 and 7586 DF,  p-value: < 2.2e-16

model_98 <- lm(earnings ~ degree + gender + age, data = dane_1998)
summary(model_98)

## 
## Call:
## lm(formula = earnings ~ degree + gender + age, data = dane_1998)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -16.448  -4.142  -0.890   2.978  32.109 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      3.3003     0.8729   3.781 0.000158 ***
## degreebachelor   5.3357     0.1643  32.473  < 2e-16 ***
## genderfemale    -2.4932     0.1656 -15.055  < 2e-16 ***
## age              0.3144     0.0290  10.843  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.245 on 5907 degrees of freedom
## Multiple R-squared:  0.1823, Adjusted R-squared:  0.1819 
## F-statistic:   439 on 3 and 5907 DF,  p-value: < 2.2e-16

Analiza wykazała, że wykształcenie , płeć i wiek (“age”) znacząco wpływają na wynagrodzenie w obu badanych latach. Wykształcenie jest związane z wyższymi zarobkami, a jego wpływ wzrósł między 1992 a 1998 rokiem. Luka płacowa między kobietami a mężczyznami utrzymuje się, a w 1998 roku nawet się pogłębiła. Wiek również pozytywnie wpływa na wynagrodzenie, ale jego znaczenie nieco zmalało w 1998 roku. Modele wyjaśniają tylko część zmienności wynagrodzeń, co sugeruje, że inne czynniki również odgrywają istotną rolę.

Wizualizacja danych

Wykorzystano wykres pudełkowy do porównania rozkładu wynagrodzeń w 1992 i 1998 roku.Wykres pozwala zobaczyć różnice w medianie, rozstępie międzykwartylowym oraz ekstremach wynagrodzeń między latami.

library(ggplot2)

ggplot(CPSSW9298, aes(x = factor(year), y = earnings, fill = factor(year))) +
  geom_boxplot(alpha = 0.7) +
  labs(title = "Rozkład wynagrodzeń w latach 1992 i 1998",
       x = "Rok", y = "Wynagrodzenie") +
  theme_minimal()

Model regresji kwantylowej dla 1992

Regresja kwantylowa pokazuje, jak wpływ predyktorów różni się w zależności od poziomu wynagrodzenia. Test WALDA (joint = TRUE) ocenia, czy współczynniki są istotnie różne między kwantylami.

reg_92 <- rq(earnings ~ degree + gender + age, tau = c(0.25, 0.5, 0.75), data = dane_1992)

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

summary(reg_92, se = "boot")

## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = c(0.25, 
##     0.5, 0.75), data = dane_1992)
## 
## tau: [1] 0.25
## 
## Coefficients:
##                Value    Std. Error t value  Pr(>|t|)
## (Intercept)     1.92427  0.75588    2.54572  0.01092
## degreebachelor  3.41880  0.15391   22.21324  0.00000
## genderfemale   -1.14850  0.12298   -9.33889  0.00000
## age             0.18222  0.02518    7.23795  0.00000
## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = c(0.25, 
##     0.5, 0.75), data = dane_1992)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      0.73078   0.76583    0.95423   0.34000
## degreebachelor   4.33269   0.14218   30.47263   0.00000
## genderfemale    -1.91923   0.12204  -15.72591   0.00000
## age              0.31731   0.02674   11.86515   0.00000
## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = c(0.25, 
##     0.5, 0.75), data = dane_1992)
## 
## tau: [1] 0.75
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.24944   0.84683    1.47543   0.14014
## degreebachelor   5.17262   0.18558   27.87204   0.00000
## genderfemale    -2.52525   0.14972  -16.86622   0.00000
## age              0.40965   0.02867   14.28600   0.00000

anova(reg_92, test = "Wald", joint=TRUE)

## Quantile Regression Analysis of Deviance Table
## 
## Model: earnings ~ degree + gender + age
## Joint Test of Equality of Slopes: tau in {  0.25 0.5 0.75  }
## 
##   Df Resid Df F value    Pr(>F)    
## 1  6    22764  39.282 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretacja testu: Wynik testu WALDA wskazuje, że różnice między współczynnikami dla różnych kwantyli są istotne statystycznie.Oznacza to, że wpływ wykształcenia, płci i wieku na wynagrodzenie różni się znacząco w zależności od poziomu wynagrodzeń.

Wpływ rośnie w wyższych kwantylach – osoby z wyższymi wynagrodzeniami czerpią większe korzyści z wykształcenia. Płeć: Luka płacowa na niekorzyść kobiet pogłębia się wraz z kwantylem – kobiety z wyższymi zarobkami tracą relatywnie więcej względem mężczyzn. Wiek: Wiek bardziej wpływa na wyższe wynagrodzenia.

Wizualizacja linii regresji kwantylowej dla 1992.

Linie regresji kwantylowej pokazują, jak zmienia się wpływ wieku na wynagrodzenie w różnych kwantylach.Można zauważyć różne nachylenia linii, co sugeruje, że wpływ wieku na wynagrodzenie nie jest jednorodny w całym rozkładzie.

coeffs_92 <- as.data.frame(t(coef(reg_92)))
coeffs_92$quantile <- c("0.25", "0.50", "0.75")

ggplot(data = dane_1992, aes(x = age, y = earnings, color = gender)) +
  geom_point(alpha = 0.6) +
  geom_abline(data = coeffs_92, aes(intercept = `(Intercept)`, slope = age, color = quantile), size = 1) +
  theme_minimal() +
  labs(
    title = "Regresja kwantylowa dla danych z 1992 roku",
    x = "Wiek",
    y = "Wynagrodzenie",
    color = "Kwantyl"
  )

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Model regresji kwantylowej dla 1998

Analogiczna analiza regresji kwantylowej została wykonana dla danych z 1998 roku.

reg_98 <- rq(earnings ~ degree + gender+ age, tau = c(0.25, 0.5, 0.75), data = dane_1998)

## Warning in rq.fit.br(x, y, tau = tau, ...): Solution may be nonunique

summary(reg_98, se = "boot")

## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = c(0.25, 
##     0.5, 0.75), data = dane_1998)
## 
## tau: [1] 0.25
## 
## Coefficients:
##                Value    Std. Error t value  Pr(>|t|)
## (Intercept)     3.12432  0.88048    3.54843  0.00039
## degreebachelor  3.80303  0.18307   20.77315  0.00000
## genderfemale   -1.64320  0.17011   -9.65942  0.00000
## age             0.18444  0.03042    6.06394  0.00000
## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = c(0.25, 
##     0.5, 0.75), data = dane_1998)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      2.82967   0.92690    3.05281   0.00228
## degreebachelor   4.98626   0.18442   27.03690   0.00000
## genderfemale    -2.40385   0.16475  -14.59068   0.00000
## age              0.30220   0.03121    9.68145   0.00000
## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = c(0.25, 
##     0.5, 0.75), data = dane_1998)
## 
## tau: [1] 0.75
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.44231   1.29400    1.11461   0.26506
## degreebachelor   6.25000   0.32293   19.35414   0.00000
## genderfemale    -3.36538   0.24401  -13.79194   0.00000
## age              0.48077   0.04383   10.96905   0.00000

anova(reg_98, test = "Wald", joint=TRUE)

## Quantile Regression Analysis of Deviance Table
## 
## Model: earnings ~ degree + gender + age
## Joint Test of Equality of Slopes: tau in {  0.25 0.5 0.75  }
## 
##   Df Resid Df F value    Pr(>F)    
## 1  6    17727  30.753 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interpretacja testu: Wynik testu WALDA wskazuje, że różnice między współczynnikami dla różnych kwantyli są istotne statystycznie. Oznacza to, że wpływ wykształcenia, płci i wieku na wynagrodzenie różni się w zależności od poziomu wynagrodzeń.

Wykształcenie: Największy wpływ na wyższe wynagrodzenia (τ=0.75). Płeć: Luka płacowa pogłębia się w miarę wzrostu wynagrodzeń, co wskazuje na większe nierówności płacowe w przypadku wysokich zarobków. Wiek: Wiek bardziej istotny dla osób o wyższych wynagrodzeniach. Ogólne różnice między kwantylami: Regresja kwantylowa uwidacznia istotne różnice w działaniu predyktorów na różnych poziomach wynagrodzeń, czego nie ujawnia regresja liniowa (KMNK).

Wizualizacja linii regresji kwantylowej dla 1998

Linie regresji kwantylowej dla 1998 roku można porównać z wynikami z 1992 roku, aby zobaczyć, czy wpływ wieku i innych predyktorów zmienił się na różnych poziomach wynagrodzeń.

coeffs_98 <- as.data.frame(t(coef(reg_98)))
coeffs_98$quantile <- c("0.25", "0.50", "0.75")

# Wykres
ggplot(data = dane_1998, aes(x = age, y = earnings, color = gender)) +
  geom_point(alpha = 0.6) +
  geom_abline(data = coeffs_98, aes(intercept = `(Intercept)`, slope = age, color = quantile), size = 1) +
  theme_minimal() +
  labs(
    title = "Regresja kwantylowa dla danych z 1998 roku",
    x = "Wiek",
    y = "Wynagrodzenie",
    color = "Kwantyl"
  )