Paired T-Test als Spezialfall der linearen Mixed-Effect-Modellierung

20

Wir wissen, dass ein gepaarter t- Test nur ein Sonderfall von Einweg-ANOVA mit wiederholten Messungen (oder innerhalb des Subjekts) sowie eines linearen Mischeffektmodells ist, das mit der Funktion lme () des nlme-Pakets in R demonstriert werden kann Wie nachfolgend dargestellt.

#response data from 10 subjects under two conditions
x1<-rnorm(10)
x2<-1+rnorm(10)

# Now create a dataframe for lme
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Wenn ich den folgenden gepaarten T-Test durchführe:

t.test(x1, x2, paired = TRUE)

Ich habe folgendes Ergebnis erhalten (Sie erhalten aufgrund des Zufallsgenerators ein anderes Ergebnis):

t = -2.3056, df = 9, p-value = 0.04657

Mit dem ANOVA-Ansatz können wir das gleiche Ergebnis erzielen:

summary(aov(y ~ x + Error(subj/x), myDat))

# the F-value below is just the square of the t-value from paired t-test:
          Df  F value Pr(>F)
x          1  5.3158  0.04657

Jetzt kann ich dasselbe Ergebnis mit dem folgenden Modell erhalten, wenn ich eine positiv-definite symmetrische Korrelationsmatrix für die beiden Bedingungen annehme:

summary(fm1 <- lme(y ~ x, random=list(subj=pdSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.3142115  9 -0.7918878  0.4488
# xx2          1.3325786 0.5779727  9  2.3056084  0.0466

Oder ein anderes Modell, das eine zusammengesetzte Symmetrie für die Korrelationsmatrix der beiden Bedingungen annimmt:

summary(fm2 <- lme(y ~ x, random=list(subj=pdCompSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.4023431  9 -0.618428  0.5516
# xx2          1.3325786 0.5779727  9  2.305608  0.0466

Mit dem gepaarten t-Test und der Einweg-ANOVA mit wiederholten Messungen kann ich das traditionelle Zellenmittelwertmodell wie folgt aufschreiben

Yij = μ + αi + βj + εij, i = 1, 2; j = 1, ..., 10

wobei i die Bedingung indiziert, j das Subjekt indiziert, Y ij die Antwortvariable ist, μ für den festen Effekt für den Gesamtmittelwert konstant ist, α i der feste Effekt für die Bedingung ist, β j der zufällige Effekt für das Subjekt ist, das N (0, σ folgt p 2 ) (σ p 2 ist die Populationsvarianz) und ε ij ist der Rest nach N (0, σ 2 ) (σ 2 ist die Varianz innerhalb des Subjekts).

Ich dachte, dass das obige Modell des Zellenmittelwerts nicht für die IME-Modelle geeignet wäre, aber das Problem ist, dass ich kein vernünftiges Modell für die beiden IME () -Ansätze mit der Annahme der Korrelationsstruktur finden kann. Der Grund ist, dass das IME-Modell mehr Parameter für die zufälligen Komponenten zu haben scheint, als das oben angegebene Zellenmittelwertmodell bietet. Zumindest liefert das IME-Modell genau den gleichen F-Wert, die gleichen Freiheitsgrade und den gleichen p-Wert, wie es gls nicht kann. Insbesondere liefert gls falsche DFs, da nicht berücksichtigt wird, dass jedes Subjekt zwei Beobachtungen hat, was zu stark aufgeblasenen DFs führt. Das Modell ist höchstwahrscheinlich bei der Angabe der Zufallseffekte überparametrisiert, aber ich weiß nicht, was das Modell ist und was die Parameter sind. Daher ist das Problem für mich immer noch ungelöst.

Bluepole
quelle
2
Nicht ganz sicher, was Sie fragen. Das Modell, das Sie notiert haben, ist genau das Modell für das Zufallseffektmodell. Die Korrelationsstruktur wird durch den Zufallseffekt induziert.
Aaron - Wiedereinsetzung von Monica
@Aaron: Der Zufallseffekt βj im Zellmittelwertmodell soll N (0, σp2) folgen. Meine Verwirrung ist, wie dieser Term (mit nur einem Parameter σp2) mit der Korrelationsstruktur assoziiert ist, die entweder durch eine zusammengesetzte Symmetrie oder eine einfache symmetrische Matrix im IME-Modell spezifiziert ist?
Bluepole
Wenn Sie die Korrelation zwischen den beiden Beobachtungen zu demselben Thema berechnen, lautet die Korrelation sigma_p ^ 2 / (sigma_p ^ 2 + sigma ^ 2), da sie dasselbe beta_j aufweisen. Siehe Pinheiro / Bates S.8. Außerdem entspricht das Zufallseffektmodell, wie Sie es geschrieben haben, der zusammengesetzten Symmetrie. andere Korrelationsstrukturen sind komplexer.
Aaron - Reinstate Monica
@ Aaron: Danke! Ich habe bereits das Pinheiro / Bates-Buch darüber gelesen und konnte die Einzelheiten zu den zufälligen Effekten immer noch nicht herausfinden. Die relevanteren Seiten scheinen das Beispiel auf S. 160-161 zu sein. Auch die Zufallseffekte, die von lme () mit der Annahme einer zusammengesetzten Symmetrie ausgegeben werden, scheinen nicht mit der Korrelation von σp2 / (σp2 + σ2) im Zellenmittelwertmodell übereinzustimmen. Immer noch ratlos über die Modellstruktur.
Bluepole
Nun, fast gleichbedeutend mit zusammengesetzter Symmetrie; In CS kann die Korrelation negativ sein, jedoch keine zufälligen Effekte. Vielleicht entsteht dort Ihr Unterschied. Weitere Informationen finden Sie unter stats.stackexchange.com/a/14185/3601 .
Aaron - Reinstate Monica

Antworten:

16

Die Äquivalenz der Modelle kann durch Berechnung der Korrelation zwischen zwei Beobachtungen desselben Individuums wie folgt beobachtet werden:

Y.ichj=μ+αich+βj+ϵichjβjN(0,σp2)ϵichjN(0,σ2)COv(yichk,yjk)=COv(μ+αich+βk+ϵichk,μ+αj+βk+ϵjk)=COv(βk,βk)=σp2Veinr(yichk)=Veinr(yjk)=σp2+σ2σp2/(σp2+σ2)

Beachten Sie, dass die Modelle nicht ganz gleichwertig sind, da das Zufallseffektmodell eine positive Korrelation erzwingt. Das CS-Modell und das T-Test / Anova-Modell nicht.

EDIT: Es gibt auch zwei andere Unterschiede. Erstens gehen das CS- und das Zufallseffektmodell von einer Normalität für den Zufallseffekt aus, das t-Test / Anova-Modell jedoch nicht. Zweitens werden die CS- und Zufallseffektmodelle mit maximaler Wahrscheinlichkeit angepasst, während die Anova mit mittleren Quadraten angepasst wird. Wenn alles im Gleichgewicht ist, werden sie übereinstimmen, aber nicht unbedingt in komplexeren Situationen. Schließlich wäre ich vorsichtig, wenn ich F / df / p-Werte aus den verschiedenen Anpassungen als Maß für die Übereinstimmung der Modelle verwenden würde. Weitere Einzelheiten finden Sie in Doug Bates 'berühmtem Estrich auf den DFs. (END EDIT)

Das Problem mit Ihrem RCode ist, dass Sie die Korrelationsstruktur nicht richtig angeben. Sie müssen glsmit dem verwendencorCompSymm Korrelationsstruktur verwenden.

Generieren Sie Daten, damit ein Betreff angezeigt wird:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Dann sehen Sie, wie Sie die Zufallseffekte und die zusammengesetzten Symmetriemodelle anpassen.

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

Die Standardfehler aus dem Zufallseffektmodell sind:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

Und die Korrelation und Restvarianz aus dem CS-Modell ist:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

Und sie entsprechen dem, was erwartet wird:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Andere Korrelationsstrukturen passen normalerweise nicht zu zufälligen Effekten, sondern nur durch Angabe der gewünschten Struktur. Eine häufige Ausnahme ist das AR (1) + Zufallseffektmodell, das einen Zufallseffekt und eine AR (1) -Korrelation zwischen Beobachtungen desselben Zufallseffekts aufweist.

EDIT2: Wenn ich die drei Optionen anpasse, erhalte ich genau die gleichen Ergebnisse, außer dass gls nicht versucht, den df für den gewünschten Begriff zu erraten.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(Der Achsenabschnitt ist hier anders, da bei der Standardcodierung nicht der Mittelwert aller Subjekte, sondern der Mittelwert des ersten Subjekts angegeben wird.)

Es ist auch interessant festzustellen, dass das neuere lme4Paket die gleichen Ergebnisse liefert, aber nicht einmal versucht, einen p-Wert zu berechnen.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283
Aaron - Setzen Sie Monica wieder ein
quelle
Nochmals vielen Dank für die Hilfe! Ich kenne diesen Teil aus der Perspektive des Cell Mean Model. Mit dem folgenden Ergebnis von lme () mit zusammengesetzter Symmetrie: Zufällige Effekte: Formel: ~ x - 1 | subj Structure: Compound Symmetry StdDev xx1 1.1913363 xx2 1.1913363 Corr: -0.036 Residual 0.4466733. Ich kann diese Zahlen immer noch nicht mit dem Zellenmittelwertmodell in Einklang bringen. Vielleicht können Sie mir weiter helfen, diese Zahlen zu sortieren?
Bluepole
Auch irgendwelche Gedanken über die Modellformulierung mit anderen Korrelationsstrukturen wie der einfachen symmetrischen Matrix?
Bluepole
Aha! Ich hätte deine Antwort im anderen Thread genauer lesen sollen. Ich habe vorher über die Verwendung von gls () nachgedacht, aber die Korrelationsspezifikation nicht herausgefunden. Es ist interessant, dass lme () mit zusammengesetzter Symmetriestruktur für den Zufallseffekt immer noch den gleichen t-Wert liefert, aber es scheint, dass die Varianzen für die Zufallseffekte nicht direkt interpretierbar sind. Ich schätze deine Hilfe sehr!
Bluepole
Nach einigem Überlegen habe ich das Gefühl, dass meine ursprüngliche Verwirrung immer noch ungelöst ist. Ja, gls kann verwendet werden, um die Korrelationsstruktur und den Mittelwert der quadratischen Rums zu demonstrieren, aber das Modell darunter ist nicht genau dasselbe wie der Paired-t-Test (oder die Einweg-ANOVA mit wiederholten Messungen im Allgemeinen), und eine solche Bewertung ist weiter unterstützt durch die falschen DFs und p-Werte von gls. Im Gegensatz dazu liefert mein Befehl lme mit zusammengesetzter Symmetrie die gleichen Werte für F, DFs und p. Das einzige, worüber ich mich wundere, ist die Parametrisierung des Modells, wie in meinem ursprünglichen Beitrag angegeben. Hilfe da draußen?
Bluepole
Ich weiß nicht, wie ich Ihnen helfen soll. Könntest du aufschreiben, was deiner Meinung nach die beiden verschiedenen Modelle sind? Irgendwas stimmt nicht, wenn Sie an einen von ihnen denken.
Aaron - Wiedereinsetzung von Monica
3

Sie können auch die Verwendung von function mixedin package afexin Betracht ziehen , um p-Werte mit der Kenward-Roger-df-Näherung zurückzugeben, die identische p-Werte als gepaarten t-Test zurückgibt:

library(afex)
mixed(y ~ x + (1|subj), type=3,method="KR",data=myDat) 

Oder

library(lmerTest)
options(contrasts=c('contr.sum', 'contr.poly'))
anova(lmer(y ~ x + (1|subj),data=myDat),ddf="Kenward-Roger")
Tom Wenseleers
quelle