Zastosowanie modelu selekcji Heckmana do modelowania aktywności zawodowej kobiet
1 Wstęp
Celem niniejszej pracy jest empiryczna analiza wpływu selekcji próby na oszacowania modelu płacowego wśród kobiet, z wykorzystaniem klasycznego modelu selekcji Heckmana. W badaniu wykorzystano dane z przekroju panelowego PSID (zbiór Mroz87), zawierające informacje o zamężnych kobietach w USA w 1975 roku. Analiza obejmuje porównanie wyników estymacji klasycznej regresji OLS oraz modelu Heckmana, który uwzględnia możliwość systematycznego różnicowania się populacji aktywnej i nieaktywnej zawodowo. Główna hipoteza badawcza zakłada, że nieuwzględnienie mechanizmu selekcji prowadzi do istotnie odmiennych wniosków niż model z korektą selekcyjną. Estymacja została przeprowadzona metodą dwustopniową, zgodnie z oryginalną propozycją Heckmana. Co istotne, procedura ta została w całości zaimplementowana samodzielnie w języku R – bez użycia gotowych funkcji pakietów do modeli selekcji, co pozwoliło na pełną kontrolę nad obliczeniami i przeprowadzenie niezależnej diagnostyki każdego z etapów modelu. Praca opiera się na założeniu, że część czynników wpływa jedynie na decyzję o uczestnictwie w rynku pracy, nie oddziałując bezpośrednio na poziom wynagrodzenia. Weryfikacja hipotezy pozwala ocenić zasadność stosowania bardziej złożonych narzędzi ekonometrycznych w analizach ekonomii pracy.
2 Przegląd literatury
Zagadnienie uczestnictwa w rynku pracy było przedmiotem intensywnych analiz ekonomicznych od lat 60. XX wieku, jednak przełomowe podejścia pojawiły się w kolejnych dekadach, wraz z rosnącym udziałem kobiet w zatrudnieniu. Wczesne modele mikroekonometryczne, takie jak te opracowane przez Mincera (1974), koncentrowały się na zależności między edukacją, doświadczeniem zawodowym a płacą. Mincer zaproponował tzw. funkcję zarobków, która jednak zakładała dostępność danych płacowych dla wszystkich jednostek. Takie podejście ignorowało problem selekcji próby, istotny w sytuacji, gdy płace obserwuje się tylko dla osób pracujących.
Gronau (1974) jako jeden z pierwszych badaczy wskazał na potrzebę jednoczesnego modelowania decyzji o uczestnictwie oraz poziomu zarobków, zauważając, że kobiety aktywne zawodowo mogą różnić się od tych nieaktywnych w sposób systematyczny (np. pod względem preferencji, umiejętności czy sytuacji rodzinnej). To podejście zostało rozwinięte przez Heckmana (1979), który w klasycznym artykule zaproponował dwustopniowy model selekcji próby, znany dziś jako model Heckmana. Jego podejście pozwalało uwzględnić zależność pomiędzy składnikami zakłócającymi decyzję o zatrudnieniu i poziomem płacy. Kluczową innowacją było wprowadzenie tzw. „korekty selekcji” (lambda Heckmana), która umożliwia uzyskanie nieskorygowanych estymatorów parametrów w równaniu płacowym.
W kolejnych dekadach model Heckmana był szeroko stosowany w analizach rynku pracy. Killingsworth i Heckman (1986) przedstawili przegląd teoretyczny i empiryczny dotyczący uczestnictwa kobiet w rynku pracy, zwracając uwagę na rolę edukacji, liczby dzieci, statusu małżeńskiego oraz warunków instytucjonalnych. Badania te potwierdzały, że ignorowanie selekcji próby prowadzi do zaniżenia lub zawyżenia oszacowań współczynników w modelach płacowych.
Ekonometrzy tacy jak Wooldridge (2010) czy Cameron i Trivedi (2005) w swoich podręcznikach metod ilościowych podkreślają, że model Heckmana ma zastosowanie nie tylko w analizie płac, ale w każdej sytuacji, w której obserwujemy zmienną zależną tylko dla wybranej części populacji (tzw. cenzurowanie lub selekcja próby). Obaj autorzy omawiają też ograniczenia tego podejścia, w tym silne założenie o normalności rozkładu błędów i potrzebę dobrego instrumentu w równaniu selekcji – zmiennej, która wpływa na decyzję o pracy, ale nie bezpośrednio na płacę (np. liczba małych dzieci). Empiryczne badania oparte na danych PSID (Panel Study of Income Dynamics), takich jak w niniejszym opracowaniu, są powszechnie cytowane w literaturze. Na przykład, Blau i Kahn (2007) analizowali różnice w partycypacji zawodowej kobiet między USA a Europą, wykorzystując dane z PSID oraz Europejskiego Panelu Gospodarstw Domowych. Udowodnili, że różnice w politykach podatkowych, urlopach macierzyńskich i dostępności usług opiekuńczych wpływają na poziom aktywności zawodowej kobiet.
Ponadto, najnowsze badania empiryczne podkreślają, że tradycyjne modele regresji – nawet po kontroli zmiennych demograficznych, nie są wystarczające do uchwycenia mechanizmów decyzyjnych związanych z podjęciem pracy (Bertrand, Goldin, i Katz 2010; Goldin 2014). Modele strukturalne, takie jak Heckman, pozostają zatem niezbędnym narzędziem analizy w sytuacjach, gdy dane są ograniczone do aktywnych zawodowo, a zainteresowanie badawcze dotyczy całej populacji.
3 Metodologia
3.1 Model selekcji Heckmana
W badaniach empirycznych, szczególnie dotyczących rynku pracy, często występuje problem selekcji próby (sample selection bias). Typowym przypadkiem jest analiza płac, gdzie dane o wynagrodzeniach dostępne są tylko dla osób aktywnych zawodowo. Jeżeli osoby pracujące różnią się systematycznie od tych, które nie pracują, estymacja modeli przy użyciu klasycznych metod (np. OLS) prowadzi do obciążonych wyników. Aby temu przeciwdziałać, James Heckman zaproponował dwustopniowy model selekcji próby, znany jako model Heckmana (Heckman selection model) (Heckman 1976).
Model Heckmana składa się z dwóch równań:
- Równanie nieobserwowalnego procesu selekcji (selection equation – np. decyzja o uczestnictwie w rynku pracy), estymowane za pomocą modelu probitowego:
\[ \mathbf{z^*} = \mathbf{w}\gamma + \mathbf{u} \]
w którym obserwujemy zmienną zero-jedynkową \(z_i\):
\[ z_i = \begin{cases} 1 & \text{dla } z^*_i > 0 \\ 0 & \text{dla } z^*_i \leq 0 \end{cases} \]
- Równanie wyniku (outcome equation):
\[ \mathbf{y} = \mathbf{x\beta} + \mathbf{\varepsilon} \]
gdzie:
- \(\mathbf{u} \sim N(0,1)\)
- \(\mathbf{\varepsilon} \sim N(0,\sigma)\)
- \(\operatorname{corr}(\mathbf{u}, \mathbf{\varepsilon}) = \rho\)
3.2 Estymacja dwustopniowa (two-step estimation)
W celu przeprowadzenia pogłębionej diagnostyki obu regresji, postanowiliśmy zaimplementować model Heckmana samodzielnie, nie używając wbudowanych funkcji lub dodatkowych bibliotek w pakietach statystycznych. Ręczne przeprowadzenie tej procedury pozwala na zastosowanie popularnych narzędzi/testów diagnostycznych dla obu modeli oddzielnie. Jest to istotne szczególnie w języku R, gdyż dedykowany pakiet do modeli selekcji sampleSelection (Toomet i Henningsen 2008) ma dość ograniczony zbiór funkcji.
Większość implementacji modelu Heckmana w popularnych pakietach statystycznych – zarówno polecenie rheckman w STATA jak i pakiet sampleSelection w R – przyjmuje MNW jako domyślną metodę estymacji. W kontekście samodzielnej implementacji modelu prostszym rozwiązaniem jest estymacja dwustopniowa (two-step estimation), będąca rozwiązaniem oryginalnie zaproponowanym przez Heckmana.
Kolejne etapy procedury prezentuja się następująco:
- Model probit równania selekcji dostarcza oszacowania prawdopodobieństw:
\[ Pr(y_i \mid \mathbf{w_i}) = \Phi(\mathbf{w_i}\gamma) \]
- Z tych oszacowań liczone jest prawdopodobieństwo, że obserwacja nie zostanie włączona do próby, co Heckman nazwał odwrotnością współczynnika Millsa (inversed Mills ratio - IMR)
\[ IMR = \frac{\phi(\mathbf{w_i}\hat{\gamma})}{\Phi(\mathbf{w_i}\hat{\gamma})} \]
Równanie wynikowe rozbudowujemy o dodatkową zmienną \(IMR\) oznaczającą ryzyko niewybrania (nonselection hazard). Otrzymujemy macierz zmiennych objaśniających \([X \quad m]\), a co za tym idzie – dodatkowe oszacowanie \(\beta_{IMR}\).
Obliczenie zgodnego estymatora wariancji błędu w regresji uwzględniającej IMR jest możliwe, ale wykracza poza ramy tej pracy.
4 Dane
4.1 Opis zbioru danych
W niniejszym badaniu wykorzystano znany i szeroko używany w ekonometrii do celów edukacyjnych zbiór danych Mroz87, oryginalnie opracowany na potrzeby artykułu autorstwa Mroza (1987). Zbiór jest dołączony do pakietu sampleSelection, skąd został załadowany do środowiska R.
Źródłem danych jest Panel Study of Income Dynamics (PSID) – jeden z najdłużej prowadzonych paneli gospodarstw domowych na świecie, prowadzony przez University of Michigan. Dane użyte w tym badaniu pochodzą z przekroju za rok 1976, który odnosi się do sytuacji respondentów w roku 1975. Zbiór zawiera 753 obserwacje dotyczące zamężnych kobiet i składa się z 22 zmiennych demograficznych, społecznych i ekonomicznych, które opisuje Tabela 1.
Mroz87
| Zmienna | Opis |
|---|---|
| lfp | Zmienna zero-jedynkowa wskazująca na uczestnictwo żony w rynku pracy |
| hours | Liczba godzin przepracowanych przez żonę w 1975 roku |
| kids5 | Liczba dzieci w wieku 5 lat lub młodszych |
| kids618 | Liczba dzieci w wieku od 6 do 18 lat |
| age | Wiek żony |
| educ | Liczba lat edukacji żony |
| wage | Średnie godzinowe wynagrodzenie żony w dolarach z 1975 roku |
| repwage | Wynagrodzenie żony zadeklarowane podczas wywiadu w 1976 roku |
| hushrs | Liczba godzin przepracowanych przez męża w 1975 roku |
| husage | Wiek męża |
| huseduc | Liczba lat edukacji męża |
| huswage | Wynagrodzenie męża w dolarach z 1975 roku |
| faminc | Dochód rodziny w dolarach z 1975 roku |
| mtr | Krańcowa stawka podatkowa dla żony |
| motheduc | Liczba lat edukacji matki żony |
| fatheduc | Liczba lat edukacji ojca żony |
| unem | Stopa bezrobocia w hrabstwie zamieszkania (w punktach procentowych) |
| city | Zmienna zero-jedynkowa: 1 = mieszka w dużym mieście, 0 = inaczej |
| exper | Faktyczna liczba lat doświadczenia zawodowego żony |
| nwifeinc | Dochód spoza pracy żony |
| wifecoll | Zmienna zero-jedynkowa: czy żona uczęszczała na studia |
| huscoll | Zmienna zero-jedynkowa: czy mąż uczęszczał na studia |
4.2 Wstępna analiza danych
Przed przystąpieniem do modelowania, kluczowa jest wstępna analiza danych (Exploratory Data Analysis – EDA). Jej celem jest zrozumienie struktury zbioru danych, identyfikacja potencjalnych problemów (braki danych, wartości odstające), a także odkrycie zależności między zmiennymi. Tabela 2 przedstawia podstawowe statystyki opisowe dla wszystkich zmiennych w zbiorze danych.
| variable | Obs | Mean | Std. dev. | Min | Max |
|---|---|---|---|---|---|
| age | 753 | 42.54 | 8.07 | 30.00 | 60.00 |
| city | 753 | 0.64 | 0.48 | 0.00 | 1.00 |
| educ | 753 | 12.29 | 2.28 | 5.00 | 17.00 |
| exper | 753 | 10.63 | 8.07 | 0.00 | 45.00 |
| faminc | 753 | 23080.59 | 12190.20 | 1500.00 | 96000.00 |
| fatheduc | 753 | 8.81 | 3.57 | 0.00 | 17.00 |
| hours | 753 | 740.58 | 871.31 | 0.00 | 4950.00 |
| husage | 753 | 45.12 | 8.06 | 30.00 | 60.00 |
| huscoll | 753 | 0.39 | 0.49 | 0.00 | 1.00 |
| huseduc | 753 | 12.49 | 3.02 | 3.00 | 17.00 |
| hushrs | 753 | 2267.27 | 595.57 | 175.00 | 5010.00 |
| huswage | 753 | 7.48 | 4.23 | 0.41 | 40.51 |
| kids5 | 753 | 0.24 | 0.52 | 0.00 | 3.00 |
| kids618 | 753 | 1.35 | 1.32 | 0.00 | 8.00 |
| motheduc | 753 | 9.25 | 3.37 | 0.00 | 17.00 |
| mtr | 753 | 0.68 | 0.08 | 0.44 | 0.94 |
| nwifeinc | 753 | 20.13 | 11.63 | -0.03 | 96.00 |
| repwage | 753 | 1.85 | 2.42 | 0.00 | 9.98 |
| unem | 753 | 8.62 | 3.11 | 3.00 | 14.00 |
| wage | 753 | 2.37 | 3.24 | 0.00 | 25.00 |
| wifecoll | 753 | 0.28 | 0.45 | 0.00 | 1.00 |
Pierwszym etapem budowy modelu Heckmana jest równanie selekcji. Warto wobec tego porównać zmienne w zbiorze w podziale na kobiety pracujące i niepracujące, aby wykryć cechy najbardziej różnicujące obie te grupy.
Rysunek 1 ilustruje rozkład wybranych zmiennych dyskretnych w podziale na status aktywności zawodowej. Na osi Y przedstawiono odsetek kobiet z danej grupy (pracujących lub niepracujących) charakteryzujących się wartością cechy określonej na osi X. Najbardziej wyraźną różnicą między rozkładami charakteryzuje się zmienna kids5. Znacznie większy odsetek kobiet nieaktywnych zawodowo posiada przynajmniej jedno dziecko w wieku do 5 lat. Pozostałe rozkłady wyglądają podobnie.
Rozkłady zmiennych ciągłych i quasi-ciągłych w podziale na status aktywności zawodowej prezentuje Rysunek 2. Kobiety pracujące charakteryzują się wyższymi wartościami zmiennych takich jak: educ, exper, faminc, huseduc.
Rysunek 3 przedstawia macierz korelacji dla wszystkich zmiennych w zbiorze danych. Zmienne są uporządkowane hierarchicznie na podstawie podobieństwa korelacji, co ułatwia dostrzeżenie skupień zmiennych o silnych wzajemnych zależnościach. Macierz korelacji może również pomóc w identyfikacji potencjalnych problemów współliniowości w danych.
Grupy zmiennych o wysokich korelacjach wewnętrznych:
- Zmienne edukacyjne:
educ,huseduc,motheduc,fatheduc– wszystkie zmienne edukacyjne są dodatnio skorelowane (0.28-0.61)- Wysoka korelacja między edukacją żony a męża (0.61) wskazuje na homogamię edukacyjną
- Zmienne dochodowe
faminc,nwifeinc,huswage– silna dodatnia korelacja- interpretacja: zazwyczaj większość dochodu rodziny stanowi dochód męża
hours,wage,repwage– silna dodatnia korelacja- krańcowa stopa opodatkowania
mtrjest ujemnie skorelowane ze zmiennymi dochodowymi, co może odzwierciedlać progresywny charakter systemu podatkowego.
- Zmienne demograficzne:
ageihusage(0.89) – małżonkowie mają podobny wiekageiexper(0.33) – umiarkowana korelacja- interpretacja: starsze kobiety mają większą szansę mieć większe doświadczenie
5 Hipoteza badawcza
Celem niniejszego badania jest empiryczne sprawdzenie, czy klasyczny model regresji liniowej (OLS) prowadzi do różnych wyników w porównaniu z modelem Heckmana, który uwzględnia selekcję próby. W tym celu estymowane będą dwa modele:
- Model regresji OLS dla logarytmu płacy godzinowej,
- Model Heckmana, w którym pierwszy etap obejmuje decyzję o uczestnictwie w rynku pracy, a drugi etap – poziom płacy.
Hipoteza badawcza:
- \(H_0\): Oszacowania modelu OLS nie różnią się istotnie od oszacowań modelu Heckmana.
- \(H_1\): Oszacowania modelu Heckmana różnią się istotnie od oszacowań OLS – co oznacza, że występuje selekcja próby.
Weryfikacja hipotezy opierać się będzie na porównaniu współczynników w obu modelach oraz na ocenie istotności statystycznej odwrotności współczynnika Millsa (IMR). Odrzucenie hipotezy zerowej będzie świadczyć o konieczności stosowania modelu selekcji zamiast klasycznego OLS.
6 Wyniki
6.1 Równanie selekcji
6.1.1 Model
Poniższy model został wybrany w procedurze ręcznej selekcji zmiennych, uwzględniającej kryteria statystyczne (istotność zmiennych), diagnostyczne, zgodność z teorią i interpretowalność.
Kluczową rolę w identyfikacji modeli selekcji pełnią zmienne, które wpływają na prawdopodobieństwo selekcji, ale niekoniecznie na wysokość płacy. Umożliwia to odróżnienie decyzji o uczestnictwie od wyniku płacowego. Choć teoretycznie identyfikacja może opierać się wyłącznie na założeniach dotyczących funkcji postaci modelu, brak zmiennych identyfikujących podważa wiarygodność wyników, ponieważ takie założenia są często arbitralne i trudne do obrony empirycznej (Heckman 1979).
W poniższym modelu celowo zostały zatem uwzględnione zmienne, które wpływają na decyzję o podjęciu pracy, ale nie przekładają się bezpośrednio na wysokość zarobków. Zmienna I(age - exper) reprezentuje wpływ wielkości “luki” w życiorysie zawodowym na prawdopodobieństwo rozpoczęcia pracy.
W modelu probitowym możemy bezpośrednio interpretować znak współczynników. Z wydruku można wywnioskować, że posiadanie jednego lub dwójki dzieci wpływa negatywnie na decyzję o wstąpieniu na rynek pracy w porównaniu do scenariusza bazowego – braku małych dzieci. Większa liczba godzin pracy i wysokość zarobków męża również negatywnie wpływa na aktywność zawodową. Co ciekawe, posiadanie wykształcenia wyższego przez męża ma pozytywny wpływ na decyzję o podjęciu pracy przez żonę.
\[ \begin{aligned} P( \operatorname{lfp} = \operatorname{1} ) &= \Phi[\beta_{0} + \beta_{1}(\operatorname{kids5}) + \beta_{2}(\operatorname{huswage}) + \beta_{3}(\operatorname{hushrs})\ + \\ &\qquad\ \beta_{4}(\operatorname{age\ -\ \exp}) + \beta_{5}(\operatorname{huscoll}_{\operatorname{TRUE}})] \end{aligned} \]
Call:
glm(formula = select.formula, family = binomial(link = "probit"),
data = data.df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 3.086e+00 3.156e-01 9.780 <2e-16 ***
kids5 -9.084e-01 1.071e-01 -8.481 <2e-16 ***
huswage -2.292e-02 1.351e-02 -1.696 0.0899 .
hushrs -1.985e-04 8.908e-05 -2.228 0.0259 *
I(age - exper) -6.742e-02 6.144e-03 -10.974 <2e-16 ***
huscollTRUE 2.728e-01 1.161e-01 2.350 0.0188 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1029.75 on 752 degrees of freedom
Residual deviance: 835.55 on 747 degrees of freedom
AIC: 847.55
Number of Fisher Scoring iterations: 46.1.2 Diagnostyka
W celu oceny poprawności przyjętej specyfikacji modelu probitowego przeprowadzono szereg testów diagnostycznych.
6.1.2.1 Jakość dopasowania
W pierwszej kolejności oceniono jakość dopasowania za pomocą testu Hosmera-Lemenshowa, który porównuje obserwowane i oczekiwane liczby sukcesów w podgrupach na podstawie przewidywanych prawdopodobieństw (Mycielski 2010). Wysokie p-value w tym teście świadczy o braku istotnych rozbieżności i tym samym sugeruje, że model został dobrze dopasowany.
Hosmer and Lemeshow goodness of fit (GOF) test
data: data.df$lfp, select.fitted.prob
X-squared = 3.8942, df = 8, p-value = 0.8665Ogólną jakość dopasowania modelu za pomocą wskaźnika pseudo-\(R^2\) McFaddena. Choć wartości tego wskaźnika są z reguły niższe niż w klasycznych modelach liniowych, wartości powyżej 0.2 wskazują na względnie dobre dopasowanie modelu binarnego. W przypadku omawianego modelu selekcji wartość pseudo-\(R^2\) McFaddena wynosi \(0.189\) i można ją uznać za akceptowalną.
fitting null model for pseudo-r2 llh llhNull G2 McFadden r2ML r2CU
-417.7733291 -514.8732046 194.1997508 0.1885899 0.2273286 0.3050312 Ocenę uzupełniono analizą krzywej ROC Rysunek 4, która obrazuje zdolność modelu do odróżniania obserwacji z różnych klas. Wartość AUC (area under curve) równa \(0.779\), oznacza akceptowalną trafność klasyfikacyjną modelu.
Area under the curve: 0.778
Jako dodatkowy punkt odniesienia w diagnostyce dopasowania modelu zamieszczono tabelę klasyfikacyjną, ze standardowym progiem odcięcia równym \(0.5\).
actual
predicted 0 1
0 199 96
1 126 3326.1.2.2 Forma funkcyjna
Do testowania formy funkcyjnej zastosowano test typu związku (linktest). Najpierw utworzono zmienną \(\widehat{y}_i^* = x_i\widetilde{\beta}\), a następnie przeprowadzono regresję probitową \(y_i\) na stałej, \(\widehat{y}_i^*\) i \(\left(\widehat{y}_i^*\right)^2\). Niestotny współczynnik przy zmiennej I(y.hat^2) nie pozwala nam odrzucić hipotezy zerowej o poprawności funkcyjnej.
Call:
glm(formula = lfp ~ y.hat + I(y.hat^2), family = binomial(link = "probit"),
data = data.df)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.003561 0.063896 0.056 0.956
y.hat 1.002669 0.083324 12.033 <2e-16 ***
I(y.hat^2) -0.009338 0.096771 -0.096 0.923
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 1029.75 on 752 degrees of freedom
Residual deviance: 835.54 on 750 degrees of freedom
AIC: 841.54
Number of Fisher Scoring iterations: 46.2 Równanie wynikowe
6.2.1 Model
Specyfikacja modelu wynikowego opiera się na regresji metodą MNK logarytmu godzinowej stawki płacy (log(wage)) dla kobiet aktywnych zawodowo (\(\text{lfp} = 1\)), z uwzględnieniem poprawki na błąd wynikający z selekcji próby poprzez dołączenie zmiennej IMR (inversed Mills ratio). Model ten został estymowany na próbie ograniczonej do pracujących kobiet, co wymagało korekty selekcyjnej, by zapewnić nieobciążone oszacowania parametrów.
W modelu przyjęto nieliniową postać zmiennych exper i educ, które są typowymi determinantami płac w teorii kapitału ludzkiego. Doświadczenie zawodowe zostało ujęte w postaci logarytmu skorygowanego o jeden (log(exper + 1)), co pozwala uchwycić malejące przyrosty płacy wraz ze wzrostem doświadczenia. Dodanie \(1\) zapobiega problemom z logarytmem zera, w przypadku gdy exper przyjmuje wartość \(0\). Wpływ edukacji natomiast został odwzorowana za pomocą funkcji ortogonalnego wielomianu drugiego stopnia (poly(educ, 2)), co umożliwia modelowi wychwycenie nieliniowych efektów wykształcenia na płace bez współliniowości między składnikami.
Istotnym elementem estymacji jest zastosowanie odpornej macierzy wariancji (HC3), w celu wyeliminowania problemu heteroskedastyczności.
\[ \begin{aligned} \operatorname{log(wage)} &= \beta_{0} + \beta_{1}(\operatorname{\log(\exp\ +\ 1)}) + \beta_{2}(\operatorname{educ}) + \beta_{3}(\operatorname{educ^2}) + \beta_{4}(\operatorname{mtr})\ + \\ &\quad \beta_{5}(\operatorname{motheduc}) + \beta_{6}(\operatorname{wifecoll}_{\operatorname{TRUE}}) + \beta_{7}(\operatorname{IMR}) + \epsilon \end{aligned} \]
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.715130 0.404935 6.7051 6.511e-11 ***
log(exper + 1) 0.171255 0.059932 2.8575 0.004482 **
poly(educ, 2)1 5.872026 1.170443 5.0169 7.765e-07 ***
poly(educ, 2)2 1.543472 0.666920 2.3143 0.021131 *
mtr -2.463867 0.423175 -5.8223 1.154e-08 ***
motheduc -0.018420 0.010683 -1.7242 0.085404 .
wifecollTRUE -0.250368 0.126681 -1.9764 0.048766 *
IMR -0.066835 0.152577 -0.4380 0.661582
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wyniki modelu regresji są zgodne z teorią ekonomiczną – większe doświadczenie i wykształcenie sprzyjają wyższym zarobkom, natomiast wyższy poziom opodatkowania działa negatywnie.
Najważniejszym wnioskiem metodologicznym jest brak statystycznej istotności współczynnika dla odwróconego współczynnika Millsa (IMR), co sugeruje, że w tym modelu nie ma istotnego obciążenia wynikającego z selekcji próby. Oznacza to, że po uwzględnieniu zmiennych z równania selekcji, czynniki nienamierzalne wpływające na decyzję o podjęciu pracy nie są już istotnie skorelowane z czynnikami nienamierzalnymi wpływającymi na wysokość płacy.
6.2.2 Diagnostyka
6.2.2.1 Normalność reszt
W pierwszym kroku zbadano rozkład reszt, aby ocenić założenie normalności. Wykonano histogram reszt oraz wykres Q-Q (kwantyl-kwantyl, Rysunek 5), który umożliwia wizualną ocenę zgodności rozkładu reszt z rozkładem normalnym. Dodatkowo zastosowano test Shapiro-Wilka, który formalnie weryfikuje hipotezę normalności. P-value bliskie zera każe odrzucić hipotezę zerową o normalności rozkładu. Próba jest jednak na tyle duża, że niespełnienie warunku normalności nie jest problematyczne w praktyce (ze wzlędu na Centralne Twierdzenie Graniczne).
Shapiro-Wilk normality test
data: outcome.res
W = 0.93548, p-value = 1.176e-12
6.2.2.2 Założenie homoskedastyczności
Następnie zbadano założenie homoskedastyczności. W tym celu wygenerowano wykres reszt względem wartości dopasowanych Rysunek 6. Nie wykazuje on wyraźnych systematycznych wzorców ani struktury. Punkty są rozproszone względnie losowo wokół osi poziomej (reszta ≈ 0). Wizualna inspekcja nie dostarcza silnych dowodów na heteroskedastyczność w tym modelu. Natomiast niskie \(p-value\) testu Breuscha-Pagana pozwala odrzucić hipotezę zerową o homoskedastycznej wariancji błędów. Wskazuje to na konieczność zastosowania standardowych błędów odpornych na heteroskedastyczność, co zostało zasygnalizowane w sekcji 6.2.1
studentized Breusch-Pagan test
data: outcome.model
BP = 24.945, df = 7, p-value = 0.0007761
6.2.2.3 Poprawność formy funkcyjnej
W celu wykrycia ogólnych problemów ze specyfikacją modelu, w tym pominiętych zmiennych istotnych (omitted variable bias) oraz niewłaściwej formy funkcyjnej, przeprowadzono test RESET. Wysoka wartość p-value wskazuje na brak podstaw do odrzucenia hipotezy zerowej o poprawnej formie funkcyjnej modelu.
RESET test
data: outcome.model
RESET = 0.43295, df1 = 2, df2 = 418, p-value = 0.64896.2.2.4 Współliniowość
Na końcu oszacowano wartości współczynników VIF (Variance Inflation Factor) dla modelu wynikowego. Na wydruku widać, że wszystkie wartości GVIF^(1/(2*Df)) są znacznie poniżej powszechnie stosowanego progu 5 (a nawet 2), co oznacza, że zmienne nie są ze sobą współliniowe.
GVIF Df GVIF^(1/(2*Df))
log(exper + 1) 1.621670 1 1.273448
poly(educ, 2) 3.526653 2 1.370379
mtr 1.232930 1 1.110374
motheduc 1.233301 1 1.110541
wifecoll 3.246395 1 1.801775
IMR 1.617476 1 1.2718007 Weryfikacja hipotezy
Na wydruku równania wynikowego (poniżej) w modelu Heckmana widać, że współczynnik przy IMR nie jest statystycznie istotny. Oznacza to, że nie ma dowodów na selekcję próby, a model Heckmana można zastąpić modelem OLS. Rysunek 7 prezentuje wizualne porównanie wielkości współczynników w obu modelach.
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.715130 0.404935 6.7051 6.511e-11 ***
log(exper + 1) 0.171255 0.059932 2.8575 0.004482 **
poly(educ, 2)1 5.872026 1.170443 5.0169 7.765e-07 ***
poly(educ, 2)2 1.543472 0.666920 2.3143 0.021131 *
mtr -2.463867 0.423175 -5.8223 1.154e-08 ***
motheduc -0.018420 0.010683 -1.7242 0.085404 .
wifecollTRUE -0.250368 0.126681 -1.9764 0.048766 *
IMR -0.066835 0.152577 -0.4380 0.661582
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
8 Podsumowanie
Mimo zastosowania dwuetapowego modelu Heckmana, kluczowym wnioskiem jest brak statystycznie istotnego obciążenia wynikającego z selekcji próby.
Wnioski z równania płacowego są zgodne z teorią kapitału ludzkiego, potwierdzając istotny i pozytywny wpływ doświadczenia zawodowego oraz wykształcenia na wysokość stawki godzinowej, przy czym edukacja wykazuje nieliniowe, rosnące korzyści krańcowe. Wyższa krańcowa stopa podatkowa (mtr) negatywnie koreluje z płacą, co jest zgodne z oczekiwaniami ekonomicznymi. Pomimo pewnego odchylenia od normalności reszt, duży rozmiar próby pozwala na ufność w asymptotyczną wiarygodność wyników.
Wyniki te nie oznaczają jednak uniwersalnej przewagi OLS nad modelami selekcji. Potencjalne rozwinięcia badania powinny obejmować:
- analizę współczesnych danych, gdzie mechanizmy selekcji mogą być bardziej złożone ze względu na zmiany społeczno-ekonomiczne,
- zastosowanie alternatywnych zmiennych instrumentalnych w równaniu selekcji,
- porównanie wyników z metodą maksymalnej wiarygodności zamiast estymacji dwustopniowej, oraz
- rozszerzenie analizy o inne grupy demograficzne i rynki pracy.