REML Pavla Dokoupilová
Restringovaná (reziduální) metoda maximální věrohodnosti
- REML
Jednou z vlastností maximálně věrohodného odhadu varianční složky je, že nebere v
úvahu stupně volnosti potřebné při odhadu pevného efektu. Výsledný odhad varianční
složky pak není nestranný (až na speciální případy).
Tento problém je odstraněn při použití restringované (reziduální) metody maximální
věrohodnosti, jejíž základní myšlenkou je odhad varianční složky založený na reziduích
získaných při odhadu pevných efektů metodou nejmenších čtverců.
Věrohodnostní funkce
Uvažujme n rozměrný normálně rozložený vektor y se střední hodnotou Xβ a kovarianční
maticí V (θ), tedy y ∼ Nn(Xβ, V (θ)). Předpokládejme, že X je n × p rozměrná matice s
hodností p a V je n × n rozměrná pozitivně definitní matice.
REML je založeno na lineární kombinaci A y, která neobsahuje žádné pevné efekty,
tj. A X = 0, kde A je n × (n − p) rozměrná matice. Přičemž každá lineární funkce L y
taková, že L X = 0, může být vyjádřena jako funkce A y. Tedy A y ∼ N(n−p)(0, A V A) a
LREML = (2π)−(n−p)/2
|A V A|−1/2
e− 1
2
y A(A V A)−1A y
. (1)
Běžně je však REML věrohodnostní funkce uváděná ve tvaru
LREML = (Konst.)|V |−1/2
|X V −1
X|−1/2
e− 1
2
(y−X ˆβ) V −1(y−X ˆβ)
, (2)
kde ˆβ = (X V −1
X)−1
X V −1
y. Dále bude ukázáno, že tyto zápisy jsou ekvivalentní.
Pomocná tvrzení a označení
1. Pokud M a N jsou m × m matice, pak |MN| = |M||N|.
2. Uvažujme regulární symetrickou matici
M =
M11 M12
M12 M22
.
Pak
|M| = |M22||M11 − M12M−1
22 M12|.
Důkaz. Platí
I −M12M−1
22
0 I
M11 M12
M12 M22
=
M11 − M12M−1
22 M12 0
M12 M22
,
I −M12M−1
22
0 I
= 1 a
M11 − M12M−1
22 M12 0
M12 M22
= |M22||M11−M12M−1
22 M12|,
čímž dostáváme tvrzení.
1
REML Pavla Dokoupilová
3.
∂ln|B(t)|
∂t
= Tr B−1 ∂B(t)
∂t
4.
∂B(t)−1
∂t
= −B−1 ∂B(t)
∂t
B−1
Tvrzení 1. Uvažujme matice V , X a A výše uvedených vlastností, pak
V −1
= V −1
X(X V −1
X)−1
X V −1
+ A(A V A)−1
A .
Důkaz. Platí
(V −1
X, A)−1
V −1
(V −1
X, A) −1
= [(V −1
X, A) V (V −1
X, A)]−1
=
X V −1
X 0
0 A V A
−1
=
(X V −1
X)−1
0
0 (A V A)−1
Dále
V −1
= (V −1
X, A)[(V −1
X, A) V (V −1
X, A)]−1
(V −1
X, A)
= (V −1
X, A)
(X V −1
X)−1
0
0 (A V A)−1 (V −1
X, A)
= V −1
X(X V −1
X)−1
X V −1
+ A(A V A)−1
A .
Tedy
y A(A V A)−1
A y = y [V −1
− V −1
X(X V −1
X)−1
X V −1
]y
= y V −1
y − y V −1
X(X V −1
X)−1
X V −1
y
= y V −1
y − y V −1
X(X V −1
X)−1
(X V −1
X)(X V −1
X)−1
X V −1
y
= (y − y V −1
X(X V −1
X)−1
X )V −1
(y − X(X V −1
X)−1
X V −1
y)
= (y − X(X V −1
X)−1
X V −1
y) V −1
(y − X(X V −1
X)−1
X V −1
y)
= (y − X ˆβ) V −1
(y − X ˆβ), (3)
čímž jsme ukázali ekvivalenci exponentů výrazů (1) a (2).
Označme dále W = V −1
− V −1
X(X V −1
X)−1
X V −1
= A(A V A)−1
A a derivaci Vi =
∂V/∂θi. Úpravami obdržíme
∂ ln |A V A|
∂θi
= Tr (A V A)−1 ∂A V A
∂θi
= Tr (A V A)−1
A ViA =
= Tr A(A V A)−1
A Vi = Tr[WVi] (4)
2
REML Pavla Dokoupilová
a dále
∂[ln |V | + ln |X V −1
X|
∂θi
= Tr[V −1
Vi] − Tr (XV −1
X)−1 ∂(X V −1
X)
∂θi
= Tr[V −1
Vi] − Tr (X V −1
X)−1
X
∂V −1
∂θi
X
= Tr[V −1
Vi] − Tr (X V −1
X)−1
X V −1
ViV −1
X
= Tr[V −1
Vi] − Tr V −1
X(X V −1
X)−1
X V −1
Vi
= Tr[(V −1
− V −1
X(X V −1
X)−1
X V −1
)Vi]
= Tr[WVi]. (5)
Z výrazů (4) a (5) vyplývá
ln |A V A| = (Konst.) + ln |V | + ln |X V −1
X|
|A V A| = (Konst.)|V ||X V −1
X|. (6)
Celkem jsme tedy pomocí (3) a (6) dokázali ekvivalenci mezi výrazy (1) a (2).
Restringované maximálně věrohodné odhady
Uvažujme dále náhodný vektor y ∼ N(Xβ, V ) ve tvaru y = Xβ + s+1
i=1 Zibi, kde cov(y) =
s+1
i=1 σ2
i ZiZi := V . Potom lineární kombinace A y bude tvaru A y = A Xβ + s+1
i=1 A Zibi.
ML odhady
Jak je ukázáno například v [2], maximálně věrohodné rovnice jsou
∂l
∂β
= X V −1
y − X V −1
Xβ = 0 (7)
∂l
∂σ2
i
= −
1
2
Tr(V −1
ZiZi) +
1
2
(y − Xβ) V −1
ZiZiV −1
(y − Xβ) = 0, (8)
pro i = 1, . . . , s + 1.
REML odhady
Substitucí y = A y, X = A X = 0, Z = A Z a V = A V A dostaneme restringované
maximálně věrohodné rovnice
Tr[(A V A)−1
A ZiZiA] = y A(A V A)−1
A ZiZiA(A V A)−1
A y, (9)
pro i = 1, . . . , s + 1.
Příklad
Uvažujme nejjednodušší případ y ∼ N(Xβ, V = σ2
I). Tedy i = 1 a Zi = I a podívejme
se, jak se liší ML a REML odhady σ2
.
3
REML Pavla Dokoupilová
ML odhady
Dosazením do rovnice (8) a použitím ML odhadu ˆβ = (X V −1
X)−1
X V −1
y obdržíme
Tr
1
σ2
In×n = (y − X ˆβ)
1
σ2
I
1
σ2
I(y − X ˆβ),
n
σ2
=
(y − X ˆβ) (y − X ˆβ)
(σ2)2
ˆσ2
=
(y − X ˆβ) (y − X ˆβ)
n
.
REML odhady
Dosazením do rovnice (9) a využitím rovnosti (3) dostaneme
Tr[
1
σ2
(A A)−1
A A] = y A
1
σ2
(A A)−1
A
1
σ2
A(A A)−1
A y,
Tr[
1
σ2
I(n−p)×(n−p)] =
y A(A A)−1
A y
(σ2)2
,
n − p
σ2
=
(y − X ˆβ) (y − X ˆβ)
(σ2)2
,
ˆσ2
=
(y − X ˆβ) (y − X ˆβ)
n − p
.
Literatura
[1] LaMotte, L. R. A Direct Derivation of the REML Likelihood Function. Statistical
Papers, 48, 2007, s. 321–327.
[2] Wimmer, G. Bioštatistika.
4