stat_gemischtemodelle.Rmd

---
Title: "Gemischte Modelle"
output: 
  html_document:
    includes:
      after_body: footer.html 
---

# Was ist überhaupt ein zufälliger Effekt ?
Modelle mit zufälligen Effekten sind eine neuere Entwicklung ([Fisher, 1919](https://www.cambridge.org/core/journals/earth-and-environmental-science-transactions-of-royal-society-of-edinburgh/article/xvthe-correlation-between-relatives-on-the-supposition-of-mendelian-inheritance/A60675052E0FB78C561F66C670BC75DE)) als Modelle mit festen Effekten ([Legendre, 1805](https://books.google.de/books/about/Nouvelles_m%C3%A9thodes_pour_la_d%C3%A9terminati.html?id=FRcOAAAAQAAJ&redir_esc=y); [Gauss, 1809](https://www.beck-shop.de/gauss-theoria-motus-corporum-coelestium-sectionibus-conicis-solem-ambientium/product/10712705)). Auch beim Erlernen statistischer Methoden folgt das Anwenden zufälliger Effekte - wenn überhaupt - nachdem einfache Regressionen angepasst und ANOVAs durchgeführt wurden. Noch heute kommt es daher vor, dass mit "einem Effekt" in der Regel ein fester Effekt gemeint ist, der zufällige Effekt also etwas besonderes ist und auch expizit als socher betitelt wird.

Wo immer zufällige Effekte erläutert werden, wird zumeist damit begonnen, dass *sie zufällige Stichproben aus einer Population* darstellen oder, dass es sich um Faktoren handelt, *deren Stufen nicht wie bei festen Effekten bestimmt, sondern zufällig aus einer Gesamtheit ausgewählt wurden*.

Ein anderer Punkt, der oft genannt wird, ist die Mindestanzahl an Faktorstufen, die vorhanden sein sollte bevor der Faktor als zufällig ins Modell genommen werden darf. Dabei schwanken die Empfehlungen meist zwischen 5 und 12.

Es soll übrigens gesagt sein, dass der Fehlerterm ($e$) in allen Modellen - also auch in Modellen mit sonst nur festen Effekten - eine Zufallsvariable, also eine zufällige, stochastische Komponente ist. Streng genommen sind also so gut wie alle Modelle *gemischte Modelle*, doch der Fehler nimmt hier eine besondere Rolle ein und wird somit nicht als zufälliger Effekt gezählt. 

### Faustregeln: Ein Faktor ist (eher) zufällig wenn

 * seine Stufen zufällig aus einer größeren Population gezogen wurden. *Beispiel: Ein paar Orte aus einer Zielregion*
 
 * er mindestens ~8 Stufen aufweist. *Beispiel: Blöcke$^1$ in einem Versuch.*
   + [Beispielkapitel Alpha Design](1f_alpha.html)
   + [Beispielkapitel Augmented Design](1f_augmented_blockfixorrandom.html)
 
 * er eine zusätzliche Randomisationseinheit in (dem Design) eines Versuchs darstellt. *Beispiel: Mainplots/Großteilstücke in Split-Plot-Designs/Spaltanlagen. Sonst ist meist nur die Parzelle Randomisationseinheit und diese wird durch den Fehler ja auch quasi als zufällig abgebildet.*
   + [Beispielkapitel Split-Plot-Design](2f_splitplot.html)
 
 * er unvollständige Blöcke$^1$ in einem Versuch darstellt.
   + [Beispielkapitel Alpha Design](1f_alpha.html)
   + [Beispielkapitel Augmented Design](1f_augmented_blockfixorrandom.html)
 
 * er Mehrfachmessungen darstellt. *Beispiel: Mehrere Pflanzen in derselben Parzelle wurden beprobt. Es gibt also Noise nicht nur aufgrund der Parzellen, sondern auch aufgrund der Pflanzen innerhalb der Parzellen.*
   + [Beispielkapitel Sub-sampling](1f_subsampling.html)
 
 * er mit einem anderen zufälligen Effekt gekreuzt ist. *Beispiel: Der Effekt Ort sei definitiv zufällig. Dann ist der Effekt Jahr:Ort auch zufällig.*
   + [Beispielkapitel Versuchsserie](3f_met_regions.html)
 
 * es nötig ist Korrelationen/Varianzstrukturen zwischen seinen Stufen anzunehmen:
   + zwischen Genotypen aufgrund von genetischer Ähnlichkeit
   + zwischen Messwiederholungen am selben Objekt aufgrund von zeitlicher Nähe
     + [Beispielkapitel Messwiederholung](1f_rcbd_messwdh.html)
   + zwischen Proben/Parzellen aufgrund von räumlicher Nähe
   
>  $^1$**Blöcke als feste oder zufällige Effekte?**: <br>Bei balancierten Daten und ausschließlich vollständigen Blöcken macht es keinen Unterschied ob die Blockeffekte als fest oder zufällig ins Modell genommen werden. Nur bei unbalancierten Daten bzw. unvollständigen Blöcken **kann** es sinnvoll sein die Blockeffekte als zufällig ins Modell zu nehmen, da nur so sowohl *Inter-Block* als auch *Intra-Block Informationen* in der Analyse berücksichtigt werden können. Um zu prüfen ob Blockeffekte als fest oder zufällig ins Modell genommen werden sollten, sollten beide Modelle aufgestellt werden und die s.e.d. (standard error of a difference) der Vergleiche zwischen den Behandlungsstufen verglichen werden [siehe Kapitel zum Augmented Design](1f_augmented_blockfixorrandom.html). Das Modell mit den durchschnittlich kleineren  s.e.d. ist vorzuziehen. Es muss dabei aber unbedingt eine Approximierung der Freiheitsgrade stattfinden, z.B. die Methode von [Kenward & Roger (1997)](https://www.ncbi.nlm.nih.gov/pubmed/9333350) oder von Satterthwaite. Dies wird von der `emmeans()` Funktion allerdings [standardmäßig getan](https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html) (im Gegensatz zu `PROC MIXED` in SAS, wo die [Option](http://support.sas.com/documentation/cdl/en/statug/65328/HTML/default/viewer.htm#statug_glimmix_details39.htm) `/DDFM=KR` angegeben werden muss). Ohne Approximation wird das zufällige Modell stets - fälschlicherweise - die kleineren s.e.d. aufweisen.

# Die generelle Schreibweise
Oft werden gemischte Modelle generell ausgedrückt als
$$ y = X\beta + Zu + e $$

## Was man direkt sieht
In der obigen Formel ist $y$ ein Vektor mit den Beobachtungen (also z.B. die Spalte mit den Ertragswerten). $X$ ist die Designmatrix der festen Effekte und $\beta$ ist der Vektor mit allen festen Effekten. Analog ist $Z$ die Designmatrix der zufälligen Effekte und $u$ der Vektor mit den zufälligen Effekten. Deshalb kann ein einfaches lineares Modell ohne zufällige Effekte auch als $y = X\beta + e$ ausgedrückt werden. 

Bevor ein Modell angepasst wird, sind $\beta$ und $u$ natürlich noch unbekannt in dem Sinne, dass wir noch keine Zahlen für die Effekte haben - sie müssen also noch geschätzt werden. $y$, $X$ und $Z$ sind aber sehr wohl bekannt - schließlich stecken in ihnen der Messwert selbst, aber einfach ausgedrückt auch welche Behandlungen o.ä. dieser Messwert erfahren hat. Man kann auch sagen, dass die Designmatrizen $X$ und $Z$ die Effekte mit den Beobachtungswerten zusammenbringen.

Schließlich haben wir noch den Zufallsvektor der Fehler $e$. Es muss klar sein, dass die Länge von $y$ und $e$, sowie die Zeilenanzahl von $X$ und $Z$ der Anzahl Beobachtungen im Datensatz entsprechen. Die Länge von $\beta$ hängt stattdessen davon ab wie viele feste Effekte im Modell angepasst werden sollen (also mindestens $\mu$) und entspricht der Spaltenanzahl von $X$.  Für $u$ gilt das analog für die zufälligen Effekten und $Z$.

## Was man nicht direkt sieht
Allein mit den Infos in den Absätzen hierüber hat man ein gemischtes Modell noch nicht vollständig beschrieben. In der Regel sollte man nämlich auch die Verteilungen und Kovarianzstrukturen erwähnen:

$$ u \sim MVN(0,G)  $$
$$ e \sim MVN(0,R)  $$
$$ y \sim MVN(X\beta,V)  $$

### Verteilung
Wenn nichts anderes gesagt ist, wird wie so oft von einer multivariaten Normalverteilung ausgegangen - man kann sie aber auch explizit erwähnen so wie hier durch $MVN()$. Würde es sich hier nicht um die Normalverteilung handeln, so sprächen wir von generalisierten gemischten Modellen.

### Erwartungswert
Innerhalb der Klammer steht dann der erste Eintrag für den Erwartungswert (*mean vector*). Man sieht, dass sowohl die zufälligen Effekte $u$, als auch die Fehler $e$ einen Erwartungswert von 0 haben. Das ergibt in sofern Sinn, als dass beide ja ledigliche *zufällige Abweichungen* von dem festen Teil des Modells ausdrücken sollen. Demnach sind einzelne Fehler mal positiv und mal negativ aber im Schnitt ist ihr Erwartungswert eben 0. Der Erwartungswert für $y$, also für unsere Beobachtungen ist natürlich nicht 0, sondern hängt von den Daten ab bzw. wird von den festen Effekten ausdgedrückt, die anhand der Daten geschätzt werden. 

Als Beispiel kann man sich einfach ein Modell mit $\mu$ als einzigem festen Effekt vorstellen: Dann gilt $\beta=\mu$ und $X$ ist nur ein Vektor gefüllt mit $1$en dessen Länge der Anzahl Beobachtungen entspricht. Einfach ausgedrückt würde in einem so simplen Fall der gesamt feste Teil des Modells auf einen Gesamtmittelwert $\mu$ runtergebrochen werden und der Erwartungswert für unsere Beobachtungen $y$ wäre dann eben dieser Gesamtmittelwert. 

### Kovarianzmatrix
Vielleicht der interessanteste Teil an gemischten Modellen *versteckt* sich allerdings in den Kovarianzmatrizen $G$, $R$ und $V$. Diese Kovarianzmatrizen (manchmal auch einfach Varianzmatrizen genannt) enthalten, egal ob bezogen auf zufällige Effekte($u$), Fehler($e$) oder Beobachtungen($y$), jeweils Informationen zu zwei Dingen: Den Varianzen und den Kovarianzen/Korrelationen. Dabei stehen die Varianzen auf der Diagonale und die Kovarianzen/Korrelationen daneben (*off-diagonal*). Deshalb sind Kovarianzmatrizen auch immer quadratisch und ihre Breite und Länge genau so groß wie die Länge ihres zugehörigen Vektors ($u$/$e$/$y$).

Nehmen wir einen Vektor mit den Zufallsvariablen $a$, $b$ und $c$ als Beispiel und betrachten dessen Kovarianzmatrix, die wir $M$ nennen:

$$ M_{i,j} = Cov\begin{pmatrix} a \\\ b \\\ c \end{pmatrix} = \begin{pmatrix} Var(a) & Cov(a,b) & Cov(a,c) \\\ Cov(b,a) & Var(b) & Cov(b,c) \\\ Cov(c,a) & Cov(c,b) & Var(c) \end{pmatrix} = \begin{pmatrix} \sigma^2_a & \sigma_{a,b} & \sigma_{a,c} \\\ \sigma_{b,a} & \sigma^2_b & \sigma_{b,c} \\\ \sigma_{c,a} & \sigma_{c,b} & \sigma^2_c \end{pmatrix} $$

Auf der Hauptdiagonale (also wenn $i=j$ für $M_{i,j}$) befinden sich alle Varianzen der einzelnen Elemente $a$, $b$ und $c$. Varianzen sind immer positiv. Auf der *off-diagonal* (also wenn $i\ne j$ für $M_{i,j}$) befinden sich alle Kovarianzen. Die Kovarianzmatrix ist symmetrisch, da die Kovarianz zweier Zufallsvariablen symmetrisch ist ($\sigma_{i,j}=\sigma_{j,i}$). Kovarianzen können auch negativ sein. 

> **Kovarianz(matrizen) und Korrelation(smatrizen)**
> <br>
> Eine Kovarianz kann auch als Korrelation ausgedrückt werden und umgekehrt - man muss sie lediglich umrechnen:
$$ r(x,y) = \frac{Cov(x,y)}{\sqrt{Var(x)Var(y)}}$$
Demnach ist die Korrelation zwischen $x$ und $y$ gleich deren Kovarianz geteilt durch die Wurzel aus dem Produkt ihrer Varianzen. Kovarianzen bewegen sich dabei auf derselben Skala wie Varianzen, während Korrelationen ja standardisiert sind, also nur Werte zwischen -1 und 1 annehmen können. Ich persönlich bevorzuge Korrelationen, da ich sie intuitiver interpretieren kann.
> <br>
> Auch Kovarianzmatrizen können als Korrelationsmatrizen geschrieben werden: 
> <br> Die Kovarianzmatrix $\begin{pmatrix}220 & 111\\\ 111 & 241\end{pmatrix}$ entspricht der Korrelationsmatrix $\begin{pmatrix}220 & 0,48\\\ 0,48 & 241\end{pmatrix}$, weil  $\frac{111}{\sqrt{220 \times 241}} = 0,48$.

Bei der Anpassung eines gemischten Modells wird deshalb oft davon gesprochen $G$ und/oder $R$ eine bestimmte Struktur vorzugeben. Tatsächlich nehmen die beiden nämlich standardmäßig - also wenn nichts anderes angegeben wird - die einfachste Kovarianzstruktur an: Die Einheitsmatrix (*identity*), welche auch als $I_n$ geschrieben werden kann, wobei $n$ für die Anzahl Zeilen/Spalten steht:

$$ I_3 = \begin{pmatrix} 1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & 1 \end{pmatrix} $$

Demnach haben alle der jeweiligen Zufallsvariablen dieselbe Varianz und keine Kovarianz/Korrelation. Zur Erinnerung: Wenn etwas keine Korrelation hat, bzw. eine Korrelation=0, dann spricht man davon, dass es unabhängig ist. Außerdem muss klar sein, dass $I_3$ nur die Struktur angibt und die Varianzen natürlich nicht genau $1$ sein müssen, sondern $I_3$ wird noch mit einem $\sigma^2$ multipliziert, welches dann als Varianzkomponente geschätzt wird. Aber es gilt eben, dass eine gemeinsame Varianzkomponente für alle drei Elemente geschätzt wird. Mehr zu weiteren möglichen Varianzstrukturen findet sich z.B. in [diesem](appendix_kovarstrukt.html) oder [diesem Kapitel](3f_met_regions.html).

Zuletzt noch eine Verdeutlichung: Wenn also keine besondere Varianzstruktur angepasst wurde und ein Datensatz mit 50 Beobachtungen ausgewertet wurde, dann ist die Kovarianzmatrix der Fehler $e$ gleich $Cov(e)=R=I_{50}\sigma^2_e$, wobei der Schätzwert für $\sigma^2_e$ im Output der Analyse als Fehlervarianz angegeben wird.

# Weitere Informationen
## ML und REML
Im Gegensatz zu einfachen linearen Modellen, können weder generalisierte, noch gemischte Modelle nur mit der **Kleinst-Quadrat-Methode (KQ)** ([Legendre, 1805](https://books.google.de/books/about/Nouvelles_m%C3%A9thodes_pour_la_d%C3%A9terminati.html?id=FRcOAAAAQAAJ&redir_esc=y)) auskommen, d.h. all ihre Parameter schätzen.

Stattdessen muss z.B. die **Maximum-Likelihood-Methode (ML)** ([Edgeworth, 1908](https://www.jstor.org/stable/2339293?origin=crossref&seq=1#page_scan_tab_contents)) genutzt werden. Wie der Name schon sagt, schätzt diese Methode den Parameter mithilfe einer Wahrscheinlichkeitsfunktion so, dass dessen Verteilung am besten zur Realisierung der beobachteten Daten passt. Die Schätzwerte sind also *maximum likely* - am wahrscheinlichsten. Man muss aber durchaus selbst angeben für *welche* Verteilung das durchgeführt werden soll. Wenn man übrigens ein einfaches lineares Modell (z.B. eine lineare Regression unter Annahme der Normalverteiliung) mit der ML-Methode schätzt, sind die Schätzwerte der Parameter ($\beta$) identisch zu denen der KQ-Methode. 

Ein Vorteil der ML-Methode ist, dass sie auch für Modelle eingesetzt werden kann, die eben nicht annehmen, dass der Fehlerterm normalverteilt annehmen. Stattdessen können auch andere Verteilungen aus der Klasse der *Exponentialfamilien* angenommen werden - dazu gehören neben der Normalverteilung auch die Binomial-, Poisson-, Gamma- und inverse Gaußverteilung.

Ein Nachteil der ML-Methode ist, dass sie in gemischten Modellen (selbst für Normalverteilung) die wahren Varianzkomponenten unterschätzt und das besonders bei kleinen Stichproben. 

Eine unverzerrte Schätzung der Varianzkomponenten liefert hingegen die **Restricted-Maximum-Likelihood-Methode (REML)** ([Bartlett, 1937](https://royalsocietypublishing.org/doi/10.1098/rspa.1937.0109)). Es wird aber auch wirklich nur $V$ und nicht $\beta$ mit mit der REML-Methode geschätzt. Stattdessen wird in gemischten Modellen die REML-Schätzung von $V$ genutzt um mit ihr eine ML-Schätzung von $\beta$ durchzuführen. Das bedeutet demnach auch, dass die mit REML geschätzten Varianzkomponenten unabhängig der festen Effekte im Modell sind.

Ein Nachteil von REML gegenüber ML ist, dass die Modellanpassungen von verschiedenen Modellen für denselben Datensatz nur bedingt miteinander vergleichbar sind: Man darf nämlich nur solange die Anpassungsgüten verschiedener REML Modelle vergleichen, wie sie die gleichen festen Effekte aufweisen. Andersherum ausgedrückt: Welches Modell besser auf die Daten angepasst werden konnte, darf nicht zwischen Modellen verglichen werden, die mit REML geschätzt wurden und unterschiedliche feste Effekte haben.

## BLUP und BLUE
Während die REML-Methode genutzt wird um die Varianzkomponenten zu schätzen, werden die zufälligen Effekte selbst als **Best linear unbiased prediction (BLUP)** *predicted/vorhergesagt*. Die festen Effekte hingegen werden als **Best linear unbiased estimation (BLUE)** *estimated/geschätzt*.