Varianzanalyse für abhängige Stichproben
Erforderliche Pakete laden
library(afex) # ANOVA
library(lsr) # Umwandlung des Datensatzes von wide in long
library(emmeans) # Berechnung der Least–Square-Mittelwerte bzw. Rand-MittelwerteDatensatz einlesen und Variablen spezifizieren
# Beispiel mit einem Within-Faktor und zwei between Faktoren
# Eingabe in der Wide-Version anlog SPSS, d.h. pro Proband eine Zeile.
# Aufbau der Variablennamen der Within-Faktoren: Name.AV_Stufe.wsFaktor1_Stufe.wsFaktor2
Nachhilfe <- rep(c("ja","nein"),c(5,6))
Geschlecht <- rep(c("männlich","weiblich","männlich","weiblich"), c(3,2,3,3))
Note_vorher <- c(5,4,5,4,5,4,4,5,5,4,5)
Note_nachher <- c(4,2,3,4,3,3,4,5,4,5,4)
Note_nach.halbem.Jahr <- c(4,3,4,4,3,3,4,5,4,4,4)
data <- data.frame(Nachhilfe, Geschlecht, Note_vorher, Note_nachher, Note_nach.halbem.Jahr)
# Variablen spezifizieren
Spaltennummern <- c(1,2,3,4,5) # Spaltennummern der für die Analyse benötiten Variablen
between <- c("Nachhilfe","Geschlecht") # Between-Faktoren: gleich wie im Datensatz. Falls keine: between <- NULL
within <- c("Messzeitpunkt") # Within-Faktoren: Abfolge gemäss Variabennamen im Datensatz!
AV <- "Note" # Abhängige Variable: gleich wie im DatensatzDeskriptive Statistik, Signifikanztest, Effektstärke und Teststärke
# Listwise Deletion
data2 <- na.omit(data[, Spaltennummern, drop=FALSE])
# Der Datensatz muss eine Spalte "Block" mit den Probanden-IDs enthalten.
data2$Block <- as.factor(1:nrow(data2))
# Umwandlung des Datensatzes von wide in long
library(lsr)
data.long <- wideToLong(data2, within=within)
# Deskriptive Statistik
Kennwerte <- function(x) c(n=length(x), mean=mean(x), sd=sd(x), se=sd(x)/sqrt(length(x)))
if (is.null(between)) {
Formel <- as.formula(paste(".~", paste(paste(within, collapse="+"))))
} else {
Formel <- as.formula(paste(".~", paste(paste(within, collapse="+"), paste(between, collapse="+"), sep="+")))
}
deskriptive.statistik <- aggregate(Formel, data.long[,c(within, between, AV)], Kennwerte)
# Univariate Tests mit Greenhouse-Geisser und Huynh-Feldt-Korrekturen
# Für Within-Effekte mit df1=1 werden keine GG- und HF-Korrekturen gerechnet!
library(afex)
TAB <- aov_ez(id="Block", dv=AV, data=data.long, between=between, within=within, return="univariate")
# Partielles Eta-Quadrat und Power (univariat, ohne Sphärizitätskorrektur)
Effekte <- length(TAB$ univariate.tests[,1])-1
Power <- vector(length=Effekte)
eta.sq <- vector(length=Effekte)
for (i in 1:Effekte) {
df1 <- TAB$univariate.tests [,2][i+1]
df2 <- TAB$ univariate.tests[,4][i+1]
eta.sq[i] <- TAB$ univariate.tests[,5][i+1]*df1/(TAB$ univariate.tests[,5][i+1]*df1+df2)
ncp<- TAB$ univariate.tests[,5][i+1]*df1 # Analog SPSS: ncp=f^2*df2=(F*df1/df2)*df2; (G*Power: ncp=f^2*N)
Fk<- qf(1-0.05,df1,df2)
Power[i]<- 1-pf(Fk, df1=df1, df2=df2, ncp=ncp)}
# Tabelle mit EtaQuadrat und Power
eta.sq.Power <- data.frame(eta.sq, Power)
rownames(eta.sq.Power) <- rownames(TAB$univariate.tests)[-1]
# Berechnung der Least–Square-Mittelwerte bzw. Rand-Mittelwerte
library(emmeans)
if(is.null(between)){
Faktoren <- within
} else {
Faktoren <- as.formula(paste("~", paste(between, collapse="*"), "*", paste(within, collapse="*")))}
TAB2 <- aov_ez(id="Block", dv=AV, data=data.long, between=between, within=within)
ls.means <- lsmeans(TAB2, Faktoren)
list(Deskriptive.Statistik = deskriptive.statistik, ANOVA = TAB, Eta.Sq.Power=eta.sq.Power, "Geschätzte Randmittelwerte"=ls.means)## $Deskriptive.Statistik
## Messzeitpunkt Nachhilfe Geschlecht Note.n Note.mean Note.sd
## 1 nach.halbem.Jahr ja männlich 3.0000000 3.6666667 0.5773503
## 2 nachher ja männlich 3.0000000 3.0000000 1.0000000
## 3 vorher ja männlich 3.0000000 4.6666667 0.5773503
## 4 nach.halbem.Jahr nein männlich 3.0000000 4.0000000 1.0000000
## 5 nachher nein männlich 3.0000000 4.0000000 1.0000000
## 6 vorher nein männlich 3.0000000 4.3333333 0.5773503
## 7 nach.halbem.Jahr ja weiblich 2.0000000 3.5000000 0.7071068
## 8 nachher ja weiblich 2.0000000 3.5000000 0.7071068
## 9 vorher ja weiblich 2.0000000 4.5000000 0.7071068
## 10 nach.halbem.Jahr nein weiblich 3.0000000 4.0000000 0.0000000
## 11 nachher nein weiblich 3.0000000 4.3333333 0.5773503
## 12 vorher nein weiblich 3.0000000 4.6666667 0.5773503
## Note.se
## 1 0.3333333
## 2 0.5773503
## 3 0.3333333
## 4 0.5773503
## 5 0.5773503
## 6 0.3333333
## 7 0.5000000
## 8 0.5000000
## 9 0.5000000
## 10 0.0000000
## 11 0.3333333
## 12 0.3333333
##
## $ANOVA
##
## Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
##
## Sum Sq num Df Error SS den Df F value
## (Intercept) 515.56 1 7.2778 7 495.8838
## Nachhilfe 1.39 1 7.2778 7 1.3359
## Geschlecht 0.15 1 7.2778 7 0.1484
## Nachhilfe:Geschlecht 0.06 1 7.2778 7 0.0534
## Messzeitpunkt 4.49 2 3.5556 14 8.8472
## Nachhilfe:Messzeitpunkt 1.33 2 3.5556 14 2.6250
## Geschlecht:Messzeitpunkt 0.35 2 3.5556 14 0.6806
## Nachhilfe:Geschlecht:Messzeitpunkt 0.15 2 3.5556 14 0.2917
## Pr(>F)
## (Intercept) 9.311e-08 ***
## Nachhilfe 0.285683
## Geschlecht 0.711479
## Nachhilfe:Geschlecht 0.823802
## Messzeitpunkt 0.003281 **
## Nachhilfe:Messzeitpunkt 0.107617
## Geschlecht:Messzeitpunkt 0.522321
## Nachhilfe:Geschlecht:Messzeitpunkt 0.751447
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Mauchly Tests for Sphericity
##
## Test statistic p-value
## Messzeitpunkt 0.42187 0.075085
## Nachhilfe:Messzeitpunkt 0.42187 0.075085
## Geschlecht:Messzeitpunkt 0.42187 0.075085
## Nachhilfe:Geschlecht:Messzeitpunkt 0.42187 0.075085
##
##
## Greenhouse-Geisser and Huynh-Feldt Corrections
## for Departure from Sphericity
##
## GG eps Pr(>F[GG])
## Messzeitpunkt 0.63366 0.01255 *
## Nachhilfe:Messzeitpunkt 0.63366 0.13721
## Geschlecht:Messzeitpunkt 0.63366 0.46567
## Nachhilfe:Geschlecht:Messzeitpunkt 0.63366 0.65590
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## HF eps Pr(>F[HF])
## Messzeitpunkt 0.7098446 0.009467578
## Nachhilfe:Messzeitpunkt 0.7098446 0.130558517
## Geschlecht:Messzeitpunkt 0.7098446 0.479684744
## Nachhilfe:Geschlecht:Messzeitpunkt 0.7098446 0.679911971
##
## $Eta.Sq.Power
## eta.sq Power
## Nachhilfe 0.160256410 0.17084470
## Geschlecht 0.020764120 0.06299155
## Nachhilfe:Geschlecht 0.007575758 0.05465780
## Messzeitpunkt 0.558282209 0.93004974
## Nachhilfe:Messzeitpunkt 0.272727273 0.43648194
## Geschlecht:Messzeitpunkt 0.088607595 0.14198931
## Nachhilfe:Geschlecht:Messzeitpunkt 0.040000000 0.08741333
##
## $`Geschätzte Randmittelwerte`
## Nachhilfe Geschlecht Messzeitpunkt lsmean SE df lower.CL upper.CL
## ja männlich nach.halbem.Jahr 3.68 0.422 13.5 2.77 4.59
## nein männlich nach.halbem.Jahr 4.02 0.422 13.5 3.11 4.92
## ja weiblich nach.halbem.Jahr 3.52 0.485 14.8 2.48 4.55
## nein weiblich nach.halbem.Jahr 4.02 0.422 13.5 3.11 4.92
## ja männlich nachher 3.02 0.422 13.5 2.11 3.92
## nein männlich nachher 4.02 0.422 13.5 3.11 4.92
## ja weiblich nachher 3.52 0.485 14.8 2.48 4.55
## nein weiblich nachher 4.35 0.422 13.5 3.44 5.26
## ja männlich vorher 4.68 0.422 13.5 3.77 5.59
## nein männlich vorher 4.35 0.422 13.5 3.44 5.26
## ja weiblich vorher 4.52 0.485 14.8 3.48 5.55
## nein weiblich vorher 4.68 0.422 13.5 3.77 5.59
##
## Warning: EMMs are biased unless design is perfectly balanced
## Confidence level used: 0.95Multivariate Tests (Phillai)
library(afex)
TAB <- aov_ez(id="Block", dv=AV, data=data.long, between=between, within=within, return="Anova")
TAB##
## Type III Repeated Measures MANOVA Tests: Pillai test statistic
## Df test stat approx F num Df den Df
## (Intercept) 1 0.98608 495.88 1 7
## Nachhilfe 1 0.16026 1.34 1 7
## Geschlecht 1 0.02076 0.15 1 7
## Nachhilfe:Geschlecht 1 0.00758 0.05 1 7
## Messzeitpunkt 1 0.64783 5.52 2 6
## Nachhilfe:Messzeitpunkt 1 0.35714 1.67 2 6
## Geschlecht:Messzeitpunkt 1 0.39552 1.96 2 6
## Nachhilfe:Geschlecht:Messzeitpunkt 1 0.05814 0.19 2 6
## Pr(>F)
## (Intercept) 9.311e-08 ***
## Nachhilfe 0.28568
## Geschlecht 0.71148
## Nachhilfe:Geschlecht 0.82380
## Messzeitpunkt 0.04368 *
## Nachhilfe:Messzeitpunkt 0.26567
## Geschlecht:Messzeitpunkt 0.22087
## Nachhilfe:Geschlecht:Messzeitpunkt 0.83553
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Geplante Vergleiche
Die unten verwendete Funktion contrast() berechnet die Prüfgrösse t wie folgt:
# Means <- summary(M)$lsmean; Contrast <- unlist(C); t <- sum(Contrast*Means)/sqrt(QSe/(n*dfe)*sum(Contrast^2))
# Means = Geschätzte Randmittelwerte, Contrast = Kontraskoefizienten
# QSe = Fehlerquadratsumme, dfe= Fehlerfreiheitsgrade, n = Zellgrösse
# QSe und dfe können aus der ANOVA-Tabelle beim entsprechnednen Effekt abgelesen werden.
# Falls die Sphärizitätsbedingung verletzt ist sollte man den p-Werte für within-Effekte
# mit korrigierten Freiheitsgraden berechnen: pf(t^2, GG.eps, dfe*GG.eps, lower.tail=FALSE)
Beispiel 1: Haupteffektkontrast Messzeitpunkt (vorher vs. nachher und nach.halbem.Jahr)
Die Abfolge der Kontrastkoeffizienten muss derjenigen der Stufen entsprechen.
C <- list(c(-1,-1,2)) # Messzeitpunkt (vorher vs. nachher und nach.halbem.Jahr)
M <- lsmeans(TAB2, ~Messzeitpunkt) # Least–Square-Mittelwerte## NOTE: Results may be misleading due to involvement in interactionscontrast(M, C)## contrast estimate SE df t.ratio p.value
## c(-1, -1, 2) 1.58 0.378 14 4.189 0.0009
##
## Results are averaged over the levels of: Nachhilfe, Geschlecht# p-Wert mit "GG eps" korrigierten Freiheitsgraden:
pf(4.189^2, 0.634, 14*0.634, lower.tail=FALSE)## [1] 0.004640286Beispiel 2: Interaktionskontrast
Messzeitpunkt (vorher vs. nachher und nach.halbem.Jahr) * Nachhilfe (ja vs. nein) * Geschlecht (männlich vs. weiblich)
Messzeitpunkt <- c(-1,-1,2) # Messzeitpunkt (vorher vs. nachher und nach.halbem.Jahr)
Nachhilfe <- c(-1,1) # Nachhilfe (ja vs. nein)
Geschlecht <- c(-1,1) # Geschlecht (männlich vs. weiblich)
C <- list(as.numeric(Messzeitpunkt %o% Nachhilfe %o% Geschlecht)) # Interaktion
M <- lsmeans(TAB2, ~Messzeitpunkt+Nachhilfe+Geschlecht) # Least–Square-Mittelwerte
contrast(M, C)## contrast estimate SE df t.ratio p.value
## c(-1, -1, 2, 1, 1, -2, 1, 1, -2, -1, -1, 2) 1 1.51 14 0.661 0.5191# p-Wert mit "GG eps" korrigierten Freiheitsgraden:
pf(0.661^2, 0.634, 14*0.634, lower.tail=FALSE)## [1] 0.4373113Beispiel 3: Interaktionskontrast
Messzeitpunkt (vorher vs. nachher und nach.halbem.Jahr) * Nachhilfe (ja vs. nein) * Geschlecht (männlich vs. weiblich)
Messzeitpunkt <- c(-1,-1,2) # Messzeitpunkt (vorher vs. nachher und nach.halbem.Jahr)
Nachhilfe <- c(-1,1) # Nachhilfe (ja vs. nein)
C <- list(as.numeric(Messzeitpunkt %o% Nachhilfe)) # Interaktion
M <- lsmeans(TAB2, ~Messzeitpunkt+Nachhilfe) # Least–Square-Mittelwerte## NOTE: Results may be misleading due to involvement in interactionscontrast(M, C)## contrast estimate SE df t.ratio p.value
## c(1, 1, -2, -1, -1, 2) -1.5 0.756 14 -1.984 0.0672
##
## Results are averaged over the levels of: Geschlecht# p-Wert mit "GG eps" korrigierten Freiheitsgraden:
pf(1.984^2, 0.634, 14*0.634, lower.tail=FALSE)## [1] 0.08765818Post-hoc-Vergleiche nach Tukey
# M <- lsmeans(TAB2, ~Messzeitpunkt+Nachhilfe+Geschlecht)
# contrast(M, method="pairwise")Grafiken
Für Boxplots, Histogramme und Liniendiagramme (Interaktionsplots) müssen Sie den Datensatz “data.long” verwenden!
Speichern Sie Ihren Originaldatensatz unter data4: data4 <- data
Setzen Sie data <- data.long
Spezifizieren Sie die Variablen gemäss Vorlage des gewünschten Plots.
Beispiel für den Interaktionsplot:
data4 <- data # Originaldatensatz unter data4 speichern
data <- data.long
# Variablen spezifizieren
Faktor1 <- "Messzeitpunkt" # Name des ersten Faktors eingeben (Abszisse)
Faktor2 <- "Nachhilfe" # Name des zweiten Faktors eingeben (Farbe)
Faktor3 <- "Geschlecht" # Unterteilung der Plots
AV <- "Note" # Name der abhängigen Variable eingeben
CI <- .95 # gewünschtes Konfidenzintervall eingeben. Kein Konfidenzintervall: CI <- 0
# Anschliessend alle Befehle nach 'Variablen spezifizieren' ausführen.Hinweis: Interaktionsplot mit den geschätzten Randmittelwerten
# Faktornamen spezifizieren
Faktor1 <- "Messzeitpunkt" # Name des ersten Faktors eingeben (Abszisse)
Faktor2 <- "Nachhilfe" # Name des zweiten Faktors eingeben (Farbe)
Faktor3 <- "Geschlecht" # Unterteilung der Plots
AV <- "Note" # Name der abhängigen Variable eingeben
CI <- .95 # gewünschtes Konfidenzintervall eingeben. Kein Konfidenzintervall: CI <- 0
data.aggr <- summary(ls.means) # Datensatz mit Mittelwerten und KI generieren
# "Liniendiagramm für Mittelwerte" aufrufen
# Den Code für die "Schwarzweissversion" oder die "Farbige Version" des entsprechenden Designs laufen lassen.
# Es werden die Mittelwerte mit 95%-Konfidenzintervallen dargestellt.