/* Reproduktion der Analysen zum Skript "Anova.pdf" mit Stata */ * ============================================================================== /* Die Ausgabe soll als "Text" in eine Datei mit dem Namen "bsp_anova.log" geschrieben werden. Dazu wird zunächst ein eventuell aktiver Logfile geschlossen. Anschliessend wird eine neuer Logfile mit dem entsprechen- den Namen spezifiziert. ACHTUNG: Stata legt die Datei im aktuellen Laufwerk und Pfad an. Falls dabei eine Fehlermeldung ausgegeben wird, kann es sein, dass für den Pfad keine Schreibberechtigung existiert. In diesem Fall muss entweder der Dateiname "Loesung.log" noch durch eine gültige Laufwerks- und Pfad- angabe ergänzt werden, z.B. "K:/bsp_anova.log" oder so ähnlich, oder der -log- Anweisung muss eine -cd- Anweisung vorangestellt werden, mit der Stata ein gültiger Pfad zum Lesen und Schreiben von Dateien zuge- wiesen wird (siehe -h cd-)! */ cap log close log using "bsp_anova.log", text replace /* Zuvor werden einige Pakete installiert, sofern sie noch nicht installiert worden sind; die nächsten 8 Anweisungen setzen Internetanschluss voraus: */ cap which tukeyhsd if _rc net install tukeyhsd, from(http://www.ats.ucla.edu/stat/stata/ado/analysis) cap which tkcomp if _rc net install tkcomp, from(http://www.ats.ucla.edu/stat/stata/ado/analysis) cap which anovaplot if _rc net install gr0009_1, from(http://www.stata-journal.com/software/sj10-1) cap which wsanova if _rc net install sg103, from(http://www.stata.com/stb/stb47) * ---------------------------------------------------------------------- * 1. Reproduktion der Analysen aus "anova1.r" * Beispiel 1: /* Hier werden auch gleich die Daten für Beispiel 2 (y2) eingegeben, was Schreibarbeit spart */ clear input gruppe y1 y2 1 16 16 1 16 16 1 14 14 1 13 13 1 12 12 2 16 23 2 14 21 2 13 20 2 13 20 2 10 17 3 16 26 3 14 24 3 12 22 3 10 20 3 10 20 4 16 11 4 14 9 4 12 7 4 12 7 4 12 7 5 14 10 5 13 9 5 13 9 5 10 6 5 10 6 end label define grp 1 "Gruppe 1" /// 2 "Gruppe 2" /// 3 "Gruppe 3" /// 4 "Gruppe 4" /// 5 "Gruppe 5" label values gruppe grp * Mittelwerte für alle Teilgruppen: tabstat y1, statistics( mean sd count ) by(gruppe) columns(statistics) f(%4.1f) /* Grafik der individuellen Werte für alle Teilgruppen: Um die Mittelwertslinien in die Grafik einzuzeichnen, werden zunächst "quietly" (d.h. ohne Output) die Mittelwerte der Gesamtgruppe sowie der Teilgruppen berechnet und die unmittelbar nach der Berechnung zur Verfügung stehenden Werte (siehe -h return-) in sogenannte lokale Makros kopiert. Man beachte, dass ein Makro in Stata ein Skalar darstellt, dem mittels -local- (oder -global-) Werte zugewiesen werden können. Lokale Makros können mit -local- definiert werden. Sie können anschließend benutzt werden, indem man ihre Namen zwischen den Begrenzungszeichen ` und ' spe- zifiziert (das erste Zeichen ist der Accent grave, auf der deutschen Tastatur mit der Kombination gefolgt von einem Leerzeichen erzeugbar; das zweite Zeichen ist das einfache Hochkomma, auf der deutschen Tastatur mit der Kombination er- zeugbar). */ qui sum y1 local meantotal `r(mean)' qui sum y1 if gruppe==1 local meang1 `r(mean)' qui sum y1 if gruppe==2 local meang2 `r(mean)' qui sum y1 if gruppe==3 local meang3 `r(mean)' qui sum y1 if gruppe==4 local meang4 `r(mean)' qui sum y1 if gruppe==5 local meang5 `r(mean)' twoway (scatter y1 gruppe, msymbol(circle_hollow) jitter(1)) /// (scatteri `meang1' 0.8 `meang1' 1.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang2' 1.8 `meang2' 2.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang3' 2.8 `meang3' 3.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang4' 3.8 `meang4' 4.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang5' 4.8 `meang5' 5.2, recast(line) lcolor(blue) lpattern(dash)) /// , title("Variabilität 'within' und 'between' (Beispiel 1)", size(medium)) /// yline(`meantotal') ytitle(y1) yscale(range(5 26)) ylabel(5(5)25) /// xtitle(Gruppe) legend(off) name(Beispiel_1, replace) /* Die Varianzanalyse "zu Fuß" lässt sich mit Stata idealerweise mit Hilfe von Mata (einer Matrixprogrammierumgebung, die innerhalb des Programms lauffähig ist) programmieren, was aber ein eigenes Kapitel ist und deshalb hier weggelassen wird. */ * Einfaktorielle Varianzanalyse mit ANOVA: oneway y1 gruppe /* Post-hoc-Tests sind hier nicht nötig, da der globale Test nicht signifikant geworden ist. */ * -------------------------------------------------------------------- * Beispiel 2 (die Daten - y2 - wurden oben schon eingegeben): * Mittelwerte für alle Teilgruppen: tabstat y2, statistics( mean sd count ) by(gruppe) columns(statistics) f(%4.1f) * Grafik der individuellen Werte für alle Teilgruppen: qui sum y2 local meantotal `r(mean)' qui sum y2 if gruppe==1 local meang1 `r(mean)' qui sum y2 if gruppe==2 local meang2 `r(mean)' qui sum y2 if gruppe==3 local meang3 `r(mean)' qui sum y2 if gruppe==4 local meang4 `r(mean)' qui sum y2 if gruppe==5 local meang5 `r(mean)' twoway (scatter y2 gruppe, msymbol(circle_hollow) jitter(1)) /// (scatteri `meang1' 0.8 `meang1' 1.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang2' 1.8 `meang2' 2.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang3' 2.8 `meang3' 3.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang4' 3.8 `meang4' 4.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang5' 4.8 `meang5' 5.2, recast(line) lcolor(blue) lpattern(dash)) /// , title("Variabilität 'within' und 'between' (Beispiel 2)", size(medium)) /// yline(`meantotal') ytitle(y2) xtitle(Gruppe) legend(off) /// name(Beispiel_2, replace) /* Einfaktorielle Varianzanalyse mit -oneway- (als multiple post-hoc-Tests stehen in Stata der Bonferroni-, der Sidák- und der Scheffé-Test zur Verfügung. Hier wird der Scheffé-Test gewählt, der bei paarweisen Vergleichen aller Mittelwerte nicht so konservativ wie die anderen beiden Tests ist: */ oneway y2 gruppe, scheffe /* Als optimaler gilt der Holm-Test, dessen Spezifikation allerdings aufwändiger ist. Er setzt eine Varianzanalyse mit -anova- voraus und wird mit dem "post- estimation"-Befehl -test- berechnet. -test- benötigt die Angabe der Kontraste (jede der 5 Gruppen gegen jede). Dazu wird eine Matrix der Kontraste spezifi- ziert, in der pro Zeile ein Kontrast (1 vs. -1) aufgeführt wird - die Zeilen- summe muss 0 sein und die letzte Spalte steht für die Konstante. Möglich sind auch der Tukey-HSD und der Tukey-Kramer-Test, letzterer ist bei stark ungleicher Zellbesetzung dem Tukey-HSD-Test vorzuziehen, hier ist dies wegen gleicher Zell- besetzung nicht nötig - beide Tests liefern deshalb hier das gleiche Ergebnis. */ anova y2 gruppe matrix input kontraste=(1,-1,0,0,0,0\1,0,-1,0,0,0\1,0,0,-1,0,0\1,0,0,0,-1,0\ /// 0,1,-1,0,0,0\0,1,0,-1,0,0\0,1,0,0,-1,0\ /// 0,0,1,-1,0,0\0,0,1,0,-1,0\ /// 0,0,0,1,-1,0) matrix list kontraste test, test(kontraste) mtest(holm) /* Holm-Test */ tukeyhsd gruppe /* Tukey HSD-Test */ tkcomp gruppe /* Tukey-Kramer-Test */ * ==================================================================== * 2. Reproduktion der Analysen aus "anova2.r" clear input male caffeine y 1 1 6 1 1 15 1 1 12 1 1 12 1 1 13 1 2 12 1 2 10 1 2 12 1 2 13 1 2 7 1 3 10 1 3 13 1 3 15 1 3 12 1 3 10 1 4 9 1 4 10 1 4 7 1 4 12 1 4 7 0 1 10 0 1 13 0 1 4 0 1 9 0 1 5 0 2 9 0 2 7 0 2 10 0 2 7 0 2 13 0 3 12 0 3 13 0 3 15 0 3 10 0 3 13 0 4 4 0 4 7 0 4 6 0 4 9 0 4 9 end label define male 0 "female" 1 "male" label define caffeine 1 "large" /// 2 "moderate" /// 3 "small" /// 4 "zero" label values male male label values caffeine caffeine label variable male "Gender" label variable caffeine "dose" label variable y "Testscore" * Deskriptive Statistiken für "male": tabstat y, statistics( mean sd count ) by(male) columns(statistics) f(%4.1f) * Deskriptive Statistiken für "caffeine": tabstat y, statistics( mean sd count ) by(caffeine) columns(statistics) f(%4.1f) * Mittelwerte für die einzelnen Zellen: tabulate male caffeine, summarize(y) nostandard nofreq * ---------------------------------------------------------------------- * Grafik der y-Werte der Gruppen: * a) Mittelwerte nach Geschlecht: qui sum y local meantotal `r(mean)' qui sum y if male==0 local meang1 `r(mean)' qui sum y if male==1 local meang2 `r(mean)' twoway (scatter y male, msymbol(circle_hollow) jitter(1)) /// (scatteri `meang1' -0.2 `meang1' 0.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang2' 0.8 `meang2' 1.2, recast(line) lcolor(blue) lpattern(dash)) /// , title("Werte und Mittelwerte nach Geschlecht", size(medium)) /// yline(`meantotal') ytitle(testscore) /// xlabel( 0(1)1, valuelabel) xtitle(sex) /// legend(off) name(Caffeine_1, replace) * b) Mittelwerte nach Koffein-Bedingung: qui sum y local meantotal `r(mean)' qui sum y if caffeine==1 local meang1 `r(mean)' qui sum y if caffeine==2 local meang2 `r(mean)' qui sum y if caffeine==3 local meang3 `r(mean)' qui sum y if caffeine==4 local meang4 `r(mean)' twoway (scatter y caffeine, msymbol(circle_hollow) jitter(1)) /// (scatteri `meang1' 0.8 `meang1' 1.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang2' 1.8 `meang2' 2.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang3' 2.8 `meang3' 3.2, recast(line) lcolor(blue) lpattern(dash)) /// (scatteri `meang4' 3.8 `meang4' 4.2, recast(line) lcolor(blue) lpattern(dash)) /// , title("Werte und Mittelwerte nach Koffein-Bedingung", size(medium)) /// yline(`meantotal') ytitle(testscore) /// xlabel( 1(1)4, valuelabel) xtitle("caffeine dose") /// legend(off) name(Caffeine_2, replace) * ---------------------------------------------------------------------- * Zweifaktorielle Varianzanalyse: /* Ab Version 11.0 müssen Haupt- und Interaktionseffekte nicht mehr separat spezifiziert werden, sondern können mit Hilfe von Doppelkreuzen simultan spezifiziert werden, z.B. als -anova y male##caffeine-, see -h fvvarlist-. Damit die Spezifikation, wie sie in früheren Stata-Versionen gültig war, auch unter Version 11.0 läuft, wird der -anova- Anweisung -version 9:- vorangestellt, wodurch auch mit einer Stata Version 11.0 die Anweisung wie in Version 9.0 ausgeführt wird (siehe -h version-): */ version 9: anova y male caffeine male*caffeine * Grafik: anovaplot caffeine male, scatter(msym(none)) title("Interaktionsplot") /* Es geht auch komplizierter mittels aggregierten Daten (via -collapse-): Mittels -preserve- werden die Daten "eingefroren", um dann später mittels -restore- wieder "aufgetaut" zu werden, so dass Datentransformationen dazwischen nur temporär sind und die Daten nicht wirklich verändern. */ preserve collapse (mean) testscore=y, by(caffeine male) separate testscore, by(male) twoway (line testscore0 testscore1 caffeine), /// xlabel( 1(1)4, valuelabel) xtitle(Dose) ytitle(Testscore) /// legend(order(1 "female" 2 "male")) title("Interaktionsplot") /// name(Caffeine_3, replace) restore * ====================================================================== * 3. Reproduktion der Analysen aus "anova2mw.r" /* Die Daten werden im folgenden aus dem Netz geladen, falls kein Internet- anschluss besteht, muss der ensprechende Pfad zur Datei auf dem Rechner spezifiziert werden: */ clear insheet using "http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Lehre/StatIIKrim/Uebungen/anova2mw.dat" label define grp 0 "Kontrollgruppe" 1 "Experimentalgruppe" label values gruppe grp label variable gruppe "treatment" label variable y1 "y t1" label variable y2 "y t2" * Statistische Kennwerte: tabstat y1, statistics( mean sd count ) by(gruppe) columns(statistics) tabstat y2, statistics( mean sd count ) by(gruppe) columns(statistics) * Für eine Varianzanalyse mit Messwiederholung werden die Daten zunächst * in "langes" Format umgewandelt, hierfür ist eine Variable nötig, die * die Fälle eindeutig kennzeichnet (id): gen id = _n * Daten im "weiten" Format: list * Umwandlung in "langes" Format: reshape long y, i(id) j(zeitpunkt) * Daten im "langen" Format: list * Varianzanalyse mit Messwiederholung: version 9: anova y gruppe / id|gruppe zeitpunkt zeitpunkt*gruppe, repeated(zeitpunkt) /* Die Benutzung der Prozedur -wsanova- ("within subjects ANOVA") ist wesentlich intuitiver: */ wsanova y zeitpunkt, id(id) between(gruppe) /* Für diese und weitere varianzanalytische Modelle mit Messwiederholung finden sich Erläuterungen und Analysebeispiele unter http://www.stata.com/support/faqs/stat/anova2.html */ * -------------------------------------------------------------------- * Grafische Darstellung der Veränderungen: * Hier kann -anovaplot- nicht verwendet werden, deshalb also anders: preserve collapse (mean) testscore=y, by(gruppe zeitpunkt) separate testscore, by(gruppe) replace testscore0 = testscore0/12 replace testscore1 = testscore1/12 twoway (line testscore0 testscore1 zeitpunkt), /// xlabel(1 "t1" 2 "t2") xtitle(Zeitpunkt) ytitle(Testscore) /// legend(order(1 "Kontrollgruppe" 2 "Experimentalgruppe")) /// title("Interaktionsplot") name(ANOVA_MW, replace) restore * -------------------------------------------------------------------- /* Demonstration, dass die Interaktion des Messwiederholungsfaktors (hier "zeitpunkt") mit dem "between"-Faktor (hier "gruppe") für den Fall von nur zwei Messzeitpunkten auch mit Hilfe einer Varianzanalyse der Mittel- wertsdifferenzen geprüft werden kann Hierzu müssen die Daten zunächst wieder mit -reshape- in das "wide"- Format zurückverwandelt werden (die -list- Anweisung zeigt, wie die Daten im "long"- und im "wide"-Format aussehen): */ list reshape wide y, i(id) j(zeitpunkt) list * Berechnung der Veränderungswerte dy: gen dy = y2-y1 * Varianzanalyse der Veränderungswerte: oneway dy gruppe * ====================================================================== * Schließen des Logfiles mit den Ergebnissen der Analysen: log close