Estymacja parametrów modelu

Jak wspomniano w podrozdziale 3.2, do szacowania parametrów modelu stosuje się zapis przedstawiony przy pomocy równania (3.1) w postaci macierzowej. W podrozdziałach 3.2 oraz 3.3 opisano formę, jaką przyjmują składniki równania (3.1), w kolejnych etapach rozszerzania modelu wielopoziomowego, zarówno klasycznego jak i z dwoma kryteriami grupowania,. Chociaż modele te różniły się pod względem strategii ich budowy, parametry modelu szacowane są dokładnie w taki sam sposób dla obu przypadków (por. Bates, 2010).

W celu estymacji współczynników modelu wielopoziomowego najczęściej stosuje się metodę największej wiarogodności (ML maximum likelihood) bądź resztową metodę największej wiarogodności (REML lub RML restricted/residual maximum likelihood). W rzeciwieństwie do klasycznego modelu regresji liniowej, estymatory resztowej największej wiarogodności, w przypadku modeli wielopoziomowych mogą być obciążone. Uważa się jednak, że obciążenie to w większości wypadków jest mniejsze, niż przy zastosowaniu metody największej wiarogodności (por. Biecek, 2011). Z drugiej strony, estymatory parametrów szacowane przy pomocy resztowej metody największej wiarogodności, pomimo takiej nazwy, nie charakteryzują się największą wiarogodnością. Ponadto modele szacowane tą metodą nie mogą być traktowane jako rozwinięcia modeli „uboższych”, jeżeli modyfikacja dotyczyła części stałej modelu, zatem nie można ich porównywać przy pomocy testu ilorazu wiarogodności (por. Frątczak i in., 2009; Bates, 2010; Hox, 2010; Biecek, 2011). Jest tak dlatego, że estymatory REML obliczane są po przeniesieniu części stałej modelu na lewą stronę równia (3.1). W związku z tym przy każdej modyfikacji części stałej modelu, zmienia się część objaśniana modelu, Zatem kolejne modele nie mogą być traktowane jako rozszerzenia wcześniejszych. Ponadto

109

w literaturze podaje się, że wpływ wyboru metody ML lub REML na otrzymane oszacowania dla modeli wielopoziomowych jest mniejszy niż się powszechnie uważa (por. Bates, 2010; Hox, 2010; Kreft, de Leeuw, 1998). W związku z powyższymi uwagami w następnym, czwartym rozdziale prezentowanej pracy wykorzystano, przy konstruowaniu modeli metodę największej wiarogodności (ML).

Ponadto w dostępnym oprogramowaniu estymacja parametrów modelu, zarówno w przypadku klasycznej jak i resztowej metody największej wiarogodności, przeprowadzona może zostać przy pomocy różnych algorytmów. Niektóre z nich to algorytm Newtona-Raphsona (wykorzystywana w programie SAS, por. Wolfinger, Tobias, Sall, 1994) oraz metoda z wykorzystaniem operacji na macierzach rzadkich (w pakiecie lme4 programu R, por. Bates, 2010). W prezentowanej pracy obliczenia wykonane zostały przy pomocy pakietu lme4 programu R. W dalszej części tego podrozdziału opisano algorytm szacowania parametrów modelu stosowany w tym pakiecie.

ESTYMATORY NAJWIĘKSZEJ WIAROGODNOŚCI

Rozważmy ogólny wielopoziomowy model liniowy opisany równaniem (3.1). Jak zauważono wcześniej, przy pomocy tego równania zapisać można dowolny model wielopoziomowy, zarówno klasyczny jak i z dwoma kryteriami grupowania. Ponieważ estymacja parametrów modelu przebiegać będzie dla poszczególnych elementów tego zapisu, został on powtórzony, ponadto określone zostały wymiary poszczególnych składowych

m = no + pq + r , (3.1)

gdzie:

Y – wektor zmiennej objaśnianej o wymiarze n,

X – macierz zmiennych objaśniających dla efektów stałych modelu

o wymiarze n x p,

β – wektor współczynników kierunkowych o wymiarze p,

Z – macierz zmiennych objaśniających dla efektów losowych modelu

o wymiarze n x q,

e – wektor efektów losowych o wymiarze q, s~N_†vw; x) ,

r – składnik resztowy modelu o wymiarze n, y~N_Óvw; σ{_Ô) .

Należy zwrócić uwagę, że n, podobnie jak w poprzednich podrozdziałach, określa liczbę jednostek pierwszego poziomu, jednak p i q są odpowiednio wymiarami macierzy

110

efektów stałych oraz losowych. Nie można utożsamiać ich z liczbą zmiennych objaśniających czy liczbą poziomów modelu. W podrozdziałach 3.2 oraz 3.3 przedstawiono te wymiary dla modeli opisywanych w kolejnych etapach.

W równaniu (3.1) znane są wartości: wektora zmiennej objaśnianej Y (tylko z próby uczącej), macierzy zmiennych objaśniających dla efektów stałych X oraz macierzy Z. Celem rozważań jest estymacja wektora β oraz macierzy kowariancji wektora efektów losowych e, czyli x. Parametry te oszacowane zostaną przy pomocy metody największej

wiarogodności.

W pierwszej kolejności, przeanalizowana zostanie postać macierzy kowariancji x.

Załóżmy, że omawiany model charakteryzuje się strukturą P poziomową. Ponadto ze względu na kolejne poziomy uwzględnia się ki (dla i=1,…,P) efektów losowych (wliczając w to efekt losowy wyrazu wolnego oraz współczynników kierunkowych przy zmiennych z niższych poziomów). Tak więc dla każdego poziomu należy oszacować ki wartości wariancji oraz kix(ki-1)/2 kowariancji, czyli w sumie kix(ki+1)/2 wartości w macierzy kowariancji, związanych z danym poziomem. Ponieważ, na każdym poziomie jest określona liczba jednostek, z których efekty losowe dla każdej z nich charakteryzują się takim samym rozkładem normalnym, opisane powyżej wartości pojawiają się w macierzy x tyle razy, z ilu jednostek składa się dany poziom (dwukrotnie

więcej w przypadku kowariancji). Ponieważ zakłada się, że nie występują korelacje losowe pomiędzy efektami losowymi z różnych poziomów, pozostałe wartości macierzy

x są zerami. Tak więc, przy optymalnej kolejności kolumn w macierzy Z (którą przyjąć

można bez straty ogólności, a która w procedurze lmer uzyskiwana jest przy pomocy prostej transformacji), macierz x jest macierzą blokową diagonalną, wokół diagonali

której pojawiają się kolejno macierze kowariancji z kolejnych poziomów, każda tyle razy, ile jednostek określonych jest na danym poziomie. Zatem przy potencjalnie dużej liczbie elementów macierzy x (q2

), w celu jednoznacznego jej określenia konieczne jest oszacowanie relatywnie małej liczby jej elementów ∑ R£ _vR_+ 1)/2

e ^{. Wektor tych}

parametrów oznaczać będziemy θ. W związku z tym również macierz kowariancji x

powinna być oznaczana jako x_Õ.

O zmiennej objaśnianej zakłada się, że ma rozkład normalny, zatem jeżeli przyjąć konkretną wartość wektora efektów losowych, e, zmienna objaśniana charakteryzuje się następującym warunkowym rozkładem normalnym:

111

vÖ|q)~_×vno + pq; σ{_Ô) . (3.25)

Z kolei wektor efektów losowych charakteryzuje się następującym rozkładem normalnym:

s~N_†vw; x_Õ) . (3.26)

O macierzy x_Õ, jako o macierzy kowariancji, wiadomo, że jest macierzą symetryczną, ponadto jest dodatnio półokreślona, tzn:

⋀_Û∈ÜÝ/^‡_ÙÚx_ÕÙ ≥ w .

Niech macierz Þ_Õ będzie macierzą trójkątną dolną, spełniającą warunek:

x_Õ = σÞ_ÕÞ_Õ^Ú . (3.27)

Wiadomo, że macierz spełniająca powyższy warunek istnieje, ponieważ macierz x_Õ jest

macierzą symetryczną.

Następnie wektor efektów losowych e przedstawiony zostaje jako:

q = Þ_Õß . (3.28)

Łatwo można wykazać, że wektor u charakteryzuje się następującym rozkładem normalnym:

ß~_à4w; σ{_á5 . (3.29)

Jednocześnie równanie (3.1) oraz rozkład warunkowy zmiennej objaśnianej zapisać można:

m = no + pÞ_Õß + r , (3.30)

vÖ|ß)~×vno + pÞÕß; σ{_Ô) . (3.31) W wyniku powyższego zapisu, obiektem zainteresowania stała się macierz Þ_Õ zamiast

macierzy x_Õ oraz wektor u zamiast wektora e. Zamiana taka jest korzystna ze względu

na prostszy rozkład wektora u.

Aby zastosować metodę największej wiarogodności konieczna jest znajomość rozkładu łącznego wektorów u oraz y, w celu wyznaczenia funkcji wiarogodności opisanej wzorem:

âvSÕ, o, ã|Ö_äÛå) = æ ç_ÜÝ è,éQê . (3.32) Korzystając z twierdzenia Bayesa, odpowiednią gęstość warunkową można zapisać w postaci:

112

Ponieważ znane są rozkłady vÖ|ß) oraz ß (por. odpowiednio wzory (3.31) oraz (3.29)),

w celu wyznaczenia gęstości łącznej wyznaczane są odpowiednie gęstości (por. Krzyśko, 2009): çSé|è = v2ë)^¯ì|σ{_Ô|^¯@íG î−^_vÖ − no − pÞ_Õß)vσ{_Ô)¯vÖ − no − pÞ_Õß)ï = v2ë)^¯ìvσ)^¯ì@íG î−_ð^vÖ − no − pÞ_Õß)vÖ − no − pÞ_Õß)ï = v2ëσ)¯^ì@íG î−_ð^vÖ − no − pÞ_Õß)vÖ − no − pÞ_Õß)ï , (3.34) ç_ß = v2ë)^¯Ý|σ{_Ô|¯^Ý@íG î−^_ß4σ{_á5^¯ßï = v2ë)^¯Ývσ)¯^Ý@íG î−_ð^ßßï = v2ëσ)^¯Ý@íG î−_ð^ßßï . (3.35) Podstawiając tak wyznaczone gęstości do równania (3.33) wyznaczono gęstość łączną jako:

ç_ß,Ö = v2ëσ)^¯ìñÝ @íG ò−_ð^ÎvÖ − no − pÞÕß)vÖ − no − pÞÕß) + ßßÐó. (3.36)

Niech ßô będzie wartością wektora u maksymalizującą powyższą funkcję gęstości łącznej.

Ponadto niech Ö_äÛå oznacza wektor wartości (stały), które zmienna objaśniana przyjmuje w próbie.

ßô = LPõ max_ßçSß|Övß|Ö_äÛå) = LPõ max_ßç_ß,Övß, Ö_äÛå) = LPõ max_ßçSÖ|ßvÖ_äÛå|ß)ç_ßvß).

(3.37) Powyższa równość jest prawdziwa, ponieważ jest wartością funkcji w punkcie, a więc ç_évÖ_äÛå) jest stałą.

Następnie wyznaczany jest logarytm naturalny łącznej gęstości pomnożony przez (-2):

−2 ln ç_ß,Övß, Ö_äÛå) = −2 ln öv2ëσ)^¯ìñÝ @íG ò−_ð^ÎvÖ − no − pÞ_Õß)vÖ − no − pÞ_Õß) + ßßÐó÷ = v) + ˜) lnv2ëσ) +_ð^vvÖ − no − pÞ_Õß)vÖ − no − pÞ_Õß) +

ßß) . (3.38)

Ponieważ logarytm naturalny jest funkcją rosnącą w swojej dziedzinie, powyższa funkcja jest malejąca. W związku z tym argument, dla którego funkcja ç_è,é osiąga maksimum, to ten sam argument, dla którego funkcja −2 ln ç_è,é osiąga minimum. Ponadto we wzorze (3.38) od u zależy wyłącznie drugi składnik sumy, a dokładnie druga część iloczynu wchodzącego w jego skład. W związku z tym równanie (3.37) można przekształcić do postaci:

113

ßô = LPõ min_ßvvÖ − no − pÞ_Õß)vÖ − no − pÞ_Õß) + ßß) . (3.39) W powyższym wzorze wyrażenie, z którego szukane jest minimum, nazywane jest sumą kwadratów reszt z karą (PRSS penalized residual sum of squares) ). Z kolei wartość u, dla którego to minimum jest uzyskiwane, metodą najmniejszych kwadratów z karą (PLS penalized least squares solution) (por. Bates, 2013c).

Niech:

P_Õ,o = min_ßvvÖ − no + pÞÕß)vÖ − no + pÞÕß) + ßß) (3.40)

P_Õ= min_ß,ovvÖ − no + pÞ_Õß)vÖ − no + pÞ_Õß) + ßß) (3.41) Obydwa powyższe minima wyznaczone mogą zostać bez użycia procedury iteracyjnej. Ponadto wartość wektora o (zależna od parametru ø) dla której osiągane jest minimum

(3.40) oznaczana będzie dalej jako où_Õ.

W celu rozwiązania równania (3.39) minimalizowaną sumę przedstawić można jako:

ßô = LPõ min_ßîÖ − no − pÞÕß ß ^ï  îÖ − no − pÞÕß ß ^{ï =} LPõ min_ßúûîÖ − no_{w ï − û}^pÞ_ü ^Õ ý þßþ^ûîÖ − no_{w ï − û}^pÞ_ü ^Õ ý þßþ
. (3.42) Dla tak zapisanego problemy rozwiązanie znaleźć można przy pomocy klasycznej metody najmniejszych kwadratów. Jest to wartość spełniająca następujące równanie:

ûpÞ_ü ^Õ ý þ^ûpÞ_ü ^Õ ý þ ßô = ûpÞ_ü ^Õ ý þ^îÖ − no_{w ï} , 4vpÞ_Õ)ÚpÞ_Õ+ ü_ý5ßô = pÞ_ÕvÖ − no) , 4Þ_ÕÚpÚpÞ_Õ+ ü_ý5ßô = pÞ_ÕvÖ − no) ßô = 4Þ_Õ^ÚpÚpÞ_Õ+ ü_ý5^¯pÞ_ÕvÖ − no). (3.43) Rozwiązanie tego równania nie jest oczywiście wyznaczane przez obliczanie macierzy odwrotnej, ale przy pomocy metody Choleskiego. Tak więc, przy pomocy dekompozycji Choleskiego dla macierzy rzadkich wyznaczana jest macierz trójkątna dolna

_Õ spełniająca warunek:

114 Podstawiając (3.44) do wzoru (3.43) otrzymujemy:

_Õ_Õ^Úßô = pÞ_ÕvÖ − no) . (3.45) Przy pomocy metody Choleskiego możliwe jest rozwiązanie powyższego równania. Następnie rozwiązanie to podstawione zostanie do równania (3.40):

P_Õ,o = vÖ − no + pÞ_Õßô)vÖ − no + pÞ_Õßô) + ßôßô . (3.46) Ponadto dla dowolnego u prawdziwa jest równość:

vÖ − no + pÞ_Õß)vÖ − no + pÞ_Õß) + ßß = P_Õ,o + ¦_Õ^Úvß − ßô)§^¦_Õ^Úvß − ßô)§ .

(3.47) Podstawiając powyższe równanie do równania (3.36) otrzymujemy:

ç_è,é = v2ëσ)¯^ìñÝ @íG −_ð^ûP_Õ,o + ¦_Õ^Úvß − ßô)§^¦_Õ^Úvß − ßô)§þ÷ . (3.48) Na potrzeby dalszych rozważań zdefiniowano następującą zmienną pomocniczą z:

I =^Õ^Úv߯ßô)

 ^{. (3.49)}

Tak zdefiniowana zmienna ma następującą własność:

^_ß = ^Õ^Ú

 =^|Õ|

 ^{. (3.50)}

W tym momencie możliwe jest wyznaczenie wiarogodności:

âvSÕ, o, ã|Ö_Ù) = æ ç_ÜÝ è,éQß = = æ v2ëσ_ÜÝ )¯^ìñÝ @íG −_ð^ûP_Õ,o + ¦_Õ^Úvß − ßô)§^¦_Õ^Úvß − ßô)§þ÷Qß == v2ëσ)^¯ì@íG −Õ,o ð æ v2ë)_ÜÝ ^¯Ý@íG −_ð^¦_Õ^Úvß − ßô)§^¦_Õ^Úvß − ßô)§÷|Õ| |Õ| ß ð = ^{I =} Õ^Úv߯ßô)  _ß^ =^|Õ|   = v2ëσ)^¯ì@íG −Õ,o ð |_Õ|¯æ v2ë)_ÜÝ ^¯Ý@íG ¦−^C_^§Q = = v2ëσ)¯^ì@íG −Õ,o ð |_Õ|¯1 = v2ëσ)¯^ì@íG −Õ,o ð |_Õ|¯ _(3.51) Podwojony logarytm naturalny wiarogodności nazywany jest dewiacją i oznaczany:

QvSÕ, o, ã|Ö_Ù) = −2 ln âvSÕ, o, ã|Ö_Ù) = ) lnv2ëσ) + 2 ln|_Õ| +Õ,o

115

Estymatorami największej wiarogodności nazywamy te oszacowania parametrów, które minimalizują powyższą dewiację. Niech où_Õ oznacza tę wartość o, która

minimalizuje sumę kwadratów reszt z karą (PRSS) (por. Bates, 2010) dla ustalonego Õ

oraz niech P_Õ _{oznacza sumę kwadratów reszt z karą dla} o = où_Õ. Wtedy:

QvSÕ, ã|Ö_Ù) = ) lnv2ëσ) + 2 ln|_Õ| +Õ

ð . (3.53)

Ponadto niech σ_Õ będzie wartością σ_{, która minimalizuje dewiację dla ustalonego} P_Õ_. Wtedy σ_{Õ =Õ}

× ^oraz:

QvSÕ|Ö_Ù) = ) ln ¦^Õ

× § + 2 ln|_Õ| + ) = 2 ln|_Õ| + ) ¦1 + ln ¦^Õ

× §§ . (3.54) W celu wyznaczenia ßô oraz où_Õ zastosować można następujący zapis:

û ßô_où Õþ = LPõ min_ß,oîÖ − pÞÕß − no ß ^ï  îÖ − pÞÕß − no ß ^{ï =} LPõ min_ßúûîÖwï − û^pÞü_ý^{Õ n}wþ îoïþ^ß  ûîÖwï − û^pÞü_ý^{Õ n}wþ î^ßoïþ
. (3.55) Zgodnie z metodą najmniejszych kwadratów rozwiązaniem powyższego problemu są oszacowania spełniające następujący warunek:

û^pÞ_ü ^{Õ n} ý wþ  û^pÞ_ü ^{Õ n} ý wþ ûoù^ßô_Õ^{þ = ûpÞ}ü_ý^{Õ n}wþ  îÖwï , (3.56)

co na drodze elementarnych przekształceń daje:

Þ^{ÕÚpÚpÞÕ+ üý ÞÕÚpÚn} nÚpÞ_Õ nÚn ^{ û ßôoù}Õþ = ÞÕ^ÚpÚÖ nÚÖ ^{ ,} û ßô_où Õþ = Þ^{ÕÚpÚpÞÕ+ üý ÞÕÚpÚn} nÚpÞ_Õ nÚn ^ ¯ ÞÕ^ÚpÚÖ nÚÖ ^{. (3.57)}

Podobnie jak wcześniej, rozwiązanie powyższego problemu znalezione zostaje przy pomocy metody Choleskiego. Przeprowadzana zostaje następująca dekompozycja Choleskiego (dla macierzy rzadkich):

Þ^{ÕÚpÚpÞÕ+ üý ÞÕÚpÚn}

nÚpÞ_Õ nÚn ^{ = û }pn^Õ ^w_n^{þ Õ}

Ú _pn^Ú w _n^{Ú .}

116

W powyższym zapisie _Õ, tak jak wcześniej, oznacza trójkątną dolną macierz

wyznaczoną przy pomocy dekompozycji Choleskiego spełniającą warunek określony wzorem (3.44). Macierze _pn oraz _n są odpowiednio wymiaru pxq oraz pxp spełniają poniższe warunki:

_Õ_pnÚ = Þ_ÕÚpÚn (3.58)

oraz

_n_n^Ú = nÚn − _pn_pn^Ú . (3.59) Oszacowania parametrów równania (3.57) wyznaczone zostaną więc przy pomocy metody Choleskiego na podstawie równania:

û _^{Õ w}

pn _n^{þ Õ}

Ú _pn^Ú

w _n^{Ú û ßôoù}Õþ = Þ^Õ_n^Ú_Ú^p_Ö^ÚÖ . (3.60) Powyższe macierze uzyskane z dekompozycji Choleskiego zostały zapisane w formie blokowej, ponieważ jeden z wyznaczonych bloków posłuży w dalszych rozważaniach do wyznaczenia estymatorów resztowej metody największej wiarogodności.

ESTYMATORY RESZTOWEJ METODY NAJWIĘKSZEJ WIAROGODNOŚCI Funkcję resztowej największej wiarogodności zapisać można:

â_ÜvSÕ, ã|Ö_Ù) = æ âvSÕ, o, ã|Ö_Ü Ù)Qo . (3.61) Resztowa dewiacja wyznaczana jest zatem jako:

Q_ÜvSÕ, ã|Ö_Ù) = −2 ln â_ÜvSÕ, ã|Ö_Ù) = −2 ln æ âvSÕ, o, ã|Ö_Ü Ù)Qo . (3.62) Sformułować można następującą równość:

P_Õ,o = P_Õ+ ¦_n^Ú4o − où_Õ5§^¦_n^Ú4o − où_Õ5§ . (3.63) Podstawiając równanie (3.51) do równania (3.61), resztową wiarogodność można przedstawić jako:

â_ÜvSÕ, ã|ÖÙ) = æ v2ëσ)^¯ì@íG −Õ,o

ð |_Õ|¯Qo

Ü . (3.64)

117 â_ÜvSÕ, ã|Ö_Ù) = æ v2ëσ)^¯ì@íG −Õ.¦n^Ú4o¯oùÕ5§^C¦n^Ú4o¯oùÕ5§ ð  |_Õ|¯Qo Ü = @íG ¦− Õ ð§ v2ëσ)^¯ì± |Õ|¯æ v2ëσ)^¯@íG −^{¦nÚ4o¯oùÕ5§} C ¦n^Ú4o¯oùÕ5§ ð ^|n| |_n|^oð‚ Ü = ^{I =} n^Ú4o¯oùÕ5  _o^ =^|n| ð‚  = @íG ¦− Õ ð§ v2ëσ)¯^ì± |_Õ|¯|_n|¯æ v2ëσ_Ü )¯^@íG ¦−^C_^§ Q = @íG ¦− Õ ð§ v2ëσ)¯^ì± |_Õ|¯|_n|¯_{. (3.65)} Zatem resztowa dewiacja wynosi:

Q_ÜvSÕ, ã|Ö_Ù) = −2 ln ¦@íG ¦− Õ

ð§ v2ëσ)¯^ì± |_Õ|¯|_n|¯§ = v) − G& lnv2ëσ& + 2 lnv|_Õ||_n|& +Õ

ð . (3.66)

Estymatorem resztowej największej wiarogodności parametru ã jest σ_{Õ = Õ} ׯ¤^. Zatem: Q_ÜvSÕ|Ö_Ù& = v) − G& ln ¦^Õ ×¯¤§ + 2 lnv|_Õ||_n|& + ) − G = =

2

lnv|_Õ||_n|& + v) − G& ¦1 + ln ¦^Õ ×¯¤§§ . (3.67)

Oszacowanie parametru Õ, Õù_ wyznaczane jest jako wartość minimalizująca resztową

dewiację, zaś oszacowaniem ã jest σ_Õù

. Jako estymator o przyjmuje się où_Õù, gdzie où_Õ

jest estymatorem wyznaczonym metodą największej wiarogodności dla danego Õ.

W dokumencie Delimitacja lokalnych rynków pracy w Polsce z wykorzystaniem modelowania wielopoziomowego (Stron 108-117)