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.

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.

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 wczytujemy dane:

data("CPSSW9298")
#?CPSSW9298 

Zestaw danych CPSSW9298 jest używany do analizy rynku pracy. Obejmuje następujące zmienne:

-year: Rok, w którym dane zostały zebrane -earnings: Wynagrodzenie badanej osoby (wyrażone w tysiącach) -degree: Poziom wykształcenia badanej osoby, np. “highschool” (szkoła średnia), “bachelor” (tytuł licencjata) -gender: Płeć osoby, np. “male” (mężczyzna), “female” (kobieta) -age: Wiek badanej osoby

Wizualizacja danych

Teraz zostanie przeprowadzona analiza kilku wykresów, aby zbadać i lepiej zrozumieć zależności pomiędzy różnymi zmiennymi.

ggplot(CPSSW9298, aes(x = earnings)) +
  geom_histogram(bins = 30, fill = "mediumpurple", color = "black") +
  labs(
    title = "Rozkład wynagrodzeń",
    x = "Wynagrodzenie (w tysiącach)",
    y = "Liczba obserwacji"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, color = "black"),
    axis.title = element_text(size = 12, color = "black"),
    axis.text = element_text(size = 10)
  )

Rozkład wynagrodzeń jest asymetryczny, z wyraźnym prawostronnym ogonem, co oznacza, że większa liczba osób zarabia niższe wynagrodzenia, podczas gdy wyższe wynagrodzenia są rzadsze. Kształt sugeruje, że wynagrodzenia mogą być rozkładem skośnym (skośność prawostronna).Najwięcej osób zarabia około 10 tysięcy (zauważalny najwyższy słupek histogramu w tej okolicy).

ggplot(CPSSW9298, aes(x = degree, y = earnings, fill = degree)) +
  geom_boxplot() +
  scale_fill_manual(values = c("mediumpurple", "gold")) + 
  labs(
    title = "Wynagrodzenie względem wykształcenia",
    x = "Poziom wykształcenia",
    y = "Wynagrodzenie (w tysiącach)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, color = "black"),  
    axis.title = element_text(size = 12, color = "black"),  
    axis.text = element_text(size = 10),
    legend.position = "none"
  )

Osoby z wykształceniem “bachelor” (licencjackim) mają wyższą medianę wynagrodzeń w porównaniu do osób z wykształceniem “highschool” (średnim).Wynagrodzenia dla obu grup są zróżnicowane, ale rozstęp międzykwartylowy (odległość między dolnym i górnym kwartylem) jest większy dla osób z wykształceniem licencjackim, co wskazuje na większe zróżnicowanie wynagrodzeń w tej grupie.W obu grupach widać obecność wartości odstających (outliers), czyli osób zarabiających znacznie więcej niż większość. Liczba wartości odstających jest wyższa w grupie z wykształceniem “bachelor”.

ggplot(CPSSW9298, aes(x = year, y = earnings)) +
  stat_summary(fun = "mean", geom = "line", aes(group = 1), color = "darkorchid", linewidth = 1) +
  stat_summary(fun = "mean", geom = "point", color = "purple", size = 2) +
  labs(
    title = "Średnie wynagrodzenie w czasie",
    x = "Rok",
    y = "Średnie wynagrodzenie (w tysiącach)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, color = "black"),
    axis.title = element_text(size = 12, color = "black"),
    axis.text = element_text(size = 10)
  )

Średnie wynagrodzenie systematycznie rosło w czasie w analizowanym okresie (1992–1998). Taki wzrost może wynikać z ogólnej poprawy gospodarki, inflacji, lub zmian w strukturze rynku pracy.

ggplot(CPSSW9298, aes(x = gender, y = earnings, fill = gender)) +
  geom_boxplot() +
  scale_fill_manual(values = c("mediumpurple", "gold")) +
  labs(
    title = "Wynagrodzenie względem płci",
    x = "Płeć",
    y = "Wynagrodzenie (w tysiącach)"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, color = "black"),
    axis.title = element_text(size = 12, color = "black"),
    axis.text = element_text(size = 10),
    legend.position = "none"
  )

Mężczyźni (“male”) mają wyższą medianę wynagrodzeń w porównaniu do kobiet (“female”).Wynagrodzenia zarówno mężczyzn, jak i kobiet charakteryzują się obecnością wartości odstających, ale mężczyźni mają większy rozstęp maksymalnych wynagrodzeń.

ggplot(CPSSW9298, aes(x = earnings, fill = gender)) +
  geom_histogram(alpha = 0.7, bins = 30, position = "identity") +
  scale_fill_manual(values = c("mediumpurple", "gold")) +
  labs(
    title = "Rozkład wynagrodzeń według płci",
    x = "Wynagrodzenie (w tysiącach)",
    y = "Liczba obserwacji"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12),
    legend.title = element_blank()
  )

Wynagrodzenia dla obu płci mają prawostronną skośność – większość osób zarabia w dolnym zakresie wynagrodzeń (poniżej 20 tysięcy), a liczba osób zarabiających więcej systematycznie maleje. Mężczyźni osiągają szerszy zakres wynagrodzeń, w tym najwyższe zarobki, co potwierdzają dłuższe „ogony” fioletowego rozkładu.

ggplot(CPSSW9298, aes(x = age, y = earnings)) +
  geom_point(alpha = 0.7, color = "mediumpurple") +
  geom_smooth(aes(x = age, y = earnings), method = "lm", formula = y ~ x, se = FALSE, linetype = "dashed", color = "gold") +
  labs(
    title = "Wynagrodzenie względem wieku",
    x = "wiek",
    y = "wynagrodzenie (w tysiącach)"
  ) +
  facet_grid(gender ~ degree) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16, face = "bold"),
    axis.title = element_text(size = 14),
    axis.text = element_text(size = 12)
  )

Powyższy wykres potwierdza dotychczasowe obserwacje.Osoby z wykształceniem “bachelor” (licencjackim) osiągają wyższe wynagrodzenia w porównaniu do osób z wykształceniem “highschool” (średnim). Jest to widoczne zarówno w przypadku mężczyzn, jak i kobiet.Różnica między wynagrodzeniami mężczyzn i kobiet jest widoczna zarówno w grupie osób z wykształceniem średnim, jak i licencjackim.

Model liniowy

model_lm <- lm(CPSSW9298$earnings ~ CPSSW9298$degree + CPSSW9298$gender + CPSSW9298$age)
summary(model_lm)  

## 
## Call:
## lm(formula = CPSSW9298$earnings ~ CPSSW9298$degree + CPSSW9298$gender + 
##     CPSSW9298$age)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.023  -3.735  -0.819   2.755  33.526 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               1.70477    0.52435   3.251  0.00115 ** 
## CPSSW9298$degreebachelor  4.91123    0.09938  49.418  < 2e-16 ***
## CPSSW9298$genderfemale   -2.24217    0.09902 -22.643  < 2e-16 ***
## CPSSW9298$age             0.33245    0.01740  19.106  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.665 on 13497 degrees of freedom
## Multiple R-squared:  0.1877, Adjusted R-squared:  0.1876 
## F-statistic:  1040 on 3 and 13497 DF,  p-value: < 2.2e-16

Interpretacje:

• Osoby z wykształceniem licencjackim zarabiają średnio o 4.91 tys. więcej niż osoby z wykształceniem średnim (przy stałej płci i wieku).
• Kobiety zarabiają średnio o 2.24 tys. mniej niż mężczyźni (przy stałym poziomie wykształcenia i wieku).
• Każdy dodatkowy rok życia zwiększa zarobki średnio o 0.33 tys. (przy stałej płci i wykształceniu).

Model wyjaśnia około 18.77% zmienności w zarobkach. Oznacza to, że model nie uwzględnia wielu innych czynników, które wpływają na zarobki (np. branża, lokalizacja, doświadczenie zawodowe).

Regresja kwantylowa

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

## 
## Call: rq(formula = earnings ~ age + degree + gender, tau = kwantyle, 
##     data = data)
## 
## tau: [1] 0.25
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      2.67308   0.38777    6.89339   0.00000
## age              0.17308   0.01365   12.67920   0.00000
## degreebachelor   3.59188   0.06431   55.84980   0.00000
## genderfemale    -1.32265   0.07062  -18.72785   0.00000
## 
## Call: rq(formula = earnings ~ age + degree + gender, tau = kwantyle, 
##     data = data)
## 
## tau: [1] 0.5
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.45427   0.49704    2.92588   0.00344
## age              0.31389   0.01729   18.15274   0.00000
## degreebachelor   4.70043   0.11085   42.40338   0.00000
## genderfemale    -2.12970   0.09624  -22.12868   0.00000
## 
## Call: rq(formula = earnings ~ age + degree + gender, tau = kwantyle, 
##     data = data)
## 
## tau: [1] 0.75
## 
## Coefficients:
##                Value     Std. Error t value   Pr(>|t|) 
## (Intercept)      1.53846   0.92264    1.66747   0.09544
## age              0.43269   0.02998   14.43478   0.00000
## degreebachelor   5.75962   0.15277   37.70122   0.00000
## genderfemale    -2.97009   0.12878  -23.06241   0.00000

Analiza wpływu wieku, wykształcenia i płci na zarobki w różnych kwantylach dostarcza interesujących wniosków. Wyniki wskazują, że starszy wiek ma większy wpływ na zarobki osób znajdujących się w wyższych kwantylach wynagrodzeń. W dolnym 25% zarobków każdy rok życia zwiększa wynagrodzenie o 0.17 tys., podczas gdy w grupie najlepiej zarabiających (75%) wzrost wynosi już 0.43 tys. Oznacza to, że starszy wiek bardziej korzystnie wpływa na osoby osiągające wyższe dochody.

Wykształcenie licencjackie również odgrywa istotną rolę, szczególnie w medianie i górnych kwantylach. W grupie osób z najniższymi zarobkami wykształcenie licencjackie zwiększa wynagrodzenie średnio o 3.59 tys., ale w medianie i wśród najlepiej zarabiających wpływ ten rośnie do 4.70 tys. Świadczy to o tym, że wyższe wykształcenie ma szczególną wartość wśród osób z wyższymi dochodami.

Różnice w wynagrodzeniach między płciami są widoczne we wszystkich kwantylach, ale najbardziej niepokojące są w górnym 25% zarobków. Kobiety zarabiają mniej niż mężczyźni, a różnice te rosną w zależności od poziomu dochodów: od 1.32 tys. w najniższych zarobkach do 2.97 tys. wśród najlepiej zarabiających.

plot(rq(earnings ~ age + degree + gender, tau = c(0.25, 0.5, 0.75), data = data))

kwantyle <- c(0.25, 0.50, 0.75)
model_q25 <- rq(earnings ~ age + degree + gender, tau = 0.25, data = data)
model_q50 <- rq(earnings ~ age + degree + gender, tau = 0.50, data = data)
model_q75 <- rq(earnings ~ age + degree + gender, tau = 0.75, data = data)
coef_q25 <- coef(model_q25)
coef_q50 <- coef(model_q50)
coef_q75 <- coef(model_q75)

coef_table <- data.frame(
  Kwantyl = c(0.25, 0.50, 0.75),
  Intercept = c(coef_q25[1], coef_q50[1], coef_q75[1]),
  Age = c(coef_q25["age"], coef_q50["age"], coef_q75["age"]),
  Degree = c(coef_q25["degreebachelor"], coef_q50["degreebachelor"], coef_q75["degreebachelor"]),
  Gender = c(coef_q25["genderfemale"], coef_q50["genderfemale"], coef_q75["genderfemale"])
)

print(coef_table)

##   Kwantyl Intercept       Age   Degree    Gender
## 1    0.25  2.673078 0.1730769 3.591880 -1.322649
## 2    0.50  1.454274 0.3138889 4.700427 -2.129701
## 3    0.75  1.538464 0.4326922 5.759616 -2.970085

Tabela przedstawia wyniki regresji kwantylowej dla trzech wybranych kwantyli: 0.25 (dolne 25% zarobków), 0.50 (mediana zarobków), i 0.75 (górne 25% zarobków). Jest to analiza, która pozwala zobaczyć, jak zmienne niezależne wpływają na zarobki w różnych częściach rozkładu wynagrodzeń, a nie tylko na średnią (jak w klasycznej regresji liniowej).

Weryfikacja statystyczna istotności różnic między pierwszym, drugim i trzecim kwartylem

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

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

Hipoteza zerowa: Współczynniki dla kwantyli są takie same. Hipoteza alternatywna: Współczynniki różnią się między kwantylami.

Ponieważ wartość p jest bardzo niska, odrzucamy hipotezę zerową. Oznacza to, że wpływ zmiennych (wiek, wykształcenie, płeć) różni się w różnych częściach rozkładu zarobków.

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

## Quantile Regression Analysis of Deviance Table
## 
## Model: earnings ~ age + degree + 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  72.465 < 2.2e-16 ***
## degreebachelor  2    40501 117.817 < 2.2e-16 ***
## genderfemale    2    40501  88.771 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Wszystkie zmienne (age, degreebachelor, genderfemale) mają istotnie różny wpływ na zarobki w zależności od kwantyla.

Dobroć dopasowania

# Obliczenie reszt dla modeli
reszty_lm <- resid(model_lm)  
model2 <- rq(earnings ~ 1, tau = 0.25, data = data)  
reszty_2 <- resid(model2) 
model3 <- rq(earnings ~ 1, tau = 0.75,data=data)
reszty3 <- resid(model3)

par(mfrow = c(1, 3))  # Podziel przestrzeń na 3 kolumny

# Histogram dla modelu liniowego
hist(reszty_lm,
     main = "Histogram reszt (model liniowy)",
     xlab = "Reszty",
     col = "black", border = "grey")

# Histogram dla modelu kwantylowego tau = 0.25
hist(reszty_2,
     main = "Histogram reszt (tau = 0.25)",
     xlab = "Reszty",
     col = "gold", border = "black")

# Histogram dla modelu kwantylowego tau = 0.75
hist(reszty3,
     main = "Histogram reszt (tau = 0.75)",
     xlab = "Reszty",
     col = "purple", border = "black")

coef_long <- reshape2::melt(coef_table, id.vars = "Kwantyl", 
                            variable.name = "Zmienna", value.name = "Współczynnik")
ggplot(coef_long, aes(x = Kwantyl, y = Współczynnik, color = Zmienna)) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  labs(
    title = "Współczynniki regresji kwantylowej",
    x = "Kwantyl",
    y = "Współczynnik",
    color = "Zmienna"
  ) +
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5, size = 16),
    axis.title = element_text(size = 12),
    legend.title = element_text(size = 12)
  )

## 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.

Wpływ zmiennych jest zróżnicowany w różnych częściach rozkładu zarobków, co wskazuje na to, że klasyczna regresja liniowa (która zakłada stały wpływ zmiennych na całość rozkładu) mogłaby nie uchwycić tych różnic.

-Wpływ wieku rośnie wraz ze wzrostem kwantyla. Starszy wiek ma największy wpływ na zarobki w grupie najlepiej zarabiających, co jest zgodne z oczekiwaniem, że doświadczenie zawodowe bardziej się opłaca w tej grupie. - Wykształcenie ma rosnący wpływ w miarę wzrostu kwantyla. Sugeruje to, że wartość wyższego wykształcenia jest bardziej doceniana w grupie osób osiągających najwyższe dochody. - Wartości ujemne współczynnika Gender wskazują na niższe zarobki kobiet w porównaniu do mężczyzn.