Nieklasyczne metody statystyki

Regresja kwantylowa

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.

data("CPSSW9298")
# ?CPSSW9298 

data92<-CPSSW9298 %>% filter(CPSSW9298$year==1998)
  data92%>%
  ggplot(aes(age,earnings))+
    geom_point()

data(data92) #dane 
p <- ggplot(data = data92) +
    geom_point(mapping = aes(x = age, y = earnings), color = "blue")
taus <- c(0.1, 0.25, 0.5, 0.75, 0.90, 0.95)
fits <- data.frame(
    coef(lm(earnings ~ age, data = data92)),
    sapply(taus, function(x) coef(rq(formula = earnings ~ age, data = data92, 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

Wykres przedstawiający kształtowanie się zarobków w stosunku do wieku oraz regresje kwantylowe dla kwantyli 0.1, 0.25 , 0.5, 0.75, 0.90

data(data92) #dane 
p <- ggplot(data = data92) +
    geom_point(mapping = aes(x = gender, y = earnings), color = "blue")
taus <- c(0.1, 0.25, 0.5, 0.75, 0.90, 0.95)
fits <- data.frame(
    coef(lm(earnings ~ gender, data = data92)),
    sapply(taus, function(x) coef(rq(formula = earnings ~ gender, data = data92, 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

Wykres przedstawiający kształtowanie się zarobków dla mężczyzn i kobiet oraz przedstawione regresje dla kwantyli 0.1, 0.25 , 0.5, 0.75, 0.90

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)	14.805467	7.019231	9.615385	13.461538	18.376068	24.038462	28.846153
genderfemale	-2.044285	-1.019231	-1.442308	-1.923077	-2.991453	-2.884615	-3.846153

q2 <- rq(log(earnings) ~ age+degree+gender, data = data92, tau = 0.25)
q4 <- rq(log(earnings) ~ age+degree+gender, data = data92, tau = 0.50)
q6 <- rq(log(earnings) ~ age+degree+gender, data = data92, tau = 0.75)



stargazer(q2, q4, q6, title = "Wyniki regresji kwantylowych", type = "text")

## 
## Wyniki regresji kwantylowych
## ============================================
##                     Dependent variable:     
##                -----------------------------
##                        log(earnings)        
##                   (1)       (2)       (3)   
## --------------------------------------------
## age            0.020***  0.025***  0.027*** 
##                 (0.003)   (0.002)   (0.003) 
##                                             
## degreebachelor 0.398***  0.396***  0.369*** 
##                 (0.016)   (0.014)   (0.015) 
##                                             
## genderfemale   -0.181*** -0.195*** -0.212***
##                 (0.016)   (0.014)   (0.015) 
##                                             
## Constant       1.540***  1.714***  1.949*** 
##                 (0.085)   (0.075)   (0.079) 
##                                             
## --------------------------------------------
## Observations     5,911     5,911     5,911  
## ============================================
## Note:            *p<0.1; **p<0.05; ***p<0.01

Wyniki regresji kwantylowych kształtują się w sposób przedstawiony w tabeli powyżej

my_qr <- rq(earnings ~ age + gender, data = data92, tau = seq(0.25, 0.75, 0.25))

# Przetwarzanie wyników współczynników
intercept_slope <- coef(my_qr) %>%
  t() %>%
  data.frame()

# Ustawianie poprawnych nazw kolumn
colnames(intercept_slope) <- c("intercept", "slope_age", "slope_gender_female")

# Dodanie kolumny z kwantylami
intercept_slope$quantile <- rownames(intercept_slope)

# Wykres regresji kwantylowej
ggplot() +
  geom_point(data = data92, aes(x = age, y = earnings, color = gender), alpha = 0.5) +
  geom_abline(data = intercept_slope, aes(intercept = intercept, slope = slope_age, color = quantile)) +
  theme_minimal() +
  labs(x = "Wiek", y = "Zarobki", title = "Regresje kwantylowe z tau = 0.25, 0.50 oraz 0.75")

Wykres przedstawiający kształtowanie się zarobków w zależności od wieku oraz płci
Przedstawione są tu również regresje kwantylowe z tau 0.25, 0.50 i 0.75


```r
kmnk <- lm(log(earnings) ~ age + factor(degree) + factor(gender), data = data92)
summary(kmnk)

## 
## Call:
## lm(formula = log(earnings) ~ age + factor(degree) + factor(gender), 
##     data = data92)
## 
## 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 ***
## factor(degree)bachelor  0.38277    0.01173   32.64   <2e-16 ***
## factor(gender)female   -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

kwantyle <- c(0.25, 0.50, 0.75)
reg_kwantylowa <- rq(log(earnings) ~ age+degree+gender,tau = kwantyle,data = data92)
summary(reg_kwantylowa, se = "boot")

## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = data92)
## 
## tau: [1] 0.25
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.53975   0.09329   16.50566   0.00000
## age              0.02026   0.00320    6.32616   0.00000
## degreebachelor   0.39850   0.01732   23.00156   0.00000
## genderfemale    -0.18133   0.01695  -10.69557   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = data92)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.71413   0.08190   20.92888   0.00000
## age              0.02479   0.00273    9.08393   0.00000
## degreebachelor   0.39563   0.01367   28.93719   0.00000
## genderfemale    -0.19526   0.01379  -14.16278   0.00000
## 
## Call: rq(formula = log(earnings) ~ age + degree + gender, tau = kwantyle, 
##     data = data92)
## 
## tau: [1] 0.75
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.94914   0.07226   26.97351   0.00000
## age              0.02666   0.00238   11.21714   0.00000
## degreebachelor   0.36898   0.01442   25.59046   0.00000
## genderfemale    -0.21198   0.01400  -15.14580   0.00000

Zarówno dla modelu wykonanego KMNK jak i wszystkich modeli regresji kwantylowej, wszystkie zmienne objaśniające istotnie wpływają na kształtowanie się zarobków.

#Badanie statystycznej różnicy między 25 i 50 kwantylem warunkowym
kwantyle <- c(0.25, 0.50)
reg_kwantylowa <- rq(log(earnings) ~ age+degree+gender,tau = kwantyle,data = data92)
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.25 0.5  }
## 
##   Df Resid Df F value Pr(>F)
## 1  3    11819  1.5073 0.2104

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

## Quantile Regression Analysis of Deviance Table
## 
## Model: log(earnings) ~ age + degree + gender
## Tests of Equality of Distinct Slopes: tau in {  0.25 0.5  }
## 
##                Df Resid Df F value  Pr(>F)  
## age             1    11821  3.3290 0.06809 .
## degreebachelor  1    11821  0.0416 0.83830  
## genderfemale    1    11821  0.9879 0.32028  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Test wskazuje że 25 i 50 kwantyl warunkowy nie są statystycznie różne

kwantyle <- c(0.25, 0.50, 0.75)

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

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

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

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

## Quantile Regression Analysis of Deviance Table
## 
## Model: log(earnings) ~ age + factor(degree) + factor(gender)
## Tests of Equality of Distinct Slopes: tau in {  0.25 0.5 0.75  }
## 
##                        Df Resid Df F value Pr(>F)
## age                     2    17731  2.2269 0.1079
## factor(degree)bachelor  2    17731  2.1488 0.1167
## factor(gender)female    2    17731  1.4965 0.2239

Test na kwartylach pokazuje że dalej nie występują istotne statystycznie różnice między modelami dla poszczególnych kwartyli

Na koniec sprawdzamy miary dobroci dopasowania

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

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

goodfit(reszty1, reszty2, 0.1)

## [1] 0.07776773

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

summary(model_lm)$r.squared

## [1] 0.1828058

Wartość miary dobroci dopasowania regresji kwantylowej wyniosła 0.06987542 a więc możemy stwierdzić że odchylenia bezwzględne w modelach w pełni sparametryzowanych były w mniejsze od odchyleń bezwzględnych w modelu kwantylowym bezwarunkowym. Wartość dopasowania r2 miała dość małą wartość 0.1871756, więc tylko około 18% zmienności kształtowania się wynagrodzeń zostało objaśnione przez ten model. Zarówno KMNK jak i regresja kwantylowa miały do siebie dość zbliżone wyniki jeśli chodzi o istotność zmiennych oraz ich oddziaływanie na zarobki.