Raport 4.

Regresja kwantylowa

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.

Na początku należy zaimportować dane i przedstawić ich rozkład na wykresie:

data("CPSSW9298")

wynagrodzenia<-CPSSW9298


a <- ggplot(wynagrodzenia, aes(x=earnings)) +
  
           geom_histogram(bins=20) +
  
           ggtitle("Rozkład danych earnings") +
  
           xlab("Wynagrodzenia") + ylab("Ilość")
plot(a)

Dane nie są równomiernie rozłożone, zatem można przypuszczać, że więcej wynagrodzeń jest mniejsza niż średnie wynagrodzenie.

Płace vs płeć

b <- ggplot(wynagrodzenia,
            
           aes(x=gender, y=earnings)) + 
  
           geom_boxplot() +
  
           ggtitle("Płace wg płci") +
  
           ylab("Wynagrodzenia") +
  
           labs(fill="Płeć") 
plot(b)

Widać, że kobiety zarabiają mniej niż mężczyźni, ponieważ dane na wykresie wskazują wyższe wartości dla mężczyzn.

Metoda KMNK:

kmnk <- lm(earnings ~ age + gender , data = CPSSW9298)

summary(kmnk)

## 
## Call:
## lm(formula = earnings ~ age + gender, data = CPSSW9298)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.446  -4.347  -1.020   3.072  36.507 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    4.5092     0.5664   7.961 1.85e-15 ***
## age            0.3012     0.0189  15.941  < 2e-16 ***
## genderfemale  -1.8909     0.1073 -17.619  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.156 on 13498 degrees of freedom
## Multiple R-squared:  0.04077,    Adjusted R-squared:  0.04062 
## F-statistic: 286.8 on 2 and 13498 DF,  p-value: < 2.2e-16

Z procedury KMNK widać, że na zarobki wpływa wiek i płeć.

Regresja kwantylowa rózno-poziomowa

kwantyle <- c(0.25, 0.50, 0.75)

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

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

summary(reg_kwantylowa , se = "boot")

## 
## Call: rq(formula = earnings ~ age + gender, tau = kwantyle, data = wynagrodzenia)
## 
## tau: [1] 0.25
## 
## Coefficients:
##              Value     Std. Error t value   Pr(>|t|) 
## (Intercept)    4.06469   0.57355    7.08686   0.00000
## age            0.16026   0.01939    8.26381   0.00000
## genderfemale  -1.34033   0.10620  -12.62107   0.00000
## 
## Call: rq(formula = earnings ~ age + gender, tau = kwantyle, data = wynagrodzenia)
## 
## tau: [1] 0.5
## 
## Coefficients:
##              Value     Std. Error t value   Pr(>|t|) 
## (Intercept)    4.24680   0.64262    6.60862   0.00000
## age            0.27244   0.02155   12.64074   0.00000
## genderfemale  -1.71474   0.11942  -14.35953   0.00000
## 
## Call: rq(formula = earnings ~ age + gender, tau = kwantyle, data = wynagrodzenia)
## 
## tau: [1] 0.75
## 
## Coefficients:
##              Value     Std. Error t value   Pr(>|t|) 
## (Intercept)    5.04127   1.15553    4.36273   0.00001
## age            0.39664   0.03834   10.34513   0.00000
## genderfemale  -2.47197   0.15027  -16.45043   0.00000

Jak widać, wszystkie zmienne wykazują statystyczną istotność dla każdego kwantyla (0,25, 0,5 oraz 0,75).

Badanie statystycznej różnicy

Badanie statystycznej różnicy między 25, 50 i 75 kwantylem warunkowym procedurą ANOVA (sprawdzamy istnienie różnic)

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

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

Wszystkie trzy kwartyle różnią się, co widzimy po otrzymanych wynikach.

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

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

Zmienne objaśniające są oddzielnie różne statystycznie.

Dobroć dopasowania

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

reszty1 <- resid(kmnk)

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

goodfit(reszty1, reszty2, 0.5)

## [1] 0.002287609

## bezwarunkowy (pusty) model kwantylowy
model3 <- rq( earnings~1, tau = 0.75,data=wynagrodzenia)
reszty3 <- resid(model3)

goodfit(reszty1, reszty3, 0.5)

## [1] 0.207917

## r2 modelu KMNK dla porównania
summary(kmnk)$r.squared

## [1] 0.04076514

Najwyższa wartość jest dla współczynnika dobroci dopasowania jest dla modelu trzeciego, czyli modelu kwantylowego kwantyla=0,75. Następnie najlepszy współczynnik dała metoda KMNK, a następnie regresja kwantylowa drugiego kwantyla.