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.51562   16.72360  0.00000
## disp        -0.02441  0.00853   -2.86281  0.00786
## hp          -0.01672  0.01542   -1.08383  0.28768
## factor(am)1  1.39719  1.38280    1.01041  0.32095
## 
## 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.65388   16.62588  0.00000
## disp        -0.02253  0.01620   -1.39028  0.17540
## hp          -0.02713  0.02352   -1.15343  0.25849
## factor(am)1  3.37328  1.97624    1.70692  0.09891
## 
## 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.94463   14.43146  0.00000
## disp         0.00445  0.01575    0.28265  0.77953
## hp          -0.06662  0.02054   -3.24387  0.00305
## factor(am)1  7.91402  2.47636    3.19582  0.00344

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

Teraz Wasza kolej ;-)

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.

data("CPSSW9298")
# ?CPSSW9298 

Dane 1992

dane<-CPSSW9298
dane_1992<-dane %>% filter(year==1992)
dane_1992 %>%
  ggplot(aes(age,earnings))+
  geom_point()

data(dane_1992) #dane 

## Warning in data(dane_1992): zbiór danych 'dane_1992' nie został znaleziony

p <- ggplot(data = dane_1992) +
    geom_point(mapping = aes(x = age, y = earnings), color = "blue")
taus <- c(0.20, 0.40, 0.60, 0.80, 0.95, 0.999)
fits <- data.frame(
    coef(lm(earnings ~ age, data = dane_1992)),
    sapply(taus, function(x) coef(rq(formula = earnings ~ age, data = dane_1992, tau = x))))
names(fits) <- c("OLS", sprintf("$\\tau_{%0.2f}$", taus))
nf <- ncol(fits)
colors <- colorRampPalette(colors = c("black", "red"))(nf)
p <- p + geom_abline(intercept = fits[1, 1], slope = fits[2, 1], color = colors[1], linewidth = 1.5)
for (i in seq_len(nf)[-1]) {
    p <- p + geom_abline(intercept = fits[1, i], slope = fits[2, i], color = colors[i])
}
p

data(dane_1992) #dane 

## Warning in data(dane_1992): zbiór danych 'dane_1992' nie został znaleziony

p <- ggplot(data = dane_1992) +
    geom_point(mapping = aes(x = gender, y = earnings), color = "blue")
taus <- c(0.20, 0.40, 0.60, 0.80, 0.95, 0.999)
fits <- data.frame(
    coef(lm(earnings ~ gender, data = dane_1992)),
    sapply(taus, function(x) coef(rq(formula = earnings ~ gender, data = dane_1992, tau = x))))
names(fits) <- c("OLS", sprintf("$\\tau_{%0.2f}$", taus))
nf <- ncol(fits)
colors <- colorRampPalette(colors = c("black", "red"))(nf)
p <- p + geom_abline(intercept = fits[1, 1], slope = fits[2, 1], color = colors[1], linewidth = 1.5)
for (i in seq_len(nf)[-1]) {
    p <- p + geom_abline(intercept = fits[1, i], slope = fits[2, i], color = colors[i])
}
p

q20 <- rq(earnings ~ age+degree+gender, data = dane_1992, tau = 0.20)
q40 <- rq(earnings ~ age+degree+gender, data = dane_1992, tau = 0.40)
q60 <- rq(earnings ~ age+degree+gender, data = dane_1992, tau = 0.60)
q80 <- rq(earnings ~ age+degree+gender, data = dane_1992, tau = 0.80)
q95 <- rq(earnings ~ age+degree+gender, data = dane_1992, tau = 0.95)
q999 <- rq(earnings ~ age+degree+gender, data = dane_1992, tau = 0.999)
# Tabela z porównaniem wyników trzech modeli: 

stargazer(q20, q40, q60, q80, q95, q999, title = "Wyniki regresji kwantylowych", type = "text")

## 
## Wyniki regresji kwantylowych
## =========================================================================
##                                   Dependent variable:                    
##                ----------------------------------------------------------
##                                         earnings                         
##                   (1)       (2)       (3)       (4)       (5)      (6)   
## -------------------------------------------------------------------------
## age            0.151***  0.256***  0.374***  0.451***  0.593***   0.901  
##                 (0.020)   (0.022)   (0.024)   (0.028)   (0.059)  (0.628) 
##                                                                          
## degreebachelor 3.136***  4.008***  4.690***  5.502***  6.506***   4.687  
##                 (0.135)   (0.144)   (0.154)   (0.189)   (0.362)  (3.598) 
##                                                                          
## genderfemale   -0.962*** -1.621*** -2.233*** -2.794*** -4.071***  -2.284 
##                 (0.112)   (0.124)   (0.138)   (0.158)   (0.335)  (3.555) 
##                                                                          
## Constant       2.225***   1.383**    0.214     0.918     1.923    7.812  
##                 (0.597)   (0.650)   (0.718)   (0.822)   (1.777)  (18.586)
##                                                                          
## -------------------------------------------------------------------------
## Observations     7,590     7,590     7,590     7,590     7,590    7,590  
## =========================================================================
## Note:                                         *p<0.1; **p<0.05; ***p<0.01

Interpretacja:

-Kobiety zarabiają średnio o 0,962 dolara za godzinę mniej niż mężczyźni (na podstawie pierwszego kwartylu).

-W wyższych dochodach płeć oraz wykształcenie mają jeszcze silniejszy wpływ na wysokość wynagrodzenia.

-Wyraz wolny w górnych kwartylach jest statystycznie nieistotny, co sugeruje, że warto rozważyć transformację dochodów, na przykład przez ich logarytmowanie, aby lepiej uchwycić zależności.

q20 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.20)
q40 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.40)
q60 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.60)
q80 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.80)
q95 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.95)
q999 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.999)
# Tabela z porównaniem wyników trzech modeli: 

stargazer(q20, q40, q60, q80, q95, q999, title = "Wyniki regresji kwantylowych", type = "text")

## 
## Wyniki regresji kwantylowych
## =========================================================================
##                                   Dependent variable:                    
##                ----------------------------------------------------------
##                                      log(earnings)                       
##                   (1)       (2)       (3)       (4)       (5)      (6)   
## -------------------------------------------------------------------------
## age            0.020***  0.027***  0.031***  0.032***  0.028***   0.024  
##                 (0.003)   (0.002)   (0.002)   (0.002)   (0.003)  (0.017) 
##                                                                          
## degreebachelor 0.402***  0.403***  0.376***  0.364***  0.323***   0.135  
##                 (0.017)   (0.013)   (0.012)   (0.012)   (0.017)  (0.098) 
##                                                                          
## genderfemale   -0.141*** -0.174*** -0.187*** -0.201*** -0.197***  -0.061 
##                 (0.016)   (0.013)   (0.012)   (0.012)   (0.017)  (0.094) 
##                                                                          
## Constant       1.300***  1.392***  1.484***  1.696***  2.121***  2.815***
##                 (0.086)   (0.070)   (0.062)   (0.065)   (0.093)  (0.504) 
##                                                                          
## -------------------------------------------------------------------------
## Observations     7,590     7,590     7,590     7,590     7,590    7,590  
## =========================================================================
## Note:                                         *p<0.1; **p<0.05; ***p<0.01

Wyraz wolny jest teraz istotny. Kwartyl na poziomie 0,999 okazał się nieistotny, dlatego decydujemy się na jego wykluczenie.

q20 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.20)
q40 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.40)
q60 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.60)
q80 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.80)
q95 <- rq(log(earnings) ~ age+degree+gender, data = dane_1992, tau = 0.95)

# Tabela z porównaniem wyników trzech modeli: 

stargazer(q20, q40, q60, q80, q95, title = "Wyniki regresji kwantylowych", type = "text")

## 
## Wyniki regresji kwantylowych
## ================================================================
##                               Dependent variable:               
##                -------------------------------------------------
##                                  log(earnings)                  
##                   (1)       (2)       (3)       (4)       (5)   
## ----------------------------------------------------------------
## age            0.020***  0.027***  0.031***  0.032***  0.028*** 
##                 (0.003)   (0.002)   (0.002)   (0.002)   (0.003) 
##                                                                 
## degreebachelor 0.402***  0.403***  0.376***  0.364***  0.323*** 
##                 (0.017)   (0.013)   (0.012)   (0.012)   (0.017) 
##                                                                 
## genderfemale   -0.141*** -0.174*** -0.187*** -0.201*** -0.197***
##                 (0.016)   (0.013)   (0.012)   (0.012)   (0.017) 
##                                                                 
## Constant       1.300***  1.392***  1.484***  1.696***  2.121*** 
##                 (0.086)   (0.070)   (0.062)   (0.065)   (0.093) 
##                                                                 
## ----------------------------------------------------------------
## Observations     7,590     7,590     7,590     7,590     7,590  
## ================================================================
## Note:                                *p<0.1; **p<0.05; ***p<0.01

Model KMNK

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

## 
## Call:
## lm(formula = log(earnings) ~ age + degree + gender, data = dane_1992)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.80573 -0.26371  0.02737  0.29516  1.53710 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.486349   0.053186   27.95   <2e-16 ***
## age             0.026398   0.001764   14.97   <2e-16 ***
## degreebachelor  0.375044   0.010162   36.91   <2e-16 ***
## genderfemale   -0.167263   0.010009  -16.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4305 on 7586 degrees of freedom
## Multiple R-squared:  0.1872, Adjusted R-squared:  0.1869 
## F-statistic: 582.3 on 3 and 7586 DF,  p-value: < 2.2e-16

Model wykazuje niski współczynnik determinacji R^2, co wskazuje na słabe dopasowanie modelu do danych.

kwantyle <- c(0.20, 0.40, 0.60, 0.80, 0.95)
reg_kwantylowa <- rq(log(earnings) ~ age+degree+gender,tau = kwantyle, data = dane_1992)
summary(reg_kwantylowa, se = "boot")

## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1992)
## 
## tau: [1] 0.2
## 
## Coefficients:
##                Value    Std. Error t value  Pr(>|t|)
## (Intercept)     1.30031  0.07573   17.16936  0.00000
## age             0.02018  0.00248    8.14165  0.00000
## degreebachelor  0.40238  0.01670   24.09849  0.00000
## genderfemale   -0.14129  0.01436   -9.83596  0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1992)
## 
## tau: [1] 0.4
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.39213   0.06268   22.21106   0.00000
## age              0.02659   0.00202   13.17524   0.00000
## degreebachelor   0.40268   0.01080   37.27034   0.00000
## genderfemale    -0.17407   0.01149  -15.14400   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1992)
## 
## tau: [1] 0.6
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.48369   0.06263   23.68889   0.00000
## age              0.03119   0.00207   15.04158   0.00000
## degreebachelor   0.37644   0.01242   30.29884   0.00000
## genderfemale    -0.18712   0.01233  -15.17826   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1992)
## 
## tau: [1] 0.8
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.69618   0.06778   25.02535   0.00000
## age              0.03201   0.00223   14.36281   0.00000
## degreebachelor   0.36352   0.01081   33.62071   0.00000
## genderfemale    -0.20070   0.01203  -16.68646   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1992)
## 
## tau: [1] 0.95
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      2.12112   0.09912   21.39921   0.00000
## age              0.02807   0.00323    8.67976   0.00000
## degreebachelor   0.32269   0.02067   15.61489   0.00000
## genderfemale    -0.19652   0.01916  -10.25558   0.00000

kwantyle <- c(0.20, 0.40, 0.60, 0.80, 0.95)
reg_kwantylowa <- rq(log(earnings) ~ age+degree+gender,tau = kwantyle, data = dane_1992)
anova(reg_kwantylowa, test = "Wald", joint=TRUE)

## Quantile Regression Analysis of Deviance Table
## 
## Model: log(earnings) ~ age + degree + gender
## Joint Test of Equality of Slopes: tau in {  0.2 0.4 0.6 0.8 0.95  }
## 
##   Df Resid Df F value    Pr(>F)    
## 1 12    37938  4.1132 1.829e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wartość p (1.829e-06) jest bardzo mała, co wskazuje na istotność statystyczną. Oznacza to, że istnieje statystycznie istotna różnica w współczynnikach nachylenia dla różnych kwantyli, tzn. współczynniki te nie są równe na poziomach kwantyli 0.2, 0.4, 0.6, 0.8 i 0.95.

Model kwantylowy

## model kwantylowy
model1 <- rq(log(earnings) ~ age+degree+gender, tau = 0.80, data = dane_1992)
reszty1 <- resid(model1)

## bezwarunkowy (pusty) model kwantylowy
model2 <- rq(log(earnings) ~ 1, tau = 0.80, data = dane_1992)
reszty2 <- resid(model2)

goodfit(reszty1, reszty2, 0.80)

## [1] 0.1256738

## r2 modelu KMNK dla porównania
model_lm <- lm(log(earnings) ~ age+degree+gender, data = dane_1992)

summary(model_lm)$r.squared

## [1] 0.1871756

Wartość 0.1256738 modelu kwantylowego jest relatywnie niska, co sugeruje, że dopasowanie modelu kwantylowego nie jest idealne, ale nadal może być użyteczne. Niższa wartość wskazuje na to, że model kwantylowy nie poprawił się znacząco w porównaniu do modelu pustego. Warto jednak zaznaczyć, że model kwantylowy lepiej radzi sobie z nierówną zmiennością w różnych częściach danych, a model liniowy zapewnia tylko jedno ogólne dopasowanie.

Dane 1998

dane_1998<-dane %>% filter(year==1998)
dane_1998 %>%
  ggplot(aes(age,earnings))+
  geom_point()

data(dane_1998) #dane 

## Warning in data(dane_1998): zbiór danych 'dane_1998' nie został znaleziony

p <- ggplot(data = dane_1998) +
    geom_point(mapping = aes(x = age, y = earnings), color = "blue")
taus <- c(0.20, 0.40, 0.60, 0.80, 0.90, 0.95, 0.99)
fits <- data.frame(
    coef(lm(earnings ~ age, data = dane_1992)),
    sapply(taus, function(x) coef(rq(formula = earnings ~ age, data = dane_1998, tau = x))))

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

names(fits) <- c("OLS", sprintf("$\\tau_{%0.2f}$", taus))
nf <- ncol(fits)
colors <- colorRampPalette(colors = c("black", "red"))(nf)
p <- p + geom_abline(intercept = fits[1, 1], slope = fits[2, 1], color = colors[1], linewidth = 1.5)
for (i in seq_len(nf)[-1]) {
    p <- p + geom_abline(intercept = fits[1, i], slope = fits[2, i], color = colors[i])
}
p

Tutaj zdecydowaliśmy się na dodanie kwantyla 0.99, aby uchwycić więcej informacji z danych i lepiej zrozumieć zależności w najwyższych poziomach rozkładu dochodów.

q20_98 <- rq(earnings ~ age+degree+gender, data = dane_1998, tau = 0.20)

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

q40_98 <- rq(earnings ~ age+degree+gender, data = dane_1998, tau = 0.40)
q60_98 <- rq(earnings ~ age+degree+gender, data = dane_1998, tau = 0.60)
q80_98 <- rq(earnings ~ age+degree+gender, data = dane_1998, tau = 0.80)
q90_98 <- rq(earnings ~ age+degree+gender, data = dane_1998, tau = 0.90)
q95_98 <- rq(earnings ~ age+degree+gender, data = dane_1998, tau = 0.95)
q99_98 <- rq(earnings ~ age+degree+gender, data = dane_1998, tau = 0.99)
# Tabela z porównaniem wyników trzech modeli: 

stargazer(q20_98, q40_98, q60_98, q80_98, q90_98, q95_98,q99_98, title = "Wyniki regresji kwantylowych", type = "text")

## 
## Wyniki regresji kwantylowych
## ====================================================================================
##                                         Dependent variable:                         
##                ---------------------------------------------------------------------
##                                              earnings                               
##                   (1)       (2)       (3)       (4)       (5)       (6)       (7)   
## ------------------------------------------------------------------------------------
## age            0.160***  0.253***  0.361***  0.454***  0.481***  0.466***    0.310  
##                 (0.024)   (0.029)   (0.033)   (0.046)   (0.069)   (0.087)   (0.208) 
##                                                                                     
## degreebachelor 3.526***  4.432***  5.311***  6.867***  8.761***  10.096*** 12.392***
##                 (0.142)   (0.171)   (0.192)   (0.281)   (0.419)   (0.527)   (1.104) 
##                                                                                     
## genderfemale   -1.442*** -1.976*** -2.671*** -3.633*** -4.434*** -4.677*** -5.706***
##                 (0.133)   (0.162)   (0.183)   (0.258)   (0.388)   (0.485)   (1.105) 
##                                                                                     
## Constant       3.205***  2.972***   2.408**   3.351**  5.769***  9.101***  20.755***
##                 (0.693)   (0.846)   (0.958)   (1.366)   (2.032)   (2.665)   (6.516) 
##                                                                                     
## ------------------------------------------------------------------------------------
## Observations     5,911     5,911     5,911     5,911     5,911     5,911     5,911  
## ====================================================================================
## Note:                                                    *p<0.1; **p<0.05; ***p<0.01

Wiek (age):

Wpływ wieku na wynagrodzenie rośnie w wyższych kwantylach dochodów. Na przykład, dla 20. percentyla wzrost wieku o jedną jednostkę zwiększa dochód średnio o 0.160, podczas gdy dla 80. percentyla ten wzrost wynosi już 0.454. W górnych kwantylach (np. 95. percentyl) efekt wieku pozostaje istotny statystycznie, ale nie zmienia się znacząco w stosunku do 90. percentyla. Jednak dla 99. percentyla efekt jest słabszy i staje się nieistotny.

Wykształcenie (degreebachelor):

Efekt posiadania stopnia licencjata jest wyraźny i rośnie w miarę wzrostu kwantyla. Na przykład, dla 20. percentyla wykształcenie zwiększa dochód średnio o 3.526 jednostek, podczas gdy dla 95. percentyla efekt wynosi już 10.096 jednostek, a dla 99. percentyla aż 12.392 jednostek. Wskazuje to, że w wyższych dochodach wykształcenie odgrywa coraz większą rolę, co sugeruje, że wykształcenie jest kluczowym czynnikiem w uzyskaniu wysokich dochodów.

Płeć (genderfemale):

Efekt płci (bycia kobietą) jest negatywny i wzrasta wraz z kwantylem, co oznacza, że różnica w wynagrodzeniu między kobietami a mężczyznami staje się większa w wyższych dochodach. Dla 20. percentyla kobiety zarabiają średnio o 1.442 jednostki mniej niż mężczyźni, ale w 95. percentylu różnica ta wynosi już 4.677 jednostek, a w 99. percentylu (model 7) aż 5.706 jednostek. Sugeruje to, że luka płacowa między kobietami a mężczyznami jest szczególnie wyraźna w najwyższych poziomach dochodów. Co istotne, luka płciowa (różnice w dochodach między kobietami a mężczyznami) jest bardziej widoczna w danych z 1998 roku niż w danych z 1992 roku, co może wskazywać na nasilające się nierówności w wynagrodzeniach w tym okresie.

Model KMNK

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

## 
## Call:
## lm(formula = log(earnings) ~ age + degree + gender, data = dane_1998)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.96583 -0.27644  0.02536  0.30209  1.50215 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     1.78845    0.06230   28.71   <2e-16 ***
## age             0.02142    0.00207   10.35   <2e-16 ***
## degreebachelor  0.38277    0.01173   32.64   <2e-16 ***
## genderfemale   -0.18003    0.01182  -15.23   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4457 on 5907 degrees of freedom
## Multiple R-squared:  0.1828, Adjusted R-squared:  0.1824 
## F-statistic: 440.5 on 3 and 5907 DF,  p-value: < 2.2e-16

Powyższy model, podobnie jak w przypadku modelu z danymi z 1992 roku, wykazuje niski poziom R^2, zatem model słabo dopasowuje się do danych.

kwantyle_1998 <- c(0.20, 0.40, 0.60, 0.80, 0.90, 0.95, 0.99)
reg_kwantylowa_1998 <- rq(log(earnings) ~ age+degree+gender,tau = kwantyle, data = dane_1998)

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

summary(reg_kwantylowa_1998, se = "boot")

## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1998)
## 
## tau: [1] 0.2
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.50442   0.08728   17.23633   0.00000
## age              0.01843   0.00301    6.12589   0.00000
## degreebachelor   0.41255   0.01453   28.38607   0.00000
## genderfemale    -0.16917   0.01648  -10.26764   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1998)
## 
## tau: [1] 0.4
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.66252   0.07218   23.03171   0.00000
## age              0.02255   0.00241    9.35761   0.00000
## degreebachelor   0.39379   0.01338   29.42057   0.00000
## genderfemale    -0.17863   0.01224  -14.59643   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1998)
## 
## tau: [1] 0.6
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.78327   0.06705   26.59740   0.00000
## age              0.02605   0.00225   11.59740   0.00000
## degreebachelor   0.37942   0.01340   28.31883   0.00000
## genderfemale    -0.19710   0.01360  -14.48759   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1998)
## 
## tau: [1] 0.8
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      2.04254   0.07741   26.38640   0.00000
## age              0.02569   0.00253   10.16456   0.00000
## degreebachelor   0.37692   0.01461   25.80463   0.00000
## genderfemale    -0.21061   0.01590  -13.24672   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = dane_1998)
## 
## tau: [1] 0.95
## 
## Coefficients:
##                Value    Std. Error t value  Pr(>|t|)
## (Intercept)     2.55322  0.14726   17.33780  0.00000
## age             0.01898  0.00472    4.02407  0.00006
## degreebachelor  0.40561  0.02137   18.98149  0.00000
## genderfemale   -0.19263  0.02062   -9.34298  0.00000

reg_kwantylowa_1998 <- rq(log(earnings) ~ age+degree+gender,tau = kwantyle_1998, data = dane_1998)

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

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

## Quantile Regression Analysis of Deviance Table
## 
## Model: log(earnings) ~ age + degree + gender
## Joint Test of Equality of Slopes: tau in {  0.2 0.4 0.6 0.8 0.9 0.95 0.99  }
## 
##   Df Resid Df F value   Pr(>F)   
## 1 18    41359  2.0464 0.005521 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wartość p jest mniejsza niż 0.01, co wskazuje na istotność statystyczną na poziomie 1%. Oznacza to, że istnieje statystycznie istotna różnica w nachyleniach (współczynnikach) dla różnych kwantyli w zakresie {0.2,0.4,0.6,0.8,0.9,0.95,0.99}.

Model kwantylowy

## model kwantylowy
model1_1998 <- rq(log(earnings) ~ age+degree+gender, tau = 0.70, data = dane_1998)
reszty1_1998 <- resid(model1_1998)

## bezwarunkowy (pusty) model kwantylowy
model2_1998 <- rq(log(earnings) ~ 1, tau = 0.70, data = dane_1998)
reszty2_1998 <- resid(model2_1998)

goodfit(reszty1_1998, reszty2_1998, 0.70)

## [1] 0.1116213

## r2 modelu KMNK dla porównania
model_lm_1998 <- lm(log(earnings) ~ age+degree+gender, data = dane_1998)

summary(model_lm_1998)$r.squared

## [1] 0.1828058

Wartość 0.1116213 dla miary goodfit w modelu kwantylowym (tau = 0.70) jest relatywnie niska. Oznacza to, że dopasowanie modelu kwantylowego w porównaniu do modelu pustego (bez zmiennych objaśniających) nie jest idealne, ale nadal może być uznane za użyteczne.

Niższa wartość goodfit sugeruje, że model kwantylowy nie wprowadza dużego wyjaśnienia zmienności w danych w porównaniu do modelu pustego. Jednak istotą modelu kwantylowego jest to, że koncentruje się on na specyficznej części rozkładu, co czyni go bardziej odpowiednim do analizy rozkładu dochodów niż klasyczny model liniowy.

Wartość R^2 wynosi 0.1828058, co oznacza, że model liniowy wyjaśnia około 18.28% zmienności logarytmów dochodów w danych. Mimo że R^2 jest wyższe niż wartość goodfit w modelu kwantylowym, model liniowy zakłada jednolitą zależność między zmiennymi objaśniającymi a zmienną zależną w całym rozkładzie danych, co może nie odzwierciedlać różnic w różnych częściach rozkładu.

Wnioski:

Model kwantylowy może być bardziej użyteczny, gdy interesuje nas analiza specyficznej części rozkładu, ponieważ pokazuje, jak zmienne objaśniające wpływają na dochody w tym konkretnym kwantylu. Model liniowy dostarcza bardziej ogólnej analizy, wyjaśniając globalną zmienność w danych, ale może nie oddawać precyzyjnie różnic w zależnościach między zmiennymi w różnych częściach rozkładu dochodów.