Regresja od zera — tutorial z prawdziwymi danymi

Cel tutorialu

Zbudujemy model wyjaśniający zróżnicowanie płac w USA w 1985 roku.

Pytanie badawcze: Jak wykształcenie, doświadczenie i płeć wpływają na wynagrodzenie?

Dane: Current Population Survey 1985 (CPS1985) — 534 obserwacje, 11 zmiennych.

Krok 1: Załaduj i zbadaj dane

# R
library(AER)
data(CPS1985)
str(CPS1985)

# Python
from statsmodels.datasets import cpunish
# CPS1985 nie jest w statsmodels, ale można z AER przez rpy2 lub pliku CSV
# Alternatywa: użyjemy wbudowanego wage dataset z statsmodels
df = sm.datasets.wage.load_pandas().data
df.head()
df.describe()

Zmienne w CPS1985

Krok 2: Eksploracja danych (EDA)

# Statystyki opisowe
summary(CPS1985)

# Rozkład płac
ggplot(CPS1985, aes(x = wage)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  geom_vline(xintercept = median(CPS1985$wage), 
             color = "red", linestyle = "dashed", linewidth = 1) +
  labs(title = "Rozkład płac (CPS 1985)",
       subtitle = "Pionowa linia = mediana",
       x = "Płaca (USD/godz.)", y = "Liczba obserwacji") +
  theme_minimal()

Co widzimy: Rozkład prawoskośny (skewed right) — kilka osób zarabia bardzo dużo. To typowe dla danych płacowych.

# Sprawdź log-płace
ggplot(CPS1985, aes(x = log(wage))) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  labs(title = "Rozkład log-płac", x = "log(płaca)") +
  theme_minimal()

Log-płace mają rozkład bliższy normalnemu → modelujemy log(wage).

# Korelacje
dane_num <- CPS1985[, c("wage", "education", "experience", "age")]
cor(dane_num)

# Scatter matrix
pairs(dane_num, main = "Macierz scatter plotów")

# Płace według płci
ggplot(CPS1985, aes(x = gender, y = wage, fill = gender)) +
  geom_boxplot() +
  scale_fill_manual(values = c("#0f2044", "#2563eb")) +
  labs(title = "Płace według płci",
       x = "Płeć", y = "Płaca (USD/godz.)") +
  theme_minimal()

Obserwacja: Mężczyźni zarabiają więcej — różnica wynosi ok. 2 USD/godz.

Krok 3: Budowa modeli

Model 1: Prosta regresja

# M1: tylko edukacja
m1 <- lm(log(wage) ~ education, data = CPS1985)
summary(m1)

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.55791    0.10847   5.144 3.72e-07 ***
education    0.09296    0.00813  11.433  < 2e-16 ***

Interpretacja: Każdy dodatkowy rok edukacji wiąże się z ok. 9.3% wyższą płacą.

Model 2: Regresja wieloraka

# M2: edukacja + doświadczenie
m2 <- lm(log(wage) ~ education + experience + I(experience^2),
         data = CPS1985)
summary(m2)

Dodajemy $\text{experience}^2$ bo efekt doświadczenia może być malejący (zbyt duże = przestarzałe umiejętności).

Model 3: Pełny model

# M3: pełny model z kontrolami
m3 <- lm(log(wage) ~ education + experience + I(experience^2) +
           gender + ethnicity + union + married,
         data = CPS1985)
summary(m3)

Porównanie modeli

# Tabela porównawcza
stargazer(m1, m2, m3,
          type = "text",
          title = "Determinanty płac (CPS 1985)",
          dep.var.labels = "log(płaca)",
          covariate.labels = c("Edukacja", "Doświadczenie",
                                "Dośw. kwadrat", "Płeć: kobieta",
                                "Rasa: inna", "Rasa: hispańska",
                                "Związek zawodowy", "Żonaty/zamężna"),
          omit.stat = c("f", "ser"),
          no.space = TRUE)

Krok 4: Diagnostyka modelu

4.1. Wykresy diagnostyczne

par(mfrow = c(2, 2))
plot(m3)

Sprawdzamy 4 wykresy:

1. Reszty vs. dopasowane — szukamy wzorca (heteroskedastyczność?)
2. Normal Q-Q — czy reszty są normalne?
3. Scale-Location — stała wariancja?
4. Cook’s Distance — obserwacje wpływowe?

4.2. Test heteroskedastyczności

library(lmtest)
bptest(m3)    # Test Breusha-Pagana

Jeśli $p < 0.05$ → używamy odpornych błędów standardowych:

library(sandwich)
coeftest(m3, vcov = vcovHC(m3, type = "HC3"))

4.3. Normalność reszt

shapiro.test(residuals(m3))   # Test Shapiro-Wilka (dla n <= 5000)

Przy $n = 534$ normalność jest ważna dla dokładności testów $t$.

4.4. Multikolinearność

library(car)
vif(m3)   # VIF > 10 = poważny problem

Krok 5: Interpretacja wyników

# Podsumowanie modelu 3 z odpornymi SE
coeftest(m3, vcov = vcovHC(m3, type = "HC3"))

                    Estimate Std. Error t value  Pr(>|t|)    
(Intercept)        0.514762   0.140734   3.659 0.0002737 ***
education          0.090561   0.009207   9.836  < 2e-16  ***
experience         0.020613   0.004685   4.401 1.274e-05 ***
I(experience^2)   -0.000283   0.000099  -2.856 0.0044267 ** 
genderfemale      -0.234082   0.038265  -6.117 1.644e-09 ***

Interpretacja każdego współczynnika

Edukacja (0.091): Każdy dodatkowy rok edukacji → 9.1% wyższe płace, przy kontroli pozostałych zmiennych.

Doświadczenie (0.021, -0.000283): Efekt malejący:

• 10 lat dośw. → $+20.1\% - 0.28\% \approx +19.7\%$
• 30 lat dośw. → $+61.8\% - 2.5\% \approx +56.2\%$

Punkt przegięcia: $\text{dośw}^* = \frac{0.0206}{2 \times 0.000283} \approx 36$ lat

Płeć kobieta (-0.234): Kobiety zarabiają o $e^{-0.234} - 1 \approx -20.9\%$ mniej od mężczyzn, ceteris paribus.

To gender pay gap po kontroli wykształcenia i doświadczenia.

Krok 6: Prognozowanie

# Przewidź płacę dla nowej obserwacji
nowa_obs <- data.frame(
  education = 16,          # licencjat
  experience = 10,         # 10 lat pracy
  gender = "male",
  ethnicity = "cauc",
  union = "no",
  married = "yes"
)

predict_log <- predict(m3, newdata = nowa_obs, interval = "prediction")
exp(predict_log)   # przelicz z log na USD/godz.

     fit    lwr    upr
1  10.84   5.12  22.95

Prognozowana płaca: 10.84 USD/godz. (95% PI: 5.12 – 22.95)

Szeroki przedział — to normalny dla danych płacowych ze względu na dużą wariancję nieobserwowalną.

Wnioski

1. Edukacja ma silny, istotny statystycznie efekt (+9.1%/rok)
2. Doświadczenie ma efekt malejący — przełom ok. 36 lat
3. Luka płacowa wynosi ok. 21% na niekorzyść kobiet, nawet po kontroli
4. Model wyjaśnia ok. 30% wariancji log-płac (R² adj. ≈ 0.30)

Uwaga: Wyniki są opisowe, nie przyczynowe — mogą istnieć nieobserwowalne zmienne (np. zdolności) korelujące z wykształceniem.

Następny tutorial: Szeregi czasowe — prognozowanie PKB