\[ \def\logit{\text{logit}} \def\Var{\text{Var}} \def\iid{\text{i.i.d.}} \newcommand{\vect}[1] {\left( \begin{array}{c} #1 \end{array}\right)} \def\Cor{\text{Cor}} \def\Cov{\text{Cov}} \]
lme4 und lmerTest.lmer(Y~x+(1|ID)).Abschnitte mit * sind für die Interessierten.
Das Allgemeine Lineare Modell mit kontinuierlicher abhängiger Variable \(Y\) war:
\[ \boxed{Y_i=\mathbf{x}_i^ T\mathbf{\beta}+\epsilon_i}, \quad i=1,\dots, n \]
Dabei stand \(\mathbf{x}_i^ T\mathbf{\beta}\) für den linearen Prädiktor (\(T\) steht für transponiert; wenn \(\mathbf{x}\) ein Kolonnenvektor ist, ist \(\mathbf{x}^T\) ein Zeilenvektor). Der lineare Prädiktor ist das Skalarprodukt
\[ \mathbf{x}_i^T\mathbf{\beta}=(1, x_{i1}, \ldots, x_{im}) \vect{\beta_0\\\vdots \\ \beta_{im}}=\beta_0+\beta_1x_{i1}+\dots+\beta_mx_{im} \]
Dabei nahmen wir unabhängige und identisch verteilte stochastische Fehler an, \(\epsilon_i\, \iid\), meistens \(\epsilon_i \sim \mathcal{N}(0,\sigma^2)\). Mit \(m\) Eingangsgrössen und einem Intercept, also mit \(p=m+1\) Parametern und \(n\) Beobachtungen können wir das Modell auch in Matrixform schreiben:
\[\begin{equation} \label{eq:42} \boxed{\overset{n\times 1}{\mathbf{Y}}=\overset{n\times p}{X}\times\overset{p\times 1}{\mathbf{\beta}}+\overset{n\times 1}{\mathbf{\epsilon}}} \end{equation}\]
\[\begin{equation} \small \label{eq:43} \vect{Y_1\\Y_2\\\vdots\\Y_i\\\vdots\\Y_n}= \begin{pmatrix} 1& x_{11}&\ldots&x_{1m}\\ 1&x_{21}&\ldots &x_{2m}\\ \vdots&\vdots&\vdots\\ 1&x_{i1}&\ldots& x_{im}\\ \vdots&\vdots&\vdots\\ 1&x_{n1}&\ldots &x_{nm} \end{pmatrix} \vect{\beta_0\\ \vdots \\\beta_m} +\vect{\epsilon_1\\\epsilon_2\\\vdots\\\epsilon_i\\\vdots\\\epsilon_n}. \end{equation}\]
Man nennt die Matrix \(X\) Design-Matrix. Wir haben also \(n\) Gleichungen mit \(p\) Unbekannten.
Die Methode der kleinsten Quadrate haben wir kennengelernt in Section 9.7.1. Die Optimierungsfunktion ist also für dieses Modell
\[ \mathbf{\hat{\beta}}=\arg \underset{\mathbf{\beta}}{\operatorname{min\,}}||\mathbf{Y}-X\mathbf{\beta}||^2. \]
Ableiten der Optimierungsfunktion nach \(\mathbf{\hat{\beta}}\) gibt \[\begin{equation} (-2)X^T(\mathbf{Y}-X\mathbf{\hat{\beta}})=\mathbf{0} \end{equation}\]
Dies führt uns zu den Normalgleichungen \[\begin{equation} \overset{p\times p}{X^TX}\quad\overset{p\times 1}{\mathbf{\hat{\beta}}}=\overset{p\times 1}{X^T\mathbf{Y}}. \end{equation}\]
\(X^T\) ist die transponierte Design-Matrix. Das sind dann \(p\) lineare Gleichungen für \(p\) Unbekannte (Komponenten von \(\mathbf{\hat{\beta}}\)). Wenn die Kolonnen der Design-Matrix \(X\) linear unabhängig sind, dann ist die \(p\times p\)-Matrix \(X^TX\) invertierbar und der Kleinste-Quadrate-Schätzer ist dann \[\begin{equation} \boxed{\mathbf{\hat{\beta}}=(X^TX)^{-1}X^T\mathbf{Y}} \end{equation}\] und kann aus den Daten berechnet werden. Aus den Residuen \(r_i=Y_i-\mathbf{x}_i^T\mathbf{\hat{\beta}}\) bekommt man schliesslich eine Schätzung für \(\sigma^2\), \[\begin{equation} \label{eq:s2hat} \boxed{\hat{\sigma}^2=\frac{1}{n-p}\sum_{i=1}^n r^2_i} \end{equation}\] Dabei entspricht \(n-p\), also Anzahl Beobachteten minus die Anzahl Parameter dem Freiheitsgraden.
Für die Interessierten:
Wenn \(\epsilon_i \,\iid \sim \mathcal{N}(0,\sigma^2)\). Dann gilt
mit
\[\begin{align} \Cov(\mathbf{\hat\beta})= \begin{pmatrix} \Var(\hat{\beta}_1) & \Cov(\hat{\beta}_1,\hat{\beta}_2) & \cdots & \Cov(\hat{\beta}_1,\hat{\beta}_p) \\ \\ \Cov(\hat{\beta}_2,\hat{\beta}_1) & \Var\hat{\beta}_2) & \cdots & \Cov(\hat{\beta}_2,\hat{\beta}_p) \\ \\ \vdots & \vdots & \ddots & \vdots \\ \\ \Cov(\hat{\beta}_p,\hat{\beta}_1) & \Cov(\hat{\beta}_p,\hat{\beta}_2) & \cdots & \Var(\hat{\beta}_p) \end{pmatrix} \end{align}\]
Die Wurzeln der Diagonalelemente entsprechen dann den Standardfehlern. \(\Cov(\mathbf{\hat\beta})\) kann man aus den Daten berechnen mit \(\sigma^2 (X^TX)^{-1}\).
Die Normalverteilung der Fehler ist oft nicht erfüllt. Für grosse \(n\) sind die Aussagen dann trotzdem annähernd wahr (Zentraler Grenzwertsatz). Das ist dann die herkömmliche Begründung in der Praxis für die Konstruktion von Konfidenzintervallen und Tests für die Parameter des Modells.
Wir kommen nun zu einer Erweiterung des Linearen Modells, zum Linear Mixed Model (LMM) (Lineares Gemischtes Modell). Oft haben wir in der Gesundheitsforschung zum Beispiel Messwiederholungen innerhalb einer statistischen Einheit (mehrere Beobachtungen innerhalb einer Person) oder wir haben mehrere Beobachtungen innerhalb eines Clusters (z.B. einem Spital). Die Beobachtungen innerhalb einer Person (oder eines Clusters) sind dann ähnlicher als die Beobachtungen zwischen Personen (oder Clustern). Die Beobachtungen innerhalb einer statistischen Einheit sind dann nicht mehr unabhängig. Wir haben dann korrelierte Fehler, die wir berücksichtigen müssen. Die statistische Einheit ist die Person oder das Cluster, nicht die Beobachtung.
Bedingte Unabhängigkeit spielt eine zentrale Rolle bei den Modellen, die wir unten einführen.
Bei Jugendlichen sind Körpergrösse und Wortschatz offensichtlich nicht unabhängig. Aber sie sind bedingt unabhängig, wenn wir für jeden Jugendlichen eine latente Variable wie Maturität kennen würden. Wenn man den Reifegrad (Maturität) kennt, liefert die Körpergrösse keine zusätzliche Information über den Wortschatz. Das gilt natürlich unter der Annahme, dass die Körpergrösse nur über die Maturität mit dem Wortschatz zusammenhängt.
Dieses Prinzip wird nun in Modellen mit latenten Variablen angewandt. Man nimmt an, dass die manifesten Variablen nichts mehr gemeinsam haben, nachdem man für eine latente Variable kontrolliert hat.
Wir können dies z.B. anwenden für gepaarte Beobachtungen \(Y_{i1}\) und \(Y_{i2}\) innerhalb einer \(i\)ten Einheit. Beobachtungen innerhalb einer statistischen Einheit \(i\) sind nun korreliert, weil sie einen gemeinsamen Term \(U_i\) enthalten. Dieser gemeinsame Term ist latent. Wenn wir ihn kennen würden (auf ihn bedingen), dann sind die Beobachtungen innerhalb einer Einheit unabhängig. Die bedingte Verteilung von \(Y_1\) gegeben \(Y_2\) und \(U_i\) ist dann gleich der bedingten Verteilung von \(Y_1\) gegeben nur \(U_i\). \[ \boxed{\Pr(Y_{i1}\mid Y_{i2},U_i)=\Pr(Y_{i1} \mid U_i)} \]
Wenn man also den latenten Effekt \(U_i\) kennt, dann liefert die Information über \(Y_2\) keinen zusätzlichen Erkenntnisgewinn über die Verteilung von \(Y_1\). Der gemeinsame Einfluss ist dann herausgerechnet – und die Beobachtungen \(Y_{ij}\) sind nur noch durch unabhängige Fehler beeinflusst.
Wenn wir \(U_i\) nicht kennen, dann sind \(Y_{i1}\) und \(Y_{i2}\) marginal abhängig. \[ \boxed{\Pr(Y_{i1}\mid Y_{i2})\neq \Pr(Y_{i1})} \]
Wir können die Korrelation der Beobachtungen innerhalb einer Einheit (Person, Cluster) durch das Hinzufügen von zufälligen oder latenten Variablen \(\mathbf{U}_i\) berücksichtigen. Diese zusätzliche(n) Zufallsvariable(n) induzieren eine Korrelation der wiederholten Beobachtungen und erklären dann die Korrelation der wiederholten Beobachtungen. Gegeben die latente(n) Variable(n) sind die Beobachtungen unabhängig.
Der Vorteil, die \(\mathbf{U}_i\) als zufällig zu behandeln, ist, dass wir weniger Parameter brauchen. In linearen Mixed Modells (LMM) sind die zusätzlichen Parameter Varianzkomponenten, die die Varianz der \(\mathbf{U}_i\) spezifizieren, nicht die \(\mathbf{U}_i\) selber. Zudem hätten Fixed-Effects-Parameter \(\mathbf{U}_i\) hätten keine Interpretation als Populationsparameter.
Für die Interessierten schauen wir uns zuerst den allgemeinen Fall an. Die für uns wichtigen Spezialfälle weiter unten sind dann einfacher. Neu ist nun der \(\color{red}{rote}\) Teil in der Modellformel.
Das Modell für die Einheit \(i\) (Person, Cluster) mit \(n_i\) Beobachtungen ist
\[ \mathbf{Y}_i=\underbrace{X_i\mathbf{\beta}}_{fixed}+\underbrace{\color{red}{Z_i \mathbf{U}_i}}_{random}+\underbrace{\epsilon_i}_{random}, \quad i=1,\dots,n \]
\(\mathbf{Y}_i\) ist ein \(n_i \times 1\)-Vektor der kontinuierlichen Zielgrössen für Einheit \(i\) (Individuum oder Cluster)
\(X_i\) ist eine \(n_i \times p\) Design Matrix von erklärenden Variablen (Kovariaten). Das sind die fixed effects.
\(\mathbf{\beta}\) ist ein \(p\times 1\) Vektor von Populationsparametern, fixed effects Koeffizienten
\(\color{red}{Z_i}\) ist eine \(n_i \times k\) Design Matrix von Zufallseffekten
\(\color{red}{\mathbf{U}_i}\) ist ein \(k\times 1\) Vektor von Zufallseffekten mit Durchschnitt 0 und Kovarianzmatrix \(D\)
\(\mathbf{\epsilon}_i\) ist ein \(n_i \times 1\) Fehlerterm mit unabhängigen Komponenten, jede mit Durchschnitt 0 und within-subject Varianz \(\sigma^2\)
Bezüglich Erwartungswert und Varianz gilt:
Wenn z.B. \(Z_i\) eine Kolonne aus Einsen (\(k=1\)), haben wir das Random Intercept Model \[ \mathbf{Y}_i=X_i\mathbf{\beta}+\mathbf{U}_i+\mathbf{\epsilon}_i, \quad i,=1,\dots,n \]
mit bedingter Kovarianz \(\Cov(\mathbf{Y}_i\mid \mathbf{U}_i)\)= \[\begin{bmatrix} \sigma^2 & 0 & \cdots & 0 \\ 0 & \sigma^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \sigma^2 \\ \end{bmatrix}\] und marginaler Kovarianz \(\Cov(\mathbf{Y}_i)\)= \[\begin{bmatrix} \nu^2+\sigma^2 & \nu^2 & \cdots & \nu^2 \\ \nu^2 & \nu^2+\sigma^2 & \cdots & \nu^2 \\ \vdots & \vdots & \ddots & \vdots \\ \nu^2 & \nu^2 & \cdots & \nu^2+\sigma^2. \\ \end{bmatrix}\]Dies entspricht einem Modell mit konstanter Korrelation (Intraklassenkorrelation)) \[ \rho=\frac{\nu^2}{\nu^2+\sigma^2}. \]
Es gibt viele Modelle mit nicht-konstanter Korrelation (z.B. auto-regressive Fehlern), die wir hier nicht behandeln.
Wir bleiben nun beim Random Intercept Model, welches das bei weitem häufigste LMM darstellt. Eine Zeile (eine Beobachtung \(j\) von Einheit \(i\)) schreiben wir dann \[ \boxed{Y_{ij}=\mathbf{x}^T_{ij}\mathbf{\beta}+{\color{red}{U_i}}+\epsilon_{ij}}, \quad i=1,\dots,n, \quad j=1,\dots,n_i. \tag{12.1}\]
\(n\) steht für die Anzahl Personen oder Cluster und \(n_i\) für die Anzahl Beobachtungen für Einheit \(i\). Wir haben also total \(\sum_{i=1}^n n_i\) Beobachtungen. Im Gegensatz zum linearen Modell haben wir jetzt zwei Zufallsvariablen, \(U_i\) zusätzlich zu \(\epsilon_{ij}\). Wir nehmen für beide Normalverteilungen an, nämlich \(U_i\sim \mathcal{N}(0, \nu^2)\) und \(\epsilon_{ij}\sim \mathcal{N}(0,\sigma^2)\).
Between-und within-subject Varianz: Der Parameter \(\nu^2\) stellt nun die between-subject Varianz und der Parameter \(\sigma^2\) stellt die within-subject Varianz dar.
Verglichen mit einem linearen Modell haben wir also zusätzlich ein random intercept, \[\beta_{0i}=\beta_0+U_i.\] Wir können 12.1 dann schreiben als \[ Y_{ij}=\underbrace{\beta_0+U_i}_{\beta_{0i}}+\beta_1x_{i1}+\dots+\beta_m x_{im}+\epsilon_{ij}. \]
Erwartungswert und Varianz sind dann:
Das grundlegendste LMM ist ein Modell, das wir schon kennen, nämlich das Modell, das einem gepaarten \(t\)-Test entspricht.
Mit einen dichotomen Prädiktor \(T\), der z.B. einen pre- und eine post-Zeitpunkt definiert, ist das Modell \[ Y_{ij}=\beta_0+\beta_1T_{ij}+U_i+\epsilon_{ij},\qquad i=1,\dots, n, \quad j=1,2. \qquad T_{ij}=\begin{cases} 1 &\mbox{wenn } j=2, \\ 0 & \mbox{wenn } j=1 \end{cases} \]
Wir haben wieder \(U_i\sim \mathcal{N}(0, \nu^2),\quad \epsilon_{ij}\sim \mathcal{N}(0,\sigma^2)\) und \(\Cor(Y_{ij},Y_{ij'})=\frac{\nu^2}{\nu^2+\sigma^2}\)
Die Quantität von Interesse ist wie beim gepaarten \(t\)-Test die erwartete Veränderung zwischen \(t_1\) und \(t_2\). Diese Quantität entspricht dem Parameter \(\beta_1\).
Das data.frame d.long2 beinhaltet eine kontinuierliche abhängige Variable zu zwei Zeitpunkten mit \(n=30\) und \(n_i=2\) für alle \(i\). Daten sind im long format, d.h. wir haben eine Zeile für jede Beobachtung, also zwei Zeilen für jede statistische Einheit.
d.long2<-read.csv("https://raw.githubusercontent.com/mcdr65/StatsRsource/master/Data/Rep2.csv",stringsAsFactors=TRUE)
str(d.long2)'data.frame': 60 obs. of 3 variables:
$ subject : Factor w/ 30 levels "s1","s10","s11",..: 1 1 12 12 23 23 25 25 26 26 ...
$ time : Factor w/ 2 levels "t1","t2": 1 2 1 2 1 2 1 2 1 2 ...
$ response: num 22.9 38.4 10.8 18.2 33.2 ...library(psych)
headTail(d.long2) subject time response
1 s1 t1 22.93
2 s1 t2 38.43
3 s2 t1 10.8
4 s2 t2 18.17
... <NA> <NA> ...
57 s29 t1 6.53
58 s29 t2 16.07
59 s30 t1 34.25
60 s30 t2 42.4Mit reshape() kann man Daten im wide Format in long Format umwandeln und umgekehrt. Das Studieren der Beispiele in ?reshape lohnt sich, da entsprechende Umwandlungen doch oft benötigt werden.
d.wide2<-reshape(d.long2,timevar="time",idvar="subject",direction="wide")
headTail(d.wide2) subject response.t1 response.t2
1 s1 22.93 38.43
3 s2 10.8 18.17
5 s3 33.22 37.9
7 s4 24.82 41.64
... <NA> ... ...
53 s27 29.69 36.21
55 s28 23.93 33.57
57 s29 6.53 16.07
59 s30 34.25 42.4Im wide Format haben wir nun eine Zeile pro statistische Einheit, dafür mehrere (hier zwei) Kolonnen pro Zeitpunkt.
Die beobachtete Korrelation der wiederholten Messungen bekommen wir mit
cor(d.wide2[,-1]) response.t1 response.t2
response.t1 1.0000000 0.7097983
response.t2 0.7097983 1.0000000Um ein Linear Mixed Model anzupassen, brauchen wir Daten immer im Long-Format. Eine Zeile entspricht dann einer Beobachtung, nicht einer statistischen Einheit.
Zuerst beschreiben wir die Daten
psych::describeBy(d.long2$response,group=d.long2$time,mat=TRUE) item group1 vars n mean sd median trimmed mad min
X11 1 t1 1 30 24.73446 9.147931 23.27004 24.72994 9.938112 6.528457
X12 2 t2 1 30 34.89181 9.251691 36.44469 35.72962 6.453890 9.754138
max range skew kurtosis se
X11 41.38815 34.85970 0.07279288 -1.1000752 1.670176
X12 52.97559 43.22145 -0.75938610 0.5683042 1.689120##Alternativen
tapply(d.long2$response,d.long2$time,summary)$t1
Min. 1st Qu. Median Mean 3rd Qu. Max.
6.528 17.648 23.270 24.734 33.667 41.388
$t2
Min. 1st Qu. Median Mean 3rd Qu. Max.
9.754 32.226 36.445 34.892 40.367 52.976 sapply(split(d.long2$response,d.long2$time),summary) t1 t2
Min. 6.528457 9.754138
1st Qu. 17.648453 32.225761
Median 23.270037 36.444692
Mean 24.734463 34.891809
3rd Qu. 33.666633 40.366996
Max. 41.388153 52.975589psych::describe(d.wide2[,-1]) vars n mean sd median trimmed mad min max range skew
response.t1 1 30 24.73 9.15 23.27 24.73 9.94 6.53 41.39 34.86 0.07
response.t2 2 30 34.89 9.25 36.44 35.73 6.45 9.75 52.98 43.22 -0.76
kurtosis se
response.t1 -1.10 1.67
response.t2 0.57 1.69Mit arsenal::paired() bekommen wir
tmp<-arsenal::paired(time~response,id=subject,data=d.long2,test=FALSE)
summary(tmp,text=TRUE)
| | t1 (N=30) | t2 (N=30) | Difference (N=30) |
|:------------|:--------------:|:--------------:|:-----------------:|
|response | | | |
|- Mean (SD) | 24.734 (9.148) | 34.892 (9.252) | 10.157 (7.009) |
|- Range | 6.528 - 41.388 | 9.754 - 52.976 | -3.896 - 27.889 |Ein Spaghetti plot ist bei gepaarten Daten besser geeignet als ein Boxplot, da aus einem Boxplot das Design nicht eindeutig ist. Die Paarung der Beobachtungen ist nicht ersichtlich.
with(d.long2,interaction.plot(time,subject,response,ylab="response"))
boxplot(response~time,d.long2)
Es gibt viele Alternativen zum Plotten, z.B. mit lattice oder ggplot2:
lattice::xyplot(response~time,groups=subject,data=d.long2,type="b")
##
library(ggplot2)
p <- ggplot(data = d.long2, aes(x = time, y = response, group = subject))
p <- p+geom_point()+geom_line()+stat_smooth(aes(group = 1),method="lm",se=FALSE)
p <- p + stat_summary(aes(group=1), geom = "point", fun = mean,shape = 17, size = 4)
p
Ein gepaarter Test ist bekannterweise äquivalent zu einem one-sample Test der Veränderung:
pre<-d.long2$response[d.long2$time=="t1"]
post<-d.long2$response[d.long2$time=="t2"]
change<-post-pre
print(t.test(post,pre,paired=TRUE,data=d.long2),digits=6)
Paired t-test
data: post and pre
t = 7.937, df = 29, p-value = 9.4e-09
alternative hypothesis: true mean difference is not equal to 0
95 percent confidence interval:
7.53997 12.77472
sample estimates:
mean difference
10.1573 print(t.test(change),digits=6)
One Sample t-test
data: change
t = 7.937, df = 29, p-value = 9.4e-09
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
7.53997 12.77472
sample estimates:
mean of x
10.1573 aov() dient wie lm() zum Anpassen linearer Modelle. Der Unterschied ist die Art und Weise, wie print() und summary() den Fit behandeln. Mit aov() ist der Ausdruck in der traditionellen Sprache der Varianzanalyse anstelle der linearen Modelle. Wenn die Modellformel einen einzelnen Fehlerterm enthält, wird dieser verwendet, um Fehlerstrata zu spezifizieren, und es werden geeignete Modelle innerhalb jeder Fehlerschicht angepasst.
m.aov<-aov(response~time+Error(subject),data=d.long2)summary(m.aov)
Error: subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 29 4197 144.7
Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
time 1 1547.6 1547.6 63 9.4e-09 ***
Residuals 29 712.4 24.6
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Wir können die entsprechenden Quadratsummen von aov() reproduzieren. Dazu passen wir folgende Modelle an:
mod0 <- lm(response~1,d.long2)
mods <- lm(response~subject,d.long2)
modt <- lm(response~time,d.long2)
modts <-lm(response~subject+time,d.long2)Die Residual-Quadratsummen der Modelle sind dann:
rss.0 <- sum((mod0$residuals)^2)
rss.s <- sum((mods$residuals)^2)
rss.t <- sum((modt$residuals)^2)
rss.ts<- sum((modts$residuals)^2)Wir haben dann:
rss.0 #SS_total[1] 6456.65rss.0-rss.t #SS_time[1] 1547.575rss.0-rss.s #SS_subject[1] 4196.653rss.ts #SS_residual_within[1] 712.4216rss.0-rss.t+rss.ts#SS_total_within[1] 2259.997Repeated Measures ANOVA ist meistens nicht empfehlenswert für die Analyse von wiederholten Messungen. Eine flexiblere und moderne Alternative ist ein LMM, das wir jetzt an die Daten anpassen.
Lineare Gemischte Modelle (LMM) sind die moderne Alternative zur Analyse von Messwiederholungen oder Beobachtungen mit korrelierten Fehlern. Der Vorteil gegenüber Repeated-Measures ANOVA ist, dass LMM nicht-balancierte Daten oder Daten mit fehlenden Werten besser verarbeiten. Klassische Repeated Measures ANOVA aov(...+Error(subject)) brauchen balancierte Daten ohne Missings.
Wir verwenden im Folgenden die Funktion lmer() aus dem Paket lme4 Zusätzlich nutzen wir auch das Paket lmerTest (Mehr dazu später). LMM werden mit der Maximum-Likelihood-Schätzung (siehe Section 9.7.3) angepasst (im Gegensatz zu lm() und aov(), die mit der Methode der kleinsten Quadrate angepasst werden).
Wir laden zuerst lme4 und lmerTest
library(lme4)
library(lmerTest) #immer laden nach lme4Die Syntax für ein Random Intercept Modell ist
m.2rep<-lmer(response~time+(1|subject), data=d.long2)Wir brauchen die lmer() Funktion. Das neue Element in der Modellformel ist der Term (1|subject). Dieser stellt ein Random Intercept dar (Abstand für Einheit \(i\) zum fixed Intercept). Der Grund ist inzwischen bekannt: Die Beobachtungen innerhalb eines Subjekts sind gruppiert (wiederholte Messungen innerhalb eines Subjekts). Das Random Intercept induziert dann eine Korrelation der wiederholten Messungen.
summary() ist eine generische Funktion, der auch lmer Objekte übergeben werden können:
summary(m.2rep,cor=FALSE)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: response ~ time + (1 | subject)
Data: d.long2
REML criterion at convergence: 408.5
Scaled residuals:
Min 1Q Median 3Q Max
-2.0380 -0.4875 -0.0289 0.5657 2.1045
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 60.07 7.751
Residual 24.57 4.956
Number of obs: 60, groups: subject, 30
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 24.73 1.68 38.57 14.726 < 2e-16 ***
timet2 10.16 1.28 29.00 7.937 9.4e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Im Output sind neu die Between-subject und Within-subject Varianz der Random effects ersichtlich. Das sind die Parameter, die die Verteilung der Random effects spezifizieren. Die Random Effects selber, die \(U_i\), sehen wir hier nicht. Ersichtlich ist auch die Anzahl Beobachtungen und die Anzahl statistischer Einheiten. Die Fixed effects sind analog zum Linearen Modell zu interpretieren. Die geschätzte erwartete Veränderung ist 10.1573458.
Die geschätzten Random Effects \(\hat{U}_i\) selber können wir extrahieren mit ranef(). Sie sind posterior-Schätzungen bedingt auf die Daten, \(\hat{U}_i=E(U_i\mid Y_i)\). Ihre Varianz ist kleiner als die wahre Varianzkomponente aufgrund von sogenannter Shrinkage (auf die wir hier nicht eingehen), \(\Var(\hat{U}_i)<\nu^2\).
ranef(m.2rep)$subject
(Intercept)
s1 0.7234120
s10 8.9460286
s11 -1.5461913
s12 -2.7086903
s13 2.4543710
s14 -3.6550798
s15 -0.8162158
s16 -2.1832215
s17 7.6531476
s18 -4.1532914
s19 -3.3900425
s2 -12.7239664
s20 -5.1947473
s21 8.5356755
s22 -3.8224224
s23 8.0093269
s24 6.2283789
s25 -0.4167840
s26 -3.7005501
s27 2.6031749
s28 -0.8796224
s29 -15.3696272
s3 4.7682908
s30 7.0648977
s4 2.8367827
s5 7.7111655
s6 5.3103447
s7 -15.0365021
s8 -8.8412370
s9 11.5931946
with conditional variances for "subject" mean(ranef(m.2rep)$subject$"(Intercept)")[1] -1.599999e-13var(ranef(m.2rep)$subject$"(Intercept)")[1] 49.875In gemischten linearen Modellen (LMMs) ist die Verteilung der Teststatistiken für feste Effekte komplizierter als in klassischen linearen Modellen, da die Schätzung der Varianzkomponenten (zufällige Effekte) zusätzliche Unsicherheit einführt.
Angenommen, wir testen einen festen Effekt \(\beta_j\). Die Teststatistik lautet: \[ t = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \]
Da die Standardabweichung von geschätzten Varianzkomponenten abhängt, ist die Verteilung von \(t\) nicht exakt eine \(t\)-Verteilung mit bekannten Freiheitsgraden. Die Satterthwaite-Methode approximiert die Freiheitsgrade \(\nu\) sodass:
\[ t \sim t_\nu \]
Die Approximation basiert auf der Idee, die Varianz der Schätzung selbst als Zufallsvariable zu behandeln und deren Unsicherheit zu berücksichtigen.
Die lmerTest-Erweiterung zu lme4 berechnet \(p\)-Werte und Freiheitsgrade für feste Effekte. Die Freiheitsgrade sind nicht ganzzahlig und hängen vom jeweiligen Effekt und der Struktur der zufälligen Effekte ab. Die Methode ist besonders nützlich bei unausgeglichenen Designs oder komplexen Random-Strukturen, wo klassische df-Schätzungen versagen. Das Laden von lmerTest nach dem Laden von lme4 wird empfohlen (und das haven wir auch gemacht), wenn \(p\)-Werte benötigt werden, da es Freiheitsgrade approximiert und die \(t\)-Statistiken somit für Hypothesentests interpretierbarer macht.
Den Zeiteffekt explizit als Modellvergleich machen wir mit einem Modellvergleich:
m.2rep0<-update(m.2rep,.~.-time)
anova(m.2rep0,m.2rep)##LR-TestData: d.long2
Models:
m.2rep0: response ~ (1 | subject)
m.2rep: response ~ time + (1 | subject)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
m.2rep0 3 454.16 460.44 -224.08 448.16
m.2rep 4 421.52 429.90 -206.76 413.52 34.633 1 3.98e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Alternativen wären drop1() für Wald-Tests und einfach anova(), weil wir nur eine Eingangsgrösse haben:
drop1(m.2rep) ##Wald-typeSingle term deletions using Satterthwaite's method:
Model:
response ~ time + (1 | subject)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 1547.6 1547.6 1 29 62.996 9.4e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m.2rep)Type III Analysis of Variance Table with Satterthwaite's method
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 1547.6 1547.6 1 29 62.996 9.4e-09 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Die Intraklassenkorrelation hatten wir oben schon berechnet.
print(VarCorr(m.2rep),comp="Variance") Groups Name Variance
subject (Intercept) 60.073
Residual 24.566 Die geschätzte Intraklassenkorrelation ist also
\[ \widehat{\Cor}(Y_{ij_1},Y_{ij_2})=\frac{\hat{\nu}^2}{\hat{\nu}^2+\hat{\sigma}^2}=0.7097532 \]
In diesem Fall gibt es nur einen Kontrast von Interesse, der selber schon ein Parameter der Grundparametrisierung ist:
library(emmeans)
summary(emmeans(m.2rep,revpairwise~time),infer=TRUE)$emmeans
time emmean SE df lower.CL upper.CL t.ratio p.value
t1 24.7 1.68 38.6 21.3 28.1 14.726 <.0001
t2 34.9 1.68 38.6 31.5 38.3 20.773 <.0001
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df lower.CL upper.CL t.ratio p.value
t2 - t1 10.2 1.28 29 7.54 12.8 7.937 <.0001
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95 Hier bringt dieser Schritt keine zusätzliche Information. Mit dem emmeans() kann man sich dann vor allem bei komplexeren Modellen Kontraste von Interesse berechnen lassen, wie wir das früher schon gemacht haben.
LMM sind zweistufige Modelle, man hat bedingte (auf Subject bedingt) und marginale Vorhersagen. Wir können ein Data frame erstellen mit den bedingten und marginalen Vorhersagen und diese plotten:
predictions <- data.frame(
actual = d.long2$response,
CondPred = predict(m.2rep), # PREDICT AT SUBJECT LEVEL
MargPred=predict(m.2rep,re.form=~0), # PREDICT MARGINALLY
time = d.long2$time,
subject=d.long2$subject)Wir sehen, dass die marginalen Vorhersagen nicht von subject abhängen.
headTail(predictions,8,8) actual CondPred MargPred time subject
1 22.93 25.46 24.73 t1 s1
2 38.43 35.62 34.89 t2 s1
3 10.8 12.01 24.73 t1 s2
4 18.17 22.17 34.89 t2 s2
5 33.22 29.5 24.73 t1 s3
6 37.9 39.66 34.89 t2 s3
7 24.82 27.57 24.73 t1 s4
8 41.64 37.73 34.89 t2 s4
... ... ... ... <NA> <NA>
53 29.69 27.34 24.73 t1 s27
54 36.21 37.49 34.89 t2 s27
55 23.93 23.85 24.73 t1 s28
56 33.57 34.01 34.89 t2 s28
57 6.53 9.36 24.73 t1 s29
58 16.07 19.52 34.89 t2 s29
59 34.25 31.8 24.73 t1 s30
60 42.4 41.96 34.89 t2 s30ggplot(predictions, aes(x = time, y = CondPred, color = subject)) +
geom_point(aes(y=actual,color=subject))+
geom_line(aes(y=CondPred,group=subject),lty=1)+
geom_line(aes(y=MargPred,group=subject),lty=2,lwd=1.3,col="black")+theme(legend.position="none")+ylab("prediction")
Eine Alternative für bedingte Vorhersagen wäre
lattice::xyplot(fitted(m.2rep)~time,groups=subject,data=d.long2,type="b")
##oder
d.longG <- nlme::groupedData(fitted(m.2rep) ~ time | subject, data = d.long2)
plot(d.longG)
Was passiert bei Nicht-Berücksichtigung der korrelierten Fehler? Die Punktschätzungen werden nicht beeinflusst, aber die Standardfehler, und damit die Resultate von Intervallschätzungen und Tests:
Der Effekt von between-subject oder between-cluster Variablen wird überschätzt. (Grund: Künstliche Erhöhung der Stichprobengröße).
Der Effekt von within-subject oder within-cluster Variablen wird unterschätzt. (Grund: Die Variabilität wird überschätzt, z. B. hat ein gepaarter \(t\)-Test mehr Power als ein unabhängiger \(t\)-Test).
Vergleich wir also mit den Resultate, wenn wir statt ein LMM ein LM (lm()) anpassen würden:
summary(lm(response~time,data=d.long2))$coef ## false Estimate Std. Error t value Pr(>|t|)
(Intercept) 24.73446 1.679675 14.725746 2.976139e-21
timet2 10.15735 2.375419 4.276023 7.190056e-05Wir sehen, dass der within-subject Effekt dann einen grösseren Standardfehler hat und damit einen kleineren \(t\)-Wert als im LMM:
summary(m.2rep)$coef ## correct Estimate Std. Error df t value Pr(>|t|)
(Intercept) 24.73446 1.679675 38.57025 14.725746 2.025810e-17
timet2 10.15735 1.279746 29.00000 7.936999 9.399616e-09Wir verallgemeinern nun den gepaarten \(t\)-Test auf ein within-subject Design.
Das Data frame d.long besteht 80 Beobachtungen aus \(n_i=4\) wiederholten Messungen (time) von \(n=20\) subjects (subject) auf einer kontinuierlichen abhängigen Variablen (response).
Hier haben wir \(n_i=4\) konstant für alle \(i\). Die Anzahl Beobachtungen pro Einheit muss aber nicht konstant sein. Wir könnten ein LMM auch dann anpassen, wenn wir für jede Einheit \(i\) eine unterschiedliche Anzahl Beobachtungen hätten (z.B. \(n_1=2,n_2=5,\dots,n_n=3\)).
d.long<-read.csv("https://raw.githubusercontent.com/mcdr65/StatsRsource/master/Data/ToyRepAnova.csv",
stringsAsFactors=TRUE)
str(d.long)'data.frame': 80 obs. of 3 variables:
$ subject : Factor w/ 20 levels "s1","s10","s11",..: 1 1 1 1 12 12 12 12 14 14 ...
$ time : Factor w/ 4 levels "t1","t2","t3",..: 1 2 3 4 1 2 3 4 1 2 ...
$ response: num 53 43 66.2 62.5 26.4 ...psych::headTail(d.long) subject time response
1 s1 t1 52.98
2 s1 t2 43
3 s1 t3 66.25
4 s1 t4 62.47
... <NA> <NA> ...
77 s20 t1 63.77
78 s20 t2 23.15
79 s20 t3 54.46
80 s20 t4 26.16with(d.long,xtabs(~time))time
t1 t2 t3 t4
20 20 20 20 Zuerst beschreiben wir die Daten.
with(d.long,interaction.plot(time,subject,response,ylab="response"))
## alternative with ggplot2
p <- ggplot(data = d.long, aes(x = time, y = response, group = subject))
p <- p+geom_point()+geom_line(lty=2)+stat_summary(fun=mean,geom="line",aes(group = 1),color="red",lwd=2)
p <- p + stat_summary(aes(group=1), geom = "point", fun = mean,shape = 17, size = 4)
p
Wir können jeden Zeitpunkt einzeln beschreiben:
psych::describeBy(d.long$response,group=d.long$time,mat=TRUE) item group1 vars n mean sd median trimmed mad min
X11 1 t1 1 20 26.58151 18.84398 26.57575 26.43457 19.55346 -9.2323379
X12 2 t2 1 20 39.01437 21.47763 36.43885 39.34613 19.03513 -0.5464301
X13 3 t3 1 20 44.68697 20.85515 47.43915 43.59238 23.66503 11.9344197
X14 4 t4 1 20 36.89966 21.76513 32.50413 35.17075 19.96541 6.8175026
max range skew kurtosis se
X11 63.77199 73.00433 0.10434835 -0.7370743 4.213641
X12 72.57076 73.11719 0.06215493 -1.0902151 4.802544
X13 96.31308 84.37866 0.44954202 -0.2694602 4.663354
X14 81.48171 74.66421 0.61472843 -0.6965030 4.866831tapply(d.long$response,d.long$time,summary) #alternative$t1
Min. 1st Qu. Median Mean 3rd Qu. Max.
-9.232 10.966 26.576 26.582 35.109 63.772
$t2
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.5464 24.2050 36.4389 39.0144 52.7406 72.5708
$t3
Min. 1st Qu. Median Mean 3rd Qu. Max.
11.93 27.87 47.44 44.69 55.25 96.31
$t4
Min. 1st Qu. Median Mean 3rd Qu. Max.
6.817 22.668 32.504 36.900 48.607 81.482 aggregate(response~time,data=d.long,summary) #alternative time response.Min. response.1st Qu. response.Median response.Mean
1 t1 -9.2323379 10.9655202 26.5757505 26.5815129
2 t2 -0.5464301 24.2050217 36.4388469 39.0143684
3 t3 11.9344197 27.8672784 47.4391539 44.6869727
4 t4 6.8175026 22.6681781 32.5041280 36.8996576
response.3rd Qu. response.Max.
1 35.1086532 63.7719879
2 52.7406216 72.5707588
3 55.2541992 96.3130758
4 48.6071278 81.4817148Wir können auch mit dem Wide-Format arbeiten. Manchmal ist das für die deskriptive Anlayse einfacher.
d.wide<-reshape(d.long,v.names="response",idvar="subject",timevar="time",direction="wide")
headTail(d.wide) subject response.t1 response.t2 response.t3 response.t4
1 s1 52.98 43 66.25 62.47
5 s2 26.43 26.45 57.62 47.31
9 s3 10.35 71.25 49.91 64.49
13 s4 10.88 13.85 61.11 32.41
... <NA> ... ... ... ...
65 s17 11 23.05 32.81 23.58
69 s18 26.93 33.05 11.93 18.14
73 s19 22.18 48.82 50.68 32.6
77 s20 63.77 23.15 54.46 26.16summary(d.wide[,-1]) response.t1 response.t2 response.t3 response.t4
Min. :-9.232 Min. :-0.5464 Min. :11.93 Min. : 6.817
1st Qu.:10.966 1st Qu.:24.2050 1st Qu.:27.87 1st Qu.:22.668
Median :26.576 Median :36.4389 Median :47.44 Median :32.504
Mean :26.582 Mean :39.0144 Mean :44.69 Mean :36.900
3rd Qu.:35.109 3rd Qu.:52.7406 3rd Qu.:55.25 3rd Qu.:48.607
Max. :63.772 Max. :72.5708 Max. :96.31 Max. :81.482 Die Syntax ist äquivalent zum gepaarten \(t\)-Test. Wir nennen das neue Modellobjekt m.4rep:
m.4rep<-lmer(response~time+(1|subject), data=d.long)
summary(m.4rep,cor=FALSE)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: response ~ time + (1 | subject)
Data: d.long
REML criterion at convergence: 684
Scaled residuals:
Min 1Q Median 3Q Max
-1.8422 -0.5911 -0.1200 0.6007 2.2892
Random effects:
Groups Name Variance Std.Dev.
subject (Intercept) 97.31 9.865
Residual 333.95 18.274
Number of obs: 80, groups: subject, 20
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 26.582 4.644 65.930 5.724 2.77e-07 ***
timet2 12.433 5.779 57.000 2.151 0.03569 *
timet3 18.105 5.779 57.000 3.133 0.00273 **
timet4 10.318 5.779 57.000 1.786 0.07950 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Neue sind nun zusätzliche feste Effekte, da Zeitpunkt nun vierwertig ist.
Der Globaltest for den Effekt von time mit drop1() zeigt, dass man time nicht aus dem Modell nehmen darf.
drop1(m.4rep) ##Wald-typeSingle term deletions using Satterthwaite's method:
Model:
response ~ time + (1 | subject)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
time 3430.7 1143.6 3 57 3.4244 0.02305 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Die Between-und within-subject Varianz ist im Output ersichtlich, kann aber auch extrahiert werden:
print(VarCorr(m.4rep),comp="Variance") Groups Name Variance
subject (Intercept) 97.31
Residual 333.95 Die geschätzte within-subject (within-unit) Korrelation – also die Intraklassenkorrelation – ist:
\(\widehat{\Cor}(Y_{ij_1},Y_{ij_2})=\frac{\hat{\sigma}^2_{U}}{\hat{\sigma}^2_{U}+\hat{\sigma}^2_{\epsilon}}=0.2256404\)
Wir schauen uns noch den Tukey-Anscombe Plot an und den quantile–quantile plot an. Die Varianz der Residuen sollte homogen sein und sie sollten normalverteilt sein.
plot(m.4rep,type=c("p","smooth"))
qqnorm(residuals(m.4rep))
Es gibt hier keine klare Evidenz gegen die Modellannahmen.
Im summary()-Output sehen wir wie die geschätzten Zeiteffekte relativ zur Referenzkategorie. Wenn wir alle paarweisen Kontraste wollen (mit Konfidenzintervallen), können wir wieder emmeans brauchen.
summary(emmeans(m.4rep,revpairwise~time,infer=TRUE))$emmeans
time emmean SE df lower.CL upper.CL t.ratio p.value
t1 26.6 4.64 65.9 17.3 35.9 5.724 <.0001
t2 39.0 4.64 65.9 29.7 48.3 8.402 <.0001
t3 44.7 4.64 65.9 35.4 54.0 9.623 <.0001
t4 36.9 4.64 65.9 27.6 46.2 7.946 <.0001
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df lower.CL upper.CL t.ratio p.value
t2 - t1 12.43 5.78 57 -2.86 27.73 2.151 0.1495
t3 - t1 18.11 5.78 57 2.81 33.40 3.133 0.0141
t3 - t2 5.67 5.78 57 -9.62 20.97 0.982 0.7604
t4 - t1 10.32 5.78 57 -4.98 25.61 1.786 0.2908
t4 - t2 -2.11 5.78 57 -17.41 13.18 -0.366 0.9831
t4 - t3 -7.79 5.78 57 -23.08 7.51 -1.348 0.5370
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
Conf-level adjustment: tukey method for comparing a family of 4 estimates
P value adjustment: tukey method for comparing a family of 4 estimates Wir kommen jetzt zu einem Design, das in den Gesundheitsforschung sehr häufig vorkommt. Wir haben jetzt einen between-subject Faktor group \(G\) und eine within-subject Variable time \(T\).
Es werden also zwei Gruppen \(A\) und \(B\) (Levels des between-subject Faktors group) über die Zeit (within) beobachtet. Im Folgenden behandeln wir aber \(T\) nicht als Faktor, sondern als eine kontinuierliche Eingangsgrösse. time \(T\) könnte aber auch als Faktor behandelt werden (Zeitpunkt statt Zeit).
Wir haben jetzt einen linearen Prädiktor mit vier Parametern \(\beta_0,\dots,\beta_3\), da jede Gruppe eine Lage und eine Steigung bezüglich Zeit hat. Das Random Intercept Modell ist
\[ Y_{ij}=\underbrace{\beta_0+\beta_1 G_{ij} + \beta_2 T_{ij}+\beta_3 GT_{ij}}_{fixed}+\underbrace{\overbrace{U_i}^{rand.intercept}+\overbrace{\epsilon_{ij}}^{resid.error}}_{random} \]
mit
\[G_{ij}=\begin{cases} 1 &\mbox{if }G_{ij}=B\\ 0 & \mbox{if } G_{ij}=A \end{cases}, \qquad GT_{ij}=\begin{cases} T_{ij} &\mbox{if }G_{ij}=B\\ 0 & \mbox{if } G_{ij}=A \end{cases} \]
Marginale Erwartung: Schauen wir uns die Erwartungswerte an. Für die marginale Erwartung haben wir – da alle random effects im Schnitt Null sind – einfach den linearen Prädiktor \(\beta_0+\beta_1 G_{ij} + \beta_2 T_{ij}+\beta_3 GT_{ij}\). Für Gruppe \(A\) und \(B\) gibt das:
\[ \begin{align} A: E(Y_{ij})=&\beta_0+ \beta_2 T_{ij}\\ B: E(Y_{ij})=&(\beta_0+\beta_1) + (\beta_2+\beta_3)T_{ij}\\ \end{align} \]
Wir sehen hier schön, dass jede Gruppe eine eigene Lage und eine eigene Steigung hat. Der Parameter von Interesse ist häufig der Interaktionsparameter\(\beta_3\). Dieser ist als Unterschied der Steigungen zwischen den beiden Gruppen zu interpretieren. Interaktion bedeutet: Der Effekt von \(G\) ist dann \(\beta_1+\beta_3 T\), also eine Funktion der Zeit! Oder: Der Effekt von \(T\) ist dann \(\beta_2+\beta_3 G\), also eine Funktion der Gruppe!
Die \(U_i\)’s interessieren häufig nicht, sie müssen aber im Modell sein, um der Korrelationsstruktur der Daten Rechnung zu tragen. Ein Erweiterung des Modells ist das Random Slope Modell, dieses erlaubt zusätzlich zum Random Intercept Modell noch unterschiedliche Slopes für jedes Subject. Die Interpretation der fixed effects ändert sich aber nicht. Das Random Slope Modell ist
\[ Y_{ij}=\underbrace{\beta_0+\beta_1 G_{ij} + \beta_2 T_{ij}+\beta_3 GT_{ij}}_{fixed}+\underbrace{\overbrace{U_{i0}}^{rand.intercept}+T_{ij}\overbrace{U_{i1}}^{rand.slope}+\overbrace{\epsilon_{ij}}^{resid.error}}_{random} \]
Für die Interessierten: Mit dem folgenden Code simulieren wir Daten aus dem besprochenen Design: \(n=100\) statistische Einheiten aus den Gruppen \(A\) und \(B\) werden. Jede Einheit hat vier Beobachtungen über die Zeit. Die abhängige Variable ist \(Y\). Natürlich kennen wir bei realen Daten die wahren Parameter nicht, aber es lohnt sich, Daten aus einem Modell zu simulieren und dann die Analyse zu machen. Dies ist sehr hilfreich für das Verständnis des Designs, des Modells und schliesslich der Anlayse.
simulate <- function(seed=62) {
n <-100 # number of subjects
K <- 4 # number of measurements per subject
# data frame with the design: everyone has a baseline measurement, and then measurements at random follow-up times
d.SP <- data.frame(id = rep(seq_len(n), each = K), time = rep(1:K, n), group = rep(gl(2, n/2, labels = c("A", "B")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~group * time, data = d.SP)
Z_base <- cbind(1, d.SP$time[1:K]) # intercept and time columns for one subject
Z <- kronecker(diag(n), Z_base) # Block-diagonal matrix for all subjects
beta <- c(10, 3, 1, 4) # fixed effects coefficients
##
var_intercept <- 1 # variance of random intercepts
var_slope <- 1 # variance of random slopes
correlation <- 0 # correlation between intercept and slope
D <- matrix(c(var_intercept, correlation * sqrt(var_intercept * var_slope),
correlation * sqrt(var_intercept * var_slope), var_slope), nrow=2)
sigma <- 2 # sd error terms
# Simulate random effects for intercept and slope
set.seed(seed)
B <- MASS::mvrnorm(n, rep(0, ncol(Z_base)), D)
b <- as.vector(t(B))
# linear predictor
eta <- X %*% beta + Z%*%b
# we simulate continuous longitudinal data
d.SP$y <- rnorm(n * K, mean = eta, sd = sigma)
list(data=d.SP,parms=list(random=data.frame(varInt=1,varSlope=1,Corr=0),fixed=data.frame(beta0=10,beta1=3,beta2=1,beta3=4),sigma2=sigma^2))
}
simSP<-simulate()
d.SP<-simSP$dataWir haben ein group(2) \(\times\) time(4) between-within Design, mit \(n=100\) (50 pro Gruppe mit je vier Beobachtungen). Wir behandeln Zeit als eine kontinuierlich Eingangsgrösse!
str(d.SP)'data.frame': 400 obs. of 4 variables:
$ id : int 1 1 1 1 2 2 2 2 3 3 ...
$ time : int 1 2 3 4 1 2 3 4 1 2 ...
$ group: Factor w/ 2 levels "A","B": 1 1 1 1 1 1 1 1 1 1 ...
$ y : num 12.5 12.9 13.7 18.8 10.2 ...with(d.SP,xtabs(~group+time)) time
group 1 2 3 4
A 50 50 50 50
B 50 50 50 50Da wir die Daten aus einem Modell simuliert haben, kennen wir hier die wahren Werte der Parameter. Unter $fixed sehen wir die Parameter der fixed effects. Unter $random sehen wir die Varianz der Random Intercepts, der Random Slopes und die Korrelation dieser (die wir als Null angenommen haben). Unter $sigma2 sehen wir die Residualvarianz.
simSP$parms$random
varInt varSlope Corr
1 1 1 0
$fixed
beta0 beta1 beta2 beta3
1 10 3 1 4
$sigma2
[1] 4Die (aus einem bekannten Modell simulierten) Daten im Long Format (400 Zeilen) sind dann
psych::headTail(d.SP,4,4) id time group y
1 1 1 A 12.48
2 1 2 A 12.93
3 1 3 A 13.71
4 1 4 A 18.82
... ... ... <NA> ...
397 100 1 B 16.04
398 100 2 B 22.9
399 100 3 B 29.33
400 100 4 B 32.19und im Wide Format (100 Zeilen), oft nützlich für deskriptive Statistik:
d.SPwide<-reshape(d.SP,v.names="y",timevar="time",idvar="id",direction="wide")
psych::headTail(d.SPwide,4,4) id group y.1 y.2 y.3 y.4
1 1 A 12.48 12.93 13.71 18.82
5 2 A 10.16 14.46 9.18 14.67
9 3 A 7.27 9.44 6.97 6.88
13 4 A 17.14 15.11 16.96 15.81
... ... <NA> ... ... ... ...
385 97 B 16.94 21.71 20.28 31.72
389 98 B 19.42 20.78 32.12 34.45
393 99 B 19.31 25.29 31.9 31.84
397 100 B 16.04 22.9 29.33 32.19Hier wieder mit mehreren Varianten:
with(d.SP,interaction.plot(x.factor=time,trace.factor=group,response=y))
lattice::xyplot(data = d.SP, y ~ time | group,groups = id, type = "l")
tapply(d.SP$y, list(d.SP$group,d.SP$time), mean) 1 2 3 4
A 11.10430 12.43844 13.20379 14.75693
B 18.07467 23.24057 28.65372 34.07037psych::describeBy(d.SP$y, list(d.SP$group,d.SP$time), mat = TRUE) item group1 group2 vars n mean sd median trimmed mad
X11 1 A 1 1 50 11.10430 2.714649 10.60099 10.96877 2.580885
X12 2 B 1 1 50 18.07467 2.079430 18.29829 18.16031 1.603988
X13 3 A 2 1 50 12.43844 2.354138 12.74051 12.42869 2.910096
X14 4 B 2 1 50 23.24057 3.072497 22.84238 23.09342 3.012479
X15 5 A 3 1 50 13.20379 4.126539 13.05828 13.26863 4.910530
X16 6 B 3 1 50 28.65372 3.966173 28.56873 28.51601 4.304457
X17 7 A 4 1 50 14.75693 4.251634 14.38053 14.70346 4.216739
X18 8 B 4 1 50 34.07037 4.445264 33.85766 33.86014 3.810761
min max range skew kurtosis se
X11 6.086921 17.60898 11.522060 0.46265396 -0.09545014 0.3839093
X12 13.001167 22.85525 9.854087 -0.33820447 0.13970097 0.2940758
X13 7.400169 17.32783 9.927660 0.01630782 -0.81882901 0.3329254
X14 17.528504 31.71964 14.191136 0.50705979 -0.21230877 0.4345167
X15 3.591481 21.81274 18.221256 -0.11984030 -0.78636981 0.5835807
X16 20.279260 38.70488 18.425620 0.20702813 -0.36270868 0.5609015
X17 6.877855 22.92871 16.050854 0.09403644 -0.97331145 0.6012718
X18 24.810750 45.63598 20.825229 0.43800524 0.13648284 0.6286553aggregate(y~time+group,data=d.SP,summary)##Alternative time group y.Min. y.1st Qu. y.Median y.Mean y.3rd Qu. y.Max.
1 1 A 6.086921 9.327115 10.600989 11.104296 12.848087 17.608980
2 2 A 7.400169 10.410568 12.740514 12.438442 14.324042 17.327829
3 3 A 3.591481 10.185784 13.058277 13.203787 16.477755 21.812738
4 4 A 6.877855 11.556186 14.380528 14.756929 17.703400 22.928708
5 1 B 13.001167 17.025762 18.298287 18.074673 19.320991 22.855254
6 2 B 17.528504 21.137155 22.842381 23.240574 25.285105 31.719640
7 3 B 20.279260 25.879617 28.568727 28.653717 31.480437 38.704880
8 4 B 24.810750 31.304031 33.857658 34.070374 36.346075 45.635979tapply(d.SPwide,d.SPwide$group,summary)##mit d.SPwide$A
id group y.1 y.2 y.3
Min. : 1.00 A:50 Min. : 6.087 Min. : 7.40 Min. : 3.591
1st Qu.:13.25 B: 0 1st Qu.: 9.327 1st Qu.:10.41 1st Qu.:10.186
Median :25.50 Median :10.601 Median :12.74 Median :13.058
Mean :25.50 Mean :11.104 Mean :12.44 Mean :13.204
3rd Qu.:37.75 3rd Qu.:12.848 3rd Qu.:14.32 3rd Qu.:16.478
Max. :50.00 Max. :17.609 Max. :17.33 Max. :21.813
y.4
Min. : 6.878
1st Qu.:11.556
Median :14.381
Mean :14.757
3rd Qu.:17.703
Max. :22.929
$B
id group y.1 y.2 y.3
Min. : 51.00 A: 0 Min. :13.00 Min. :17.53 Min. :20.28
1st Qu.: 63.25 B:50 1st Qu.:17.03 1st Qu.:21.14 1st Qu.:25.88
Median : 75.50 Median :18.30 Median :22.84 Median :28.57
Mean : 75.50 Mean :18.07 Mean :23.24 Mean :28.65
3rd Qu.: 87.75 3rd Qu.:19.32 3rd Qu.:25.29 3rd Qu.:31.48
Max. :100.00 Max. :22.86 Max. :31.72 Max. :38.70
y.4
Min. :24.81
1st Qu.:31.30
Median :33.86
Mean :34.07
3rd Qu.:36.35
Max. :45.64 Wir haben als fixed effects: group, time und group:time-Interaktion. Wir haven zwei verschiedene Random-Effects-Strukturen (1|subject) (für das random intercept Modell) und (time|subject) (für das random slope Modell):
m1 <- lmer(y ~ group * time + (1 | id), data = d.SP)time|idm2 <- lmer(y ~ group * time + (time | id), data = d.SP)Versuchen Sie, die Parameterschätzungen der festen Effekte unten in der rechten Graphik nachzuvollziehen. So ist z.B. \(\hat{\beta}_3=4.1677004\) der geschätzte erwartete Unterschied in der Steigung zwischen Gruppe \(A\) und \(B\), \(\hat{\beta}_2=1.1723243\) ist die geschätzte erwartete Steigung für Gruppe \(A\), usw.
summary(m1, cor = FALSE)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: y ~ group * time + (1 | id)
Data: d.SP
REML criterion at convergence: 1985.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.97793 -0.54881 0.00251 0.58228 2.78736
Random effects:
Groups Name Variance Std.Dev.
id (Intercept) 6.778 2.603
Residual 5.342 2.311
Number of obs: 400, groups: id, 100
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 9.9451 0.5439 266.3516 18.285 < 2e-16 ***
groupB 2.7147 0.7692 266.3516 3.529 0.000491 ***
time 1.1723 0.1462 298.0000 8.020 2.44e-14 ***
groupB:time 4.1677 0.2067 298.0000 20.161 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predicted <- predict(m1)
lattice::xyplot(data = d.SP, y ~ time | group,groups = id, type = "l",main="Daten")
lattice::xyplot(data = d.SP, predicted ~ time | group, groups = id, type = "l",main="Modell")
Für das random slope Modell sehen wir zusätzlich eine Variabilität in den Steigungen für jede statistische Einheit. Die Punktschätzungen für die fixed effects bleiben aber gleich.
summary(m2, cor = FALSE)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: y ~ group * time + (time | id)
Data: d.SP
REML criterion at convergence: 1908.9
Scaled residuals:
Min 1Q Median 3Q Max
-2.9705 -0.5587 0.0275 0.5456 2.2772
Random effects:
Groups Name Variance Std.Dev. Corr
id (Intercept) 1.600 1.265
time 1.015 1.007 -0.12
Residual 3.673 1.917
Number of obs: 400, groups: id, 100
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 9.9451 0.3771 97.9991 26.374 < 2e-16 ***
groupB 2.7147 0.5333 97.9991 5.091 1.73e-06 ***
time 1.1723 0.1871 98.0009 6.267 9.87e-09 ***
groupB:time 4.1677 0.2645 98.0009 15.755 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1predicted <- predict(m2)
lattice::xyplot(data = d.SP, y ~ time | group,groups = id, type = "l",main="Daten")
lattice::xyplot(data = d.SP, predicted ~ time | group, groups = id, type = "l",main="Modell")
Modell m1 ist ja verschachtelt im Modell m2. Wir können also die Hypothese testen, dass die Variabilität der random slopes Null ist:
anova(m1, m2)Data: d.SP
Models:
m1: y ~ group * time + (1 | id)
m2: y ~ group * time + (time | id)
npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq)
m1 6 1993.8 2017.7 -990.89 1981.8
m2 8 1920.8 1952.8 -952.41 1904.8 76.963 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1oder
lmerTest::ranova(m2)ANOVA-like table for random-effects: Single term deletions
Model:
y ~ group + time + (time | id) + group:time
npar logLik AIC LRT Df Pr(>Chisq)
<none> 8 -954.44 1924.9
time in (time | id) 6 -992.89 1997.8 76.908 2 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Das Random Intercept Modell wird also zugunsten des Random Slope Modells verworfen.
Mit einem Caterpillar Plot können wir die Random Effects darstellen.
dotplot.ranef.mer(ranef(m2))$id
Das Hauptinteresse liegt meistens im Testen der fixed effects. Zuerst der Interaktionseffekt:
drop1(m2)Single term deletions using Satterthwaite's method:
Model:
y ~ group * time + (time | id)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
group:time 911.66 911.66 1 98.001 248.21 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1drop1() testet richtigerweise nur die Interaktion. Haupteffekte sollten in der Gegenwart von Interaktion nicht getestet werden.
Type III Tests: Macht man es trotzdem, muss man sogenannte Type III Tests und den Zusammenhang mit der Default-Effektkodierung sehr gut verstehen. Wir gehen hier nicht näher darauf ein. Die Funktion anova.lmerModLmerTest() aus lmerTest macht diese Tests korrekt bei Aufruf von anova() und Übergabe eines lmer-Objekts und ist unabhängig von der Parametrisierung (im Gegensatz zu car::Anova(...,type=3)).
Die Effektparameterisierung gab uns für die fixed effects:
summary(m2)$coef Estimate Std. Error df t value Pr(>|t|)
(Intercept) 9.945053 0.3770772 97.99907 26.374053 2.650476e-46
groupB 2.714720 0.5332676 97.99907 5.090727 1.725057e-06
time 1.172324 0.1870577 98.00088 6.267181 9.871746e-09
groupB:time 4.167700 0.2645395 98.00088 15.754548 1.313324e-28Wir wollen jetzt eine Schätzung für die Zeiteffekte in beiden Gruppen \(A\) und \(B\) und eine Schätzung des Unterschieds des Zeiteffekts zwischen Gruppe \(A\) und \(B\). Dies bekommen wir die mit emmeans::emtrends()
summary(emmeans::emtrends(m2,trt.vs.ctrl~group,var="time"),infer=TRUE)$emtrends
group time.trend SE df lower.CL upper.CL t.ratio p.value
A 1.17 0.187 98 0.801 1.54 6.267 <.0001
B 5.34 0.187 98 4.969 5.71 28.547 <.0001
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95
$contrasts
contrast estimate SE df lower.CL upper.CL t.ratio p.value
B - A 4.17 0.265 98 3.64 4.69 15.755 <.0001
Degrees-of-freedom method: kenward-roger
Confidence level used: 0.95 Die Gruppe \(B\) verändert sich pro Zeiteinheit um 4.17 Einheiten auf \(Y\) mehr als die Gruppe \(A\). Die Steigung in der Gruppe \(A\) hatten wir schon in der Effektparameterisierung. Die Steigung in der Gruppe \(B\) ist 5.34, das entspricht in der Grundparametrisierung 1.172+4.168.
Hier gibt es keine Evidenz gegen die Modellannahmen. Das wäre auch erstaunlich, haben wir doch die Daten aus demselben Modell simuliert, welches wir dann an die Daten angepasst haben.
Streuung Residuen
plot(m2,resid(.)~fitted(.)|group)
plot(m2,resid(.)~fitted(.)|group*time)
Streuung und Verteilung Random effects
