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.65082   15.35395  0.00000
## disp        -0.02441  0.00718   -3.39750  0.00206
## hp          -0.01672  0.01661   -1.00628  0.32290
## factor(am)1  1.39719  1.35580    1.03053  0.31158
## 
## 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.85821   14.79771  0.00000
## disp        -0.02253  0.01574   -1.43150  0.16336
## hp          -0.02713  0.02403   -1.12915  0.26842
## factor(am)1  3.37328  2.10256    1.60437  0.11985
## 
## 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.76373   15.91168  0.00000
## disp         0.00445  0.01580    0.28171  0.78024
## hp          -0.06662  0.02121   -3.14053  0.00396
## factor(am)1  7.91402  2.49269    3.17489  0.00363

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")
dane = CPSSW9298

dane92 <- filter(dane, year == 1992)
dane98 <- filter(dane, year == 1998)

1. Model KMNK

lm_model92 <- lm(earnings ~ degree + gender + age, data = dane92)
summary(lm_model92)

## 
## Call:
## lm(formula = earnings ~ degree + gender + age, data = dane92)
## 
## 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

R-kwadrat dla modelu KNMK wynosi około 0,2, co oznacza, że tylko 20% zmienności w danych jest wyjaśnione przez model. Taki wynik sugeruje, że model ma ograniczoną zdolność do wyjaśniania rozkładu zmiennych zależnych i pozostaje duża część zmienności, która nie jest ujęta w analizie.

2. Regresje kwantylowe dla tau = 0.25, 0.50, 0.75

kwantyle <- c(0.25, 0.50, 0.75)
reg_kwantylowa_92 <- rq(earnings ~ degree + gender + age, tau = kwantyle, data = dane92)

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

summary(reg_kwantylowa_92, se = "boot")

## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = kwantyle, 
##     data = dane92)
## 
## tau: [1] 0.25
## 
## Coefficients:
##                Value    Std. Error t value  Pr(>|t|)
## (Intercept)     1.92427  0.70921    2.71324  0.00668
## degreebachelor  3.41880  0.13755   24.85535  0.00000
## genderfemale   -1.14850  0.12131   -9.46765  0.00000
## age             0.18222  0.02409    7.56283  0.00000
## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = kwantyle, 
##     data = dane92)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      0.73078   0.68726    1.06331   0.28767
## degreebachelor   4.33269   0.12128   35.72555   0.00000
## genderfemale    -1.91923   0.11753  -16.32935   0.00000
## age              0.31731   0.02400   13.21971   0.00000
## 
## Call: rq(formula = earnings ~ degree + gender + age, tau = kwantyle, 
##     data = dane92)
## 
## tau: [1] 0.75
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.24944   0.94223    1.32604   0.18486
## degreebachelor   5.17262   0.17113   30.22693   0.00000
## genderfemale    -2.52525   0.16859  -14.97866   0.00000
## age              0.40965   0.03108   13.18221   0.00000

kwantyle <- c(0.25, 0.50)
reg_kwantylowa <- rq(earnings ~ degree + gender + age, tau = kwantyle, data = dane92)

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

anova(reg_kwantylowa, 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  }
## 
##   Df Resid Df F value    Pr(>F)    
## 1  3    15177  48.241 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wyniki analizy wskazują, że zmienne takie jak wykształcenie, płeć i wiek mają istotny wpływ na dochody, szczególnie w różnych kwantylach dochodów, takich jak 0.25 i 0.5. Oznacza to, że te czynniki wywierają różny wpływ na dochody w zależności od poziomu dochodu, co sugeruje, że efekt tych zmiennych może się zmieniać w różnych częściach rozkładu dochodów. Ponadto, wartość p uzyskana w badaniu jest bardzo mała, co pozwala na odrzucenie hipotezy zerowej o braku różnic w nachyleniach między kwantylami. W rezultacie, możemy uznać, że istnieje statystycznie istotna różnica w sposobie, w jaki model zachowuje się w różnych kwantylach, co oznacza, że zmienne te mają zróżnicowany wpływ na dochody w różnych częściach rozkładu.

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

## Quantile Regression Analysis of Deviance Table
## 
## Model: earnings ~ degree + gender + age
## Tests of Equality of Distinct Slopes: tau in {  0.25 0.5  }
## 
##                Df Resid Df F value    Pr(>F)    
## degreebachelor  1    15179  54.171 1.932e-13 ***
## genderfemale    1    15179  49.118 2.512e-12 ***
## age             1    15179  47.431 5.920e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Te zmienne mają istotny wpływ na dochody w różnych częściach rozkładu, co sugeruje, że wpływ wykształcenia, płci i wieku na dochody może różnić się w zależności od tego, czy analizujemy niższe, czy wyższe kwantyle dochodów. W praktyce oznacza to, że te czynniki mogą oddziaływać na osoby z niższymi dochodami w inny sposób niż na osoby z wyższymi dochodami. Takie różnice w wpływie mogą wynikać z tego, że na różnych poziomach dochodów różne zmienne mogą odgrywać odmienną rolę w kształtowaniu wynagrodzeń. Na przykład, wykształcenie może mieć większy wpływ na osoby w niższych kwantylach dochodów, podczas gdy inne czynniki, takie jak doświadczenie zawodowe czy sieć kontaktów, mogą być bardziej znaczące dla osób w wyższych kwantylach. Takie różnice wskazują, że podejście do analizy dochodów powinno uwzględniać nie tylko ogólny wpływ zmiennych, ale także ich zróżnicowane efekty w zależności od poziomu dochodu.

kwantyle <- c(0.25, 0.50, 0.75)

reg_kwantylowa <- rq(earnings ~ degree + gender + age, tau = kwantyle, data = dane92)

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

anova(reg_kwantylowa, 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

Tak jak wczesniej takze i przy uwzględnieniu kwantylu 0.75 wszystkie zmienne okazały się być istotne statystyczne.

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

## Quantile Regression Analysis of Deviance Table
## 
## Model: earnings ~ degree + gender + age
## Tests of Equality of Distinct Slopes: tau in {  0.25 0.5 0.75  }
## 
##                Df Resid Df F value    Pr(>F)    
## degreebachelor  2    22768  48.950 < 2.2e-16 ***
## genderfemale    2    22768  39.262 < 2.2e-16 ***
## age             2    22768  35.296 4.441e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wyniki te są również podobne do tych uzyskanych przy analizie kwantyli, w której nie uwzględniono wartości 0.75. Oznacza to, że wpływ zmiennych takich jak wykształcenie, płeć i wiek na dochody pozostaje zasadniczo spójny, niezależnie od tego, czy analizujemy tylko kwantyle 0.25 i 0.5, czy również dodajemy 0.75. Choć uwzględnienie dodatkowego kwantyla mogłoby potencjalnie wpłynąć na wyniki, obserwowane zmiany w dochodach w niższych i średnich kwantylach wskazują na podobne zależności, co sugeruje, że główne efekty zmiennych pozostają stabilne.

## model kwantylowy dla 0.25
model1 <- rq(earnings ~ degree + gender + age, tau = 0.25, data = dane92)
reszty1 <- resid(model1)

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

goodfit(reszty1, reszty2, 0.25)

## [1] 0.08036852

## model kwantylowy dla 0.5
model1 <- rq(earnings ~ degree + gender + age, tau = 0.5, data = dane92)

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

reszty1 <- resid(model1)

## bezwarunkowy (pusty) model kwantylowy
model2 <- rq(earnings ~ 1, tau = 0.5,data=dane92)

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

reszty2 <- resid(model2)

goodfit(reszty1, reszty2, 0.5)

## [1] 0.1095095

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

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

goodfit(reszty1, reszty2, 0.75)

## [1] 0.1401108

Wartości w zakresie od 0.08 do 0.14 wskazują, że model kwantylowy rzeczywiście wprowadza pewną poprawę w stosunku do modelu pustego, jednak ta zmiana nie jest znacząca. Takie wartości sugerują, że wpływ zmiennych na dochody w analizowanych kwantylach nie jest bardzo silny, a różnica w dopasowaniu modeli nie jest ogromna. W szczególności, gdyby wartości miary dopasowania były bliskie zeru, oznaczałoby to, że zmienne miałyby bardzo silny wpływ na dochody w danym kwantylu.

Biorąc pod uwagę miarę dopasowania, model kwantylowy dla tau = 0.25 wydaje się być najbardziej odpowiedni, ponieważ oferuje najlepsze dopasowanie do danych w tym kwantylu. Oznacza to, że zmienne takie jak wykształcenie, płeć i wiek mają najsilniejszy wpływ na dochody w tym właśnie przedziale dochodów w porównaniu do innych modeli, co sugeruje, że w tym zakresie zmienne te odgrywają kluczową rolę w kształtowaniu wynagrodzeń. Choć zmiana w dopasowaniu nie jest drastyczna, model kwantylowy dla tau = 0.25 wydaje się lepiej oddawać rzeczywisty rozkład dochodów w tym konkretnym kwantylu.