Statistik - Anova

Wissenschaft Physik
1

Hallo!

Ich habe ein Verständnisproblem:

Gegeben sei ein Datensatz von Werten („Response“) mit verschiedenen Behandlungsgruppen („Gruppe“). Die Messungen erfolgte an einer Stichprobe von Individuen („Subjekt“) jeweils vor und nach der Behandlung. Es handelt sich also um ein repeated measures design. Bis dahin ist noch alles gut.

Das ließe sich in R rechnen mit

lm(Response ~ Gruppe + Subjekt) sowie mit
aov(Response ~ Gruppe + Error(Subjekt/Response)) oder auch mit
lme(Response ~ Gruppe, random = ~1|Subjekt)

Alle liefern für den Behandlungseffekt identische Ergebnisse.

Nun kommen Pseudoreplikationen dazu: Jedes Individuum wurde vor und nach der Behandlung mehrmals gemessen („Messung“).

Althergebracht und laienhaft würden die Messwiederholungen zuerst gemittelt und dann würde eine „normale“ repeated measures ANOVA gerechnet.

Nun müßte sich aber doch ein Modell formulieren lassen, welches die Fehlerquadratesummen im Hinblick auf die Messwiederholungen pro Individuum aufspaltet. Ich habe versucht, dazu was zu finden, was ich verstehe - bin aber nur soweit gekommen, als dass sowas wohl als „nested design“ bezeichnet wird. Soweit richtig?

Ich habe das in R probiert mit

lm(Response ~ Gruppe + Subjekt*Messung) bzw. mit
aov(Response ~ Gruppe + Error(Messung/Subjekt/Response)

Ich bin mir nicht sicher, ob das stimmt. Tut es das?
Vorausgesetzt es stimmt, würde ich erwarten, dass die p-Werte am Ende die selben sind, wie sie von der „einfachen“ repeated-measures ANOVA mit zuvor gemittelten Pseudoreplikationen geliefert werden (weil sich der Verlust der Freiheitsgrade durch die Mittelwertbildung in entprechend kleineren Standardfehlern niederschlagen sollte).

In Simulationen mit schön normalverteilten Daten sehe ich aber, dass das nicht der Fall ist. Beide Varianten scheinen die Typ-I-Fehlerrate zu halten, aber die p-Werte sind nicht identisch (nur so ungefähr gleich, manchmal liegen sie auch ein paar Größenordnungen auseinander, insb. bei sehr kleinen Werten, zB. 1E-11 vs. 1E-4).

Also entweder ich verstehe da was grundsätzlich falsch oder ich mache was falsch. Kann mich da jemand helfen?

LG
Jochen

2

Servus Jochen,

endlich mal wieder eine Aufgabe von dir.
also, bevor aiwendil wieder loslegt (und dann damit wieder die übliche Verdächtigen wieder zusammen hocken), geb ich mal meinen Senf dazu. :wink:

Althergebracht und laienhaft würden die Messwiederholungen
zuerst gemittelt und dann würde eine „normale“ repeated
measures ANOVA gerechnet.

Genau.

Nun müßte sich aber doch ein Modell formulieren lassen,
welches die Fehlerquadratesummen im Hinblick auf die
Messwiederholungen pro Individuum aufspaltet. Ich habe
versucht, dazu was zu finden, was ich verstehe - bin aber nur
soweit gekommen, als dass sowas wohl als „nested design“
bezeichnet wird. Soweit richtig?

Nested designs sind der klassiker für Cross-over studien. Nesting bezieht sich darauf einen Faktor in einem anderen „unterzubringen“, im X-over design wäre das z.B. Patient within treatement sequence, um die Varianz durch das nesting besser zu erklären.
Allerdings bezieht sich das nicht auf Messwerte zum selben Zeitpunkt, sondern zu verschiedenen (zumindest habe ich es noch anders gesehen).

Ich habe das in R probiert mit

lm(Response ~ Gruppe + Subjekt*Messung) bzw. mit
aov(Response ~ Gruppe + Error(Messung/Subjekt/Response)

Soweit ich weiß, wird durch / genested, also ist das lm ein nicht genestetes Model, das aov aber schon. Soweit so gut. Allerdings nestest du die Messung im Subjekt (m.E. okay) und dann im repsonse; damit erklärt der response sich zum Teil selber. Das ist so ähnlich wie ein lineares Modell der Art A ~ A+B, da kommt natürlich nicht viel bei raus. Ich würde das aov als
aov(Response ~ Gruppe + Messung/Subjekt + Messung)
aufsetzen, damit auch der repeated-Anteil drin steckt.

Vorausgesetzt es stimmt, würde ich erwarten, dass die p-Werte
am Ende die selben sind, wie sie von der „einfachen“
repeated-measures ANOVA mit zuvor gemittelten
Pseudoreplikationen geliefert werden (weil sich der Verlust
der Freiheitsgrade durch die Mittelwertbildung in entprechend
kleineren Standardfehlern niederschlagen sollte).

Das würde ich nicht erwarten. Freiheitsgrade versus SE muss sich nicht ausgleichen, vor allem wenn die Mittelwertsunterschiede klein sind, denke ich mal, dass der df-Verlust schwerer wiegt.

In Simulationen mit schön normalverteilten Daten sehe ich
aber, dass das nicht der Fall ist. Beide Varianten scheinen
die Typ-I-Fehlerrate zu halten, aber die p-Werte sind nicht
identisch (nur so ungefähr gleich, manchmal liegen sie auch
ein paar Größenordnungen auseinander, insb. bei sehr kleinen
Werten, zB. 1E-11 vs. 1E-4).

Na, dann weißt du doch eigentlich schon, was zu tun ist, gell?

Grüße,
JPL

P.s.nch eine weiter Idee, die funktioniert, wenn du genau zwei Messungen hast und einen Vorwert: durch Subtraktion des Vorwertes von allen Pseudoreplikates kannst du die zeitliche repeated-Struktur auflösen und bist das Problem schon mal los.

3

Hallo JPL,

Mitten in der Nacht noch… wann schläfst Du?
Kennst Du eigentlich Inke und Dirk?

Allerdings bezieht sich das nicht auf Messwerte zum selben
Zeitpunkt, sondern zu verschiedenen (zumindest habe ich es
noch anders gesehen).

Ich habe das in R probiert mit
aov(Response ~ Gruppe + Error(Messung/Subjekt/Response)

Das war falsch. Wie ich jetzt herausgefunden habe, muss es schlicht

aov(Response ~ Gruppe + Error(Subjekt/Gruppe))

sein. Das heißt dann wohl, dass jedes Subjekt innerhalb der Gruppen mehrmals vorkommt. Die Einteilung nach „Messung“ braucht man nicht, weil das Modell ja die Subjekte pro Gruppe zusammenzieht, und das noch nach der Nummer der Messung einzuteilen, ist redundant.

Ich würde das aov als
aov(Response ~ Gruppe + Messung/Subjekt + Messung)
aufsetzen, damit auch der repeated-Anteil drin steckt.

Tja, hab ich probiert, geht aber nicht. Da kommt was ganz anderes raus und zwar tendentiell viiiel zu kleine p’s.

Vorausgesetzt es stimmt, würde ich erwarten, dass die p-Werte
am Ende die selben sind, wie sie von der „einfachen“
repeated-measures ANOVA mit zuvor gemittelten
Pseudoreplikationen geliefert werden (weil sich der Verlust
der Freiheitsgrade durch die Mittelwertbildung in entprechend
kleineren Standardfehlern niederschlagen sollte).

Das würde ich nicht erwarten. Freiheitsgrade versus SE muss
sich nicht ausgleichen, vor allem wenn die
Mittelwertsunterschiede klein sind, denke ich mal, dass der
df-Verlust schwerer wiegt.

Genau so ist es jetzt auch. aov(Response ~ Gruppe + Error(Subjekt/Gruppe)) liefert EXAKT die selben Werte wie aov(meanResponse ~ meanTreatment + Error(meanSubject)), wobei meanResponse die vorher über die Messungen gemittelten Werte enthält.

Na, dann weißt du doch eigentlich schon, was zu tun ist, gell?

Wenn die Ergebnisse nicht der Theorie entsprechen, kann das ein schöner Anlass sein, nach Fehlern in der Theorie zu suchen. Da es mit meinen Statistik-Kenntnissen aber nicht sooo weit her ist, suche ich den Fehler lieber bei mir :wink:

LG & herzlichen Dank!!
Jochen

P.s.nch eine weiter Idee, die funktioniert, wenn du genau zwei
Messungen hast und einen Vorwert: durch Subtraktion des
Vorwertes von allen Pseudoreplikates kannst du die zeitliche
repeated-Struktur auflösen und bist das Problem schon mal los.

Ist das nicht äquivalent zur Bildung der Mittelwerte?

4

Moin Jochen,

Mitten in der Nacht noch… wann schläfst Du?

Es gibt halt interessantere postings als andere :smile:

Kennst Du eigentlich Inke und Dirk?

Kommt drauf an wie die weiter heissen (Repsilber und König?)

Das war falsch. Wie ich jetzt herausgefunden habe, muss es
schlicht

aov(Response ~ Gruppe + Error(Subjekt/Gruppe))

sein. Das heißt dann wohl, dass jedes Subjekt innerhalb der
Gruppen mehrmals vorkommt. Die Einteilung nach „Messung“
braucht man nicht, weil das Modell ja die Subjekte pro Gruppe
zusammenzieht, und das noch nach der Nummer der Messung
einzuteilen, ist redundant.

Du schätzt damit nur den Gruppeneffekt, aber keinen Messungseffekt. Durch den error term wird die Analyse nach den Ausprägungen der Kategorien im error stratifiert, in diesem Fall also Subjekt in Grupe genested?!

z.B. mit
s
Error: s
Df Sum Sq Mean Sq
g 1 3.0476 3.0476

Error: s:g
Df Sum Sq Mean Sq
g 1 1.009 1.009

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
g 1 0.5973 0.5973 1.1962 0.2893
Residuals 17 8.4889 0.4993

Ich würde das aov als
aov(Response ~ Gruppe + Messung/Subjekt + Messung)
aufsetzen, damit auch der repeated-Anteil drin steckt.

Tja, hab ich probiert, geht aber nicht. Da kommt was ganz
anderes raus und zwar tendentiell viiiel zu kleine p’s.

Ist ja klar, ist ja auch ein andere Modell. :smile: aber dass die p’s zu klein sind … hm, wenn du schon weißt, wie groß sie sein müssen, dann brauchst du doch gar nichts mehr rechnen?

Genau so ist es jetzt auch. aov(Response ~ Gruppe +
Error(Subjekt/Gruppe)) liefert EXAKT die selben Werte wie
aov(meanResponse ~ meanTreatment + Error(meanSubject)), wobei
meanResponse die vorher über die Messungen gemittelten Werte
enthält.

d.h. mit
z2
Error: s2
Df Sum Sq Mean Sq
g2 1 1.0159 1.0159

Error: s2:g2
Df Sum Sq Mean Sq
g2 1 0.33634 0.33634

Error: Within
Df Sum Sq Mean Sq F value Pr(>F)
g2 1 0.19911 0.19911 0.3992 0.5724
Residuals 3 1.49630 0.49877

Hm.
Ich würde eher summary(aov(z ~ g + Error(s/m))) verwenden, wobei mir die Info’s zum Error-term auf der R-homepage etwas zu spärlich sind. Hast du da eine bessere Ref?

P.S.noch eine weiter Idee, die funktioniert, wenn du genau zwei
Messungen hast und einen Vorwert: durch Subtraktion des
Vorwertes von allen Pseudoreplikates kannst du die zeitliche
repeated-Struktur auflösen und bist das Problem schon mal los.

Ist das nicht äquivalent zur Bildung der Mittelwerte?

Bei der Mittelwertblidung über die Messungen ersparst du dir den Faktor Messung, bei Differenzbildung zum Vorwert den Faktor Time. aber so wie es aussieht, hast du eh kein repeated model.

Grüße,
JPL

5

Huhu,

Kommt drauf an wie die weiter heissen (Repsilber und König?)

Korrekt. Dachte ich mir doch.

z.B. mit
s

s
Die beiden Analysen liefern auch exakt das gleiche Ergebnis. Auch mußt Du doch die Grouping-Variablen als Faktoren definieren, sonst wird eine multiple lineare Regression gerechnet, die in beiden Fällen natürlich nicht das selbe Ergebnis hat.

Ist ja klar, ist ja auch ein andere Modell. :smile: aber dass die
p’s zu klein sind … hm, wenn du schon weißt, wie groß sie
sein müssen, dann brauchst du doch gar nichts mehr rechnen?

Ja, ja! Ich teste das doch erstmal an simulierten Daten.

müsste mit summary(aov(z2 ~ g2 + Error(s2/g2))) dasselbe wie
oben reauskommen:
…
Hm.

Das Problem waren die nicht als Faktoren angegebenen Grupierungsvariablen.

Ich würde eher summary(aov(z ~ g + Error(s/m))) verwenden,
wobei mir die Info’s zum Error-term auf der R-homepage etwas
zu spärlich sind. Hast du da eine bessere Ref?

Nicht wirklich. http://www.psych.upenn.edu/~baron/rpsych.pdf vielleicht, Kap. 6.10. Hier stehen noch ein paar Happen, aber sehr Lückenhaft: http://zoonek2.free.fr/UNIX/48_R/13.html.

Bin aber auch sehr dankbar für gute Tipps zu guter Literatur!!

LG
Jochen

6

Hi Jochen,

Kommt drauf an wie die weiter heissen (Repsilber und König?)

Korrekt. Dachte ich mir doch.

Interessant - woher kennst du die denn?

Nee, das Beispiel ist falsch. Du hast hier gar keine repeated
measures. Hier das wäre korrekt:

Da hast du natürlich recht: Alles als Covariaten aufzuziehen, bringt natürich nicht wirklich was.

Ist ja klar, ist ja auch ein andere Modell. :smile: aber dass die
p’s zu klein sind … hm, wenn du schon weißt, wie groß sie
sein müssen, dann brauchst du doch gar nichts mehr rechnen?

Ja, ja! Ich teste das doch erstmal an simulierten Daten.

Das war eher die generelle Frage, woher du wissen kannst, wie groß die p’s sein sollen.

Nicht wirklich. http://www.psych.upenn.edu/~baron/rpsych.pdf
vielleicht, Kap. 6.10.

Aber immerhin!
Also: Der Error-term (ET) spezifiziert, wie Streuungen aufgeteilt werden. Alles was nicht im ET auftaucht, geht als Residuum in die Berechnung des F-Wertes von g ein, bzw. für alle im ET spez. Faktoren werden die Residuen berechnet und vom gesamt Residuum abgezogen, bevor der F-Value berechnet wird.
Jetzt musst du dir überlegen, was du der Streuung von g gegenüberstellen willst und alles andere in den ET packen. Aus den 3 Faktoren ergeben sich die Kombinatioen:
s*g*m = s:g:m + s:g + g:m + s:g + s + g + m
Mit summary(aov(z ~ g+Error(s*g*m))) bleibt natürlich nichts übrig.
sinnvoll würde es aber m.E. nach sein, alles drin zu behalten, was g enthält, damit ergäbe sich das das Model
summary(aov(z ~ g+Error(s*m))), was dann übrigens genau das gleiche ist wie summary(aov(z ~ g+s*m)).
Das zeigt, dass der ET eigentlich nur dazu dient, zwar residuen ebenso wie beim konvetionnellen Modell zu berechnen, aber nur die Faktoren zu testen, die nicht im ET vorkommen.

Bin aber auch sehr dankbar für gute Tipps zu guter Literatur!!

Zusammen mit der r-help ergibt sich nun ein recht klares Bild, oder?
Viele Grüße,
JPL

7

Hallo,

Interessant - woher kennst du die denn?

Dirk habe ich in Heidelberg kennengelernt bei einem Kurs zur Auswertung von Microarray-Experimenten, Inke kenne ich durch meinen Chef; ihr Mann war auch mal in hier, bevor er den Ruf nach Lübeck angenommen hat.

Also: Der Error-term (ET) spezifiziert, wie Streuungen
aufgeteilt werden. Alles was nicht im ET auftaucht, geht als
Residuum in die Berechnung des F-Wertes von g ein, bzw. für
alle im ET spez. Faktoren werden die Residuen berechnet und
vom gesamt Residuum abgezogen, bevor der F-Value berechnet
wird.

Ja, genau so stelle ich mir das auch vor.

Jetzt musst du dir überlegen, was du der Streuung von g
gegenüberstellen willst und alles andere in den ET packen. Aus
den 3 Faktoren ergeben sich die Kombinatioen:
s*g*m = s:g:m + s:g + g:m + s:g + s + g + m
Mit summary(aov(z ~ g+Error(s*g*m))) bleibt natürlich nichts
übrig.

Es war aber (z ~ g+Error(s/g)) und nicht (z ~ g+Error(s*g))

sinnvoll würde es aber m.E. nach sein, alles drin zu behalten,
was g enthält, damit ergäbe sich das das Model
summary(aov(z ~ g+Error(s*m))), was dann übrigens genau das
gleiche ist wie summary(aov(z ~ g+s*m)).

Ich verstehe immer noch nicht, warum Du das m drinbehälst. Das gibt dem Modell doch keine Information mehr, die nicht schon in der durch s %in% g gegebenen Gruppenstruktur enthalten ist.

Doch wie dem auch sei, verstehe ich das jetzt so:

z ~ g * s

rechnet alle Varianzen getrennt für g und für s und für die Interaktion g:s. Die Gruppenstruktur ist in g:s enthalten, d.h., die Varianz von g ist der Varianz von g:s gegenüberzustellen, um den Gruppeneffekt zu bekommen.

An deinem Beispiel von vorhin:

der p-Wert war ja 0.2222, gerechnet mit z ~ g + Error(s/g). Nun rechne ich z ~ g*s und erhalte

\> summary(aov(z ~ g \* s))
 Df Sum Sq Mean Sq F value Pr(\>F) 
g 1 2.7222 2.7222 9.8 0.008684 \*\*
s 2 1.7778 0.8889 3.2 0.076945 . 
g:s 2 1.7778 0.8889 3.2 0.076945 . 
Residuals 12 3.3333 0.2778

Der p-Wert für die Gruppe, der hier berechnet wird, ist nicht 0.222, aber es ist auch nicht der gesuchte Effekt! Der sollte nämlich sein:

F = var(g) / var(g:s) = 2.7222/0.8889 = 3.062437

und mit den 1 bzw. 2 d.f. ergibt sich

\> pf(3.062437, 1, 2, lower=F)
[1] 0.2222254

Bingo. (Bei deinen Daten sind die MeanSq für s und g:s gleich; mit anderen Daten nicht unbedingt, und da stimmt die Rechnung auch - ich hab’s probiert :wink:

Zusammen mit der r-help ergibt sich nun ein recht klares Bild,
oder?

Naja, um mir was klarwerden zu lassen braucht’s halt schon was…

Danke nochmal für Deine Hilfe und Diskussionsbereitschaft!

LG
Jochen

8

Hi Jochen,

Dirk habe ich in Heidelberg kennengelernt bei einem Kurs zur
Auswertung von Microarray-Experimenten, Inke kenne ich durch
meinen Chef; ihr Mann war auch mal in hier, bevor er den Ruf
nach Lübeck angenommen hat.

Und wie seit ihr dann auf mich gekommen? Die Welt ist eben ein Dorf.

Jetzt musst du dir überlegen, was du der Streuung von g
gegenüberstellen willst und alles andere in den ET packen. Aus
den 3 Faktoren ergeben sich die Kombinatioen:
s*g*m = s:g:m + s:g + g:m + s:g + s + g + m
Mit summary(aov(z ~ g+Error(s*g*m))) bleibt natürlich nichts
übrig.

Es war aber (z ~ g+Error(s/g)) und nicht (z ~ g+Error(s*g))

Das war auch nur ein Beispiel.
Du kannst natürlich die Streuung für alle möglichen Effekte ausrechnen und abziehen lassen - Ermessenssache.

Ich verstehe immer noch nicht, warum Du das m drinbehälst. Das
gibt dem Modell doch keine Information mehr, die nicht schon
in der durch s %in% g gegebenen Gruppenstruktur enthalten ist.

m beinhaltet ja die Streuung der Messwiederholungen, wenn du es in den Residuen für g behalten willst, musst du das berücksichtigen (was aber weder Modellabhängig ist), wie das zu bewerkstelligen wäre.

Doch wie dem auch sei, verstehe ich das jetzt so:
z ~ g * s
rechnet alle Varianzen getrennt für g und für s und für die
Interaktion g:s. Die Gruppenstruktur ist in g:s enthalten,

aber auch in g:m und g:s:m. Ein Modell, um g versus g:m+g:s+g:m:s zu testen wäre summary(aov(z ~ g + Error(s*m))). Da hier das g nicht im ET vorkommt, werden die ET-parts von den Residuen abgezogen.

d.h., die Varianz von g ist der Varianz von g:s
gegenüberzustellen, um den Gruppeneffekt zu bekommen.
An deinem Beispiel von vorhin:
der p-Wert war ja 0.2222, gerechnet mit z ~ g + Error(s/g).
Nun rechne ich z ~ g*s und erhalte

> summary(aov(z ~ g * s))
Df Sum Sq Mean Sq F value Pr(>F)
g 1 2.7222 2.7222 9.8 0.008684 **
s 2 1.7778 0.8889 3.2 0.076945 .
g:s 2 1.7778 0.8889 3.2 0.076945 .
Residuals 12 3.3333 0.2778

Der p-Wert für die Gruppe, der hier berechnet wird, ist nicht
0.222, aber es ist auch nicht der gesuchte Effekt!

Das liegt auch auf der Hand. Bei aov(z ~ g*s) wird jeder der Faktoren s, g und s:g gegen die Residuals getested, bei summary(aov(z ~ g + Error(s/g))) eben nur g versus s:g.
Das einfachste/sicherste ist eigentlich, sich über summary(aov(z ~ g + Error(s*g*m))) alle Fehler ausgeben zu lassen, und dann selber die Effekte auszurechnen :smile:

Danke nochmal für Deine Hilfe und Diskussionsbereitschaft!

Kein Problem, das ist ja auch mal eine spannende Frage!

Grüße,
JPL