2.3 Gemischte Lineare Modelle
2.3.1 Einf¨ uhrung
Modelle, bei denen sowohl feste als auch zuf¨allige Effekte betrachtet werden, nennt man Gemischte Modelle. Die nachfolgende Darstellung soll zur generellen Ein-f¨uhrung der Bedeutung und der Notation bei gemischten Modellen dienen. Die Theo-rie der gemischten Modelle (
”Mixed Models Theory“) geht zur¨uck auf die nachste-hend genannten Autoren. Grundlagen der Mixed-Models Theorie sind beispielsweise bei Searle et al. [50] ausf¨uhrlich beschrieben.
Das gemischte Modell schreibt sich in Anlehnung an das allgemeine lineare Modell zu
y=Xβ+Zγ mit einer oder mehreren zuf¨alligen Komponenten inγ.
E Ã γ
e
!
= Ã 0
0
!
D=COV Ã
γ e
!
= Ã
G 0
0 R
! .
Xbezeichnet wieder die Designmatix f¨ur die festen Effekteβ, undZ bezeichnet die Designmatrix der zuf¨alligen Effekte γ. Da die block-diagonale Matrix D – wie hier dargestellt – lediglich aus Diagonalelementen besteht, werden also die Faktoren mit zuf¨alligen Effekten stets unabh¨angig vom Restfehleremodelliert; jedoch k¨onnen in-nerhalb der beiden Gruppen theoretisch beliebige Kovarianzstrukturen angenommen werden. Die Varianz-Kovarianzmatrix der Beobachtungen y =yij = (y11, . . . , ypnp)0 ergibt sich somit aus den Elementen vonG und R zu
V ar(yij) = V =ZGZ0+R . (2.1) Der Vektor der Residuen ergibt sich aus
r=y−X(X0V−1X)−X0V−1y . (2.2) W¨ahrend f¨ur den Parameterraum der festen Effekte β im gemischten Modell β∈Rk gilt, k¨onnen die Parameter in γ nur solche Werte annehmen, f¨ur die G und R positiv (semi)definit sind.
Sch¨atzverfahren in gemischten Modellen
Im Vergleich zum Modell mit ausschließlich festen Effekten ist die Sch¨atzung im gemischten Modell aufw¨andiger, da hier eine gr¨oßere Anzahl unbekannter Parameter (simultan) zu sch¨atzen ist, n¨amlich nicht nur σe2 bzw. die Elemente in R, sondern auch die Elemente inG. Kleinste-Quadrate-Methoden sind hier nicht mehr ad¨aquat (vgl. z.B. Laird and Ware [33]).
Zur Sch¨atzung der Variationskomponenten in G und R wird unter Benutzung der Normalverteilung-Annahme von γ und e zumeist mittels Likelihood-basierten Sch¨atzverfahren gearbeitet. Diese Methoden basieren grunds¨atzlich auf dem An-satz, diejenigen Werte f¨ur die zu sch¨atzenden Parameter zu finden, bei denen die Wahrscheinlichkeit (”likelihood“), dass eben diese Parameter-Werte als Realisatio-nen beobachtet werden, maximal wird.
Die allgemeine Form einer Likelihood-Funktion f¨ur einen Parameter(vektor) θ lautet bei gegebenem Beobachtungsvektor y und Dichtefunktion fy:
L(θ) = YN
i,j=1
fy(yij|θ) . (2.3)
Bei der Varianz- und bei der Parametersch¨atzung in gemischten Modellen haben sich heute im Wesentlichen zwei Methoden etabliert:
• Maximum-Likelihood (ML)- und
• Restricted Maximum-Likelihood (REML)-Sch¨atzung.
Allgemein gesprochen unterscheiden sich die beiden Methoden darin, dass bei der ML-Methodealle Parameter im Modell simultan gesch¨atzt werden. Bei der REML-Methode hingegen werden nur die zuf¨alligen Effekte separat, d.h. ohne Ber¨ucksichti-gung der festen Effekte, gesch¨atzt. Die festen Effekte werden dann in einem zweiten Schritt aus den ML-Sch¨atzern bestimmt.
Patterson und Thompson [43] haben gezeigt, dass unter der Nichtbeachtung von β keine Information ¨uber γ verloren geht. Daher erf¨ullen beide Sch¨atzmethoden die w¨unschenswerten Eigenschaften der Konsistenz, Effizienz und der asymptotischen Normalit¨at. Der Vorteil der REML-Methode liegt jedoch darin, dass im Gegensatz zur nicht restringierten ML-Methode die Sch¨atzer der Varianzkomponenten (auf-grund der benutzten Freiheitsgrade) nicht negativ verzerrt sind. Andererseits liefern die ML-Sch¨atzungen der Restvarianz σe2 f¨ur rg(X) < 5 einen geringeren mittleren quadratischen Fehler (MSE) als die REML-Sch¨atzer.
Zur Modellsch¨atzung existiert dar¨uber hinaus eine nicht-iterative Sch¨atzmethode
”MIVQUE0“ bzw.
”MINQUE“ (Minimum Variance Quadratic Unbiased Estimati-on), die von Rao et al. [45] bzw. von Hartley et al. [20] vorgeschlagen wurde. Diese Methode liefert erwartungstreue Sch¨atzer der Varianzkomponenten, die jedoch von den Sch¨atzern der festen Effekte unabh¨angig sind. Wie von Swallow et al. [53] mittels Simulationen gezeigt wurde, besitzt die MINQUE-Methode zwar etwas schlechtere Eigenschaften als die Likelihood-basierten Methoden; auf sie kann aber im Falle von Konvergenz-Problemen bei ML und REML zur¨uckgegriffen werden.
Varianzsch¨atzung im Gemischten Linearen Modell
Die Sch¨atzung der Varianzkomponenten im Modell erfolgt in einem ersten Schritt
¨uber die Maximierung einer Log-Likelihood-Funktion f¨ur Gund R
lML(G, R) =− 1
2ln|V| −1
2r0V−1r− n
2ln(2π) , lREML(G, R) =− 1
2ln|V| −1
2|X0V−1X| − 1
2r0V−1r−N −p
2 ln(2π) .
Die L¨osungen ˆGund ˆR, f¨ur die die Funktionen lML bzw. lREML maximal werden, heißen Maximum-Likelihood-Sch¨atzung, wobei p den Rang von X, also die Zahl der festen Faktor-Stufen-Kombinationen, bezeichnet.
Nur in sehr einfachen Modellen (d.h. Modellen mit nur einem festen Effekt und nur dem Restfehler als zuf¨alligem Effekt) k¨onnen f¨ur die Optimierungen von lML bzw.
lREML geschlossene L¨osungen bestimmt werden; in allen anderen F¨allen m¨ussen nu-merische Algorithmen verwendet werden, die aber zur heutigen Zeit von den g¨angi-gen Rechnern leicht bew¨altigt werden k¨onnen. Weite Verbreitung findet heute der (stabilisierte) Newton-Raphson-Algorithmus (siehe z.B. bei Littell et al. [35]).
Sch¨atzung der Effekte im Gemischten Linearen Modell
Zur Parametersch¨atzung im gemischten Modell wird das Gleichungssystem Ã
X0Rˆ−1X X0Rˆ−1Z Z0Rˆ−1X Z0Rˆ−1Z + ˆG−1
! Ã βˆ ˆ γ
!
= Ã
X0Rˆ−1y Z0Rˆ−1y
!
unter Verwendung der wie oben bestimmten Sch¨atzer f¨ur die Variationskomponen-ten nach ˆβ und ˆγ aufgel¨ost (vgl. Henderson [24]). Die L¨osungen k¨onnen wie folgt dargestellt werden (vgl. Laird und Ware Henderson [33]):
βˆ= (X0Vˆ−1X)−X0Vˆ−1y , ˆ
γ = ˆGZ0Vˆ−1(y−Xβ) .ˆ
A− bezeichnet die generalisierte Inverse einer Matrix A. Wenn γ als feste Effekte gesch¨atzt w¨aren, so stelltG−1 eine Schrumpfung (
”shrinkage“) der Sch¨atzer inγ von beobachteten Mittelwerten hin zum gesch¨atzten Populationsmittelwert ˆµ dar.
Wenn die Variationskomponenten in G bekannt sind (bzw. im Spezialfall, dass keine zuf¨alligen Effekte modelliert sind), heißt der Vektor ˆβ
”best linear unbiased estimator“ (BLUE) vonβ. Der Vektor ˆγ heißt dann
”best linear unbiased predictor“
(BLUP) vonγ.G undR sind in der Praxis zumeist unbekannt und werden nach ei-ner der oben beschriebenen Methoden (wie ML, REML oder MIVQUE0/MINQUE) gesch¨atzt. Die Sch¨atzer ˆβ und ˆγ heißen dann EBLUE und EBLUP, wobei
”E“ f¨ur empirischsteht, da G aus den Daten gesch¨atzt wird.
Die Kovarianzmatrix der Sch¨atzfehler ( ˆβ−β,γˆ−γ) lautet C =
Ã
X0R−1X X0R−1Z Z0R−1X Z0R−1Z+G−1
!−
, (2.4)
deren KomponentenG und R dann ebenfalls durch ihre Sch¨atzer ersetzt werden.
Sch¨atzprobleme
Bei komplexen gemischten Modellen k¨onnen die weiter oben zitierten Iterationsal-gorithmen Konvergenzprobleme in der Weise zeigen, dass die AlIterationsal-gorithmen gar nicht bzw. zu Werten konvergieren, die außerhalb des Parameterraums liegen. In solchen F¨allen kann h¨aufig durch Wahl von geeigneten Startwerten Abhilfe geschaffen wer-den.
Eine weitere Quelle von Sch¨atzproblemen kann bedingt durch geringe Streuungen zwischen den Stufen der zuf¨alligen Effekte gegeben sein, insbesondere wenn negative Werte gesch¨atzt werden (siehe auch Seite 67). Bei Einrichtungsvergleichen kann – da Varianzen nicht negativ sein d¨urfen und so der Sch¨atzfehler minimiert werden kann – hier schlicht 0 als wahre Varianz zwischen den Stufen angenommen werden.
Die Nullhypothese, dass keine Unterschiede zwischen den Stufen (hier Einrichtun-gen) bestehen, kann somit angenommen werden. In anderen Situationen kann dieser Ausweg jedoch nicht angezeigt sein, was eine Reparametrisierung des Modells erfor-derlich machen kann. Ein Vorschlag hierzu ist bei Verbeke und Molenberghs [57], Kapitel 5.6.1, nachzulesen.
Inferenz im gemischten Modell
In gemischten Modellen werden die festen und die zuf¨alligen Effekte getrennt hin-sichtlich ihres (signifikanten) Einflusses auf die Zielgr¨oße getestet.Wald-Z-Statistiken, die auf der Normalverteilungs-Annahme der Effekte beruhen, berechnen sich aus dem Parametersch¨atzer, dividiert durch den asymptotischen Standardfehler SEcβi. Dieser wird aus der Inversen der zweiten Ableitung der Likelihood-Funktion, nach dem jeweiligen Parameterβi, gebildet.
F¨ur ein einzelnes Element βi aus dem gesamten Parametervektor β ist f¨ur die Nullhypothese
H0 :βi =βi0 vs. βi 6=βi0 die Statistik
Z = βˆi−βi0 SEˆ βi
asymptotisch standardnormalverteilt. F¨ur beliebige LinearkombinationenLmit Rang l gilt unter Nullhypothese
H0 :Lβ =β0 vs. Lβ 6=β0 :
Z = ( ˆβ−β0)0L
L Ã n
X
i=1
Xi0Vi−1(ˆγ)Xi
!−1 L0
−1
∼χ2l .
Die gesch¨atzten Standardfehler von ˆβ untersch¨atzen im gemischten Modell aller-dings die tats¨achliche Variabilit¨at inβ (siehe Dempster et al. [9]). Diese Verzerrung kann beispielsweise durch die Verwendung von F- bzw. t-Tests ¨uber die Verteilung von β umgangen werden. Die Zahl der Freiheitsgrade muss hier noch z.B. mittels Satterthwaite-Verfahren [48] aus den Daten gesch¨atzt werden. Weitere Methoden, wie nat¨urlich die Likelihood-Ratio-Tests oder Methoden zur Korrektur der Freiheits-grade in Wald-Statistiken, sind verf¨ugbar. Bei hinreichend großen Stichprobenum-f¨angen, wie es zumeist bei Registern im Gesundheitswesen der Fall ist, werden sich die unterschiedlichen Methoden jedoch hinsichtlich der Ergebnisse kaum auswirken.
Bei Einrichtungsvergleichen ist man ohnehin in der Regel weniger an dem Einfluss-grad der festen Effekte, denn an denen der Varianzkomponenten interessiert. Nach der Bestimmung der Sch¨atzer f¨ur die Zentrumsmittelwerte ist bei Einrichtungsver-gleichen die Bestimmung der Konfidenzintervalle ein zentraler Punkt, da diese in der ¨Offentlichkeit am ehesten verstanden werden und gleichzeitig auch Testcharak-ter haben.
F¨ur Tests zur Bestimmung von zuf¨alligen Effekten im gemischten Modell werden Likelihood-basierte Statistiken verwendet. Eine solche ist die Wald-Z-Statistik, die auf der Normalverteilungs-Annahme der Effekte beruht. Z berechnet sich aus dem Parametersch¨atzer, dividiert durch den gesch¨atzten asymptotischen Standardfehler SEcγi, der sich wiederum aus der Inversen der partiellen Ableitung zweiter Ordnung (nach γ) berechnet.
F¨ur nicht normalverteilte Effekte in G k¨onnen auch Likelihood-Ratio-Statistiken verwendet werden (siehe beispielsweise Self und Liang [51]).