Estymacja największej wiarogodności w modelu Gaussa-Markowa i kryteria informacyjne
Przypomnienie: gęstość wielowymiarowego rozkładu normalnego o wektorze wartości oczekiwanych µ i nieosobliwej macierzy kowariancji Σ zadana jest wzorem:
f (x) = 1
pdet(Σ) · (2π)n e−12(x−µ)0Σ−1(x−µ), x ∈ Rn. Dany jest model liniowy Gaussa-Markowa Y = Xβ + ε, w którym
Y =
Y1 Y2
...
Yn
, X =
x10 x11 . . . x1,k−1 x20 x21 . . . x2,k−1
. . . . xn0 xn1 . . . xn,k−1
, β =
β0 β1
...
βk−1
, ε =
ε1 ε2
...
εn
.
X jest macierzą deterministyczną i rz(X) = k < n. ε1, ε2, . . . , εd są niezależnymi zmiennymi losowymi o rozkładzie N (0, σ2). W takim razie ε ∼ N (0, σ2I), a co za tym idzie Y ∼ N (Xβ, σ2I). Poszukujemy estymatorów metody największej wiarogodności parametrów β (wektor) i σ2.
Gęstość wektora losowego Y ma postać:
f (y) = 1
pdet(σ2I) · (2π)n e−12(y−Xβ)0(σ2I)−1(y−Xβ)= 1
pσ2n· (2π)n e−2σ21 (y−Xβ)0(y−Xβ). Niech
L(β, σ2) = log f (Y ) = −n
2 log σ2−n
2log(2π) − 1
2σ2(Y − Xβ)0(Y − Xβ).
Niech S(β) =Pd
i=1(Yi− β0xi0− β1xi1− . . . − βk−1xi,k−1)2= (Y − Xβ)0(Y − Xβ). W takim razie L(β, σ2) = −n
2log σ2−n
2log(2π) − 1 2σ2S(β).
Niech S0(β) oznacza gradient funkcji S a S00(β) oznacza macierz drugich pochodnych funkcji S. Przypo- mnijmy, że S0(β) = −2XY + 2X0Xβ i S00(β) = 2X0X. W takim razie
∂
∂βL(β, σ2) = ∂L
∂β0
, ∂L
∂β1
, . . . , ∂L
∂βk−1
0
= − 1
2σ2S0(β) = 1
σ2(XY − X0Xβ)
∂
∂σ2L(β, σ2) = − n 2σ2 + 1
2σ4S(β).
W związku z tym
∂
∂βL(β, σ2) = 0 ⇔ β = (X0X)−1XY
∂
∂σ2L(β, σ2) = 0 ⇔ σ2= 1
nS(β) = 1
n||Y − Xβ||2.
Wobec powyższego kandydaci na estymatory metody największej wiarogodności mają postać:
β = (Xˆ 0X)−1XY, σˆ2= 1
n||Y − X ˆβ||2= 1
n||Y − ˆY ||2.
Należy jeszcze sprawdzić, czy w punkcie ( ˆβ, ˆσ2) = ( ˆβ0, ˆβ1, . . . , ˆβk−1, ˆσ2) funkcja L osiąga maksimum globalne. W tym celu wyznaczymy macierz drugich pochodnych tej funkcji. Ma ona postać:
∂2
∂β2L(β, σ2) ∂2
∂σ2∂βL(β, σ2)
∂2
∂σ2∂βL(β, σ2)
| {z }
k
∂2
∂(σ2)2L(β, σ2)
| {z }
1
=
− 1
2σ2S00(β) − 1
σ4(XY − X0Xβ)
− 1
σ4(XY − X0Xβ)0 n 2σ4− 1
σ6S(β)
1
Należy sprawdzić, czy macierz owa w punkcie ( ˆβ, ˆσ2) jest ściśle ujemnie określona. Zauważmy, że S( ˆβ) = nˆσ2 i XY − X0X ˆβ = XY − X0X(X0X)−1X0Y = XY − XY = 0. W takim razie
∂2
∂β2L( ˆβ, ˆσ2) ∂2
∂σ2∂βL( ˆβ, ˆσ2)
∂2
∂σ2∂βL( ˆβ, ˆσ2) ∂2
∂(σ2)2L( ˆβ, ˆσ2)
=
−1 ˆ
σ2X0X 0
0 − n
2ˆσ4
.
Ustalmy t = (t0, t1, . . . , tk−1, tk)0 ∈ Rk+1 takie że t 6= 0. Niech t−k= (t0, t1, . . . , tk−1)0.
t0
∂2
∂β2L( ˆβ, ˆσ2) ∂2
∂σ2∂βL( ˆβ, ˆσ2)
∂2
∂σ2∂βL( ˆβ, ˆσ2) ∂2
∂(σ2)2L( ˆβ, ˆσ2)
t =
t0−k tk
− 1 ˆ
σ2X0X 0
0 − n
2ˆσ4
t−k
tk
=
= t0−k·
−1 ˆ σ2X0X
· t−k+ tk·
− n 2σˆ4
· tk = − 1 ˆ
σ2 · t0−kX0Xt−k− n 2ˆσ2 · t2k=
= − 1 ˆ
σ2 · (X · t−k)0(X · t−k) − n
2ˆσ2· t2k = − 1 ˆ
σ2 · ||X · t−k||2− n 2ˆσ2 · t2k
Ponieważ rz(X) = k < n, więc X · t−k 6= 0, a zatem ||X · t−k||2 > 0. W takim razie wartość formy kwadratowej powyżej jest ujemna. Oznacza to, że ta forma kwadratowa jest ściśle ujemnie określona.
Wobec tego w punkcie ( ˆβ, ˆσ2) funkcja L osiąga maksimum lokalne. Jest to jedyny punkt krytyczny funkcji L, a zatem jest to maksimum globalne. W związku z tym ˆβ i ˆσ2 są odpowiednio estymatorami metody największej wiarogodności parametrów β i σ2.
Estymator metody największej wiarogodności wektora parametrów β jest estymatorem metody naj- mniejszych kwadratów tegoż parametru. Jest więc estymatorem nieobciążonym.
Estymator metody największej wiarogodności parametru σ2 różni się od estymatora dotychczas po- znanego, toteż od tego momentu przyjmijmy następujące oznaczenia (L od likelihood – wiarogodność):
ˆ
σ2= ||Y − ˆY ||2
n − k , σˆL2 = ||Y − ˆY ||2
n .
Wiemy, że E ˆσ2= σ2, a zatem E ˆσ2L=n−kn E ˆσ2=n−kn σ2. Estymator metody największej wiarogodności parametru σ2nie jest więc estymatorem nieobciążonym.
Przyjmijmy następujące oznaczenie: ˆL = supβ∈Rk,σ2>0L(β, σ2). Zgodnie z definicją estymatora naj- większej wiarogodności ˆL = L( ˆβ, ˆσ2L). Wobec tego
L = −ˆ n
2log ˆσ2L−n
2 log(2π) − 1
2ˆσ2S( ˆβ) = −n
2 log||Y − ˆY ||2
n −n
2log(2π) − n 2. Definicja 1. Kryterium informacyjnym Akaike nazywamy wyrażenie AIC = −2 ˆL + 2k.
Definicja 2. Bayesowskim kryterium informacyjnym (lub kryterium informacyjnym Schwartza) nazy- wamy wyrażenie BIC = −2 ˆL + k log n.
Zgodnie z powyższymi definicjami
AIC = n log||Y − ˆY ||2
n + n log(2π) + n + 2k, BIC = n log||Y − ˆY ||2
n + n log(2π) + n + k log n.
Ponieważ jednak n log(2π) + n nie zależy od doboru zmiennych, więc najczęściej przyjmuje się dla uprosz- czenia, że
AIC = n log||Y − ˆY ||2
n + 2k, BIC = n log||Y − ˆY ||2
n + k log n.
W zależności od tego, do czego ma służyć model i jakie są zależności między zmiennymi, inne kry- terium wyboru modelu okaże się dobre. Jeżeli pożądane są dobre właściwości predykcyjne, odpowiednie
2
jest kryterium AIC, które wybiera duże modele, ale o dobrych właściwościach predykcyjnych. Jeżeli szu- kamy modelu opisującego prawdziwą zależność, często stosowane jest kryterium BIC, które jest zgodne (czyli asymptotycznie z prawdopodobieństwem 1 wybiera dobry model), ale wybiera mniej zmiennych niż AIC (Przemysław Biecek, Analiza danych z programem R. Modele liniowe z efektami stałymi, losowymi i mieszanymi, Wydawnictwo Naukowe PWN, Warszawa 2011, str. 122).
Korzystanie z kryteriów informacyjnych wymusza zachowanie założeń modelu Gaussa-Markowa, w szcze- gólności założenia o normalności błędów, gdyż kryteria informacyjne oparte są na maksymalizacji funkcji wiarogodności a ta zależy od łącznego rozkładu błędów.
Dobór zmiennych za pomocą kryteriów informacyjnych w R realizuje funkcja step. Jeśli pozosta- wimy domyślą wartość argumentu k=2, będziemy się posługiwali AIC, jeśli zaś przyjmiemy k=log(n), to będziemy się posługiwali kryterium BIC. W argumencie object umieszczamy model początkowy. W ar- gumencie scope podajemy listę zmiennych-kandydatek do włączenia do modelu. Wartością argumentu scope możemy także uczynić listę składającą się z dwóch elementów: lower i upper będących odpowied- nio najuboższym i najbogatszym modelem, jaki chcemy rozważać. Argument direction określa kierunek przeszukiwania modeli.
Aby wykonać jedną z czterech procedur krokowych dotąd poznanych, należy postępować w sposób następujący:
jeśli jako wartość argumentu object przyjmiemy model tylko z wyrazem wolnym (tzn. lm(Y~1)), jako wartość argumentu scope tyldę i prawą stronę formuły definiującej model pełny, a ponadto direction="forward", to zostanie przeprowadzona selekcja postępująca,
jeśli jako wartość argumentu object przyjmiemy model pełny, to zostanie przeprowadzona elimi- nacja wsteczna (automatycznie zostanie przyjęte direction="backward"),
jeśli jako wartość argumentu object przyjmiemy model tylko z wyrazem wolnym (tzn. lm(Y~1)), jako wartość argumentu scope tyldę i prawą stronę formuły definiującej model pełny, a ponadto direction="both", to zostanie przeprowadzona regresja krokowa postępująca,
jeśli jako wartość argumentu object przyjmiemy model pełny a ponadto direction="both", to zostanie przeprowadzona regresja krokowa wsteczna.
Należy zaznaczyć, że regresja krokowa postępująca i regresja krokowa wsteczna zostaną przeprowadzone z drobną różnicą w stosunku do tego, co zostało uprzednio powiedziane.
Transformacje Boxa-Coxa
Brak normalności błędów (ε) może przybrać postać złego dobrania skali, w której dokonywane są po- miary. Jeśli realizacje Y1, Y2, . . . , Yn są dodatnie, to rozwiązaniem tego problemu może być odpowied- nie przekształcenie zmiennej zależnej. Popularną rodziną przekształceń stosowanych w takich sytuacjach są transformacje Boxa-Coxa:
Z(λ)=
Zλ− 1
λ , λ 6= 0 log Z, λ = 0
Wybór λ może być dokonany za pomocą metody największej wiarogodności. Przebiega on zwykle w następujących krokach:
wybiera się pewien zbiór Λ = {λ0, λ1, . . . , λl}, najczęściej o elementach postaci λi= λ0+ i∆λ,
za pomocą każdego λ ∈ Λ dokonuje się transformacji Boxa-Coxa zmiennej zależnej, otrzymując Y(λ)= (Y1(λ), Y2(λ), . . . , Yn(λ))0, a następnie metodą największej wiarogodności dokonuje się estymacji parametrów w modelu Gaussa-Markowa Y(λ) = Xβ + ε. Następnie oblicza się maksimum funkcji wiarogodności dla tego modelu, co oznaczamy jako ˆL(λ),
wybiera się λ takie że ˆL(λ) = max{ ˆL(λ0), ˆL(λ1), . . . , ˆL(λl)}.
Warto zwrócić uwagę na to, że dla λ > 1 transformacje Boxa-Coxa są funkcjami wypukłymi, podczas gdy dla λ < 1 transformacje Boxa-Coxa są funkcjami wklęsłymi.
Obliczenia pozwalające na dokonanie wyboru wartości parametru λ w R realizuje funkcja boxcox z pakietu MASS.
3
