Erforderliche Pakete laden
library(car) # MANOVA
library(afex) # ANOVA
library(heplots) # HE Plot
library(candisc) # Diskriminazfunktionsplots
Datensatz einlesen und Variablen spezifizieren
# Datensatz einlesen
Geschlecht <- as.factor(rep(c("männlich","weiblich"), each = 4, times = 3))
Medikament <- as.factor(rep(c("a1","a2","a3"), each=8))
x1 <- c(5,5,9,7,7,6,9,8,7,7,9,6,10,8,7,6,21,14,17,12,16,14,14,10)
x2 <- c(6,4,9,6,10,6,7,10,6,7,12,8,13,7,6,9,15,11,12,10,12,9,8,5)
data <- data.frame(Geschlecht, Medikament, x1, x2)
# Variablen spezifizieren
Faktoren <- c("Geschlecht", "Medikament") # Namen der Faktoren mit unabhängigen Stichproben eingeben
AV <- c("x1", "x2") # Namen der abhängigen Variablen eingeben
Deskriptive Statistik, Signifikanztest
# Listwise Deletion (Zeilen mit Missings in den berücksichtigten Variablen eliminieren)
data2 <- na.omit(data[,c(AV, Faktoren)])
# Der Datensatz muss eine Spalte "Block" mit den Probanden-IDs enthalten.
data2$Block <- as.factor(1:nrow(data2))
# Deskriptive Statistik
# Zu berechnende Kennwerte spezifizieren.
Kennwerte <- function(x) c(mean = mean(x), sd = sd(x), se=sd(x)/sqrt(length(x)))
# Tabelle
Formel <- as.formula(paste(".~", paste(Faktoren, collapse="+")))
tab <- aggregate(Formel, data2[,c(Faktoren, AV)], Kennwerte)
n <- aggregate(Formel, data2[,c(Faktoren, AV[1])], length)
deskriptive.statistik <- data.frame(tab,n=n[,ncol(n)])
# MANOVA: Quadratsummenzerlegung nach Modell III
library(car)
RS <- paste("C(", Faktoren[1], ",contr.helmert)", sep="")
if (length(Faktoren)>1){
for (i in 2:length(Faktoren)) RS <- paste(RS, "*C(", Faktoren[i], ",contr.helmert)", sep="")
}
model <- as.formula(paste("cbind(",paste(AV,collapse=","),")~",RS, sep=""))
fit <- lm(model, data=data2)
man.res <- Manova(fit, type="III")
TAB.MANOVA <- summary(man.res, multivariate=TRUE)
# Effektstärken der MANOVA
library(heplots)
eta.sq <- etasq(man.res)
# ANOVAs: Quadratsummenzerlegung nach Modell III
library(afex)
TAB.ANOVA <- list()
for (i in 1:length(AV)){
TAB.ANOVA[[i]] <- aov_ez(id="Block", dv=AV[i], data=data2, between=Faktoren, return="Anova")
}
names(TAB.ANOVA) <- AV
list("Deskriptive Statistik" = deskriptive.statistik, MANOVA=TAB.MANOVA, Effektstärke=eta.sq, ANOVAs=TAB.ANOVA)
## $`Deskriptive Statistik`
## Geschlecht Medikament x1.mean x1.sd x1.se x2.mean x2.sd
## 1 männlich a1 6.5000000 1.9148542 0.9574271 6.250000 2.061553
## 2 weiblich a1 7.5000000 1.2909944 0.6454972 8.250000 2.061553
## 3 männlich a2 7.2500000 1.2583057 0.6291529 8.250000 2.629956
## 4 weiblich a2 7.7500000 1.7078251 0.8539126 8.750000 3.095696
## 5 männlich a3 16.0000000 3.9157800 1.9578900 12.000000 2.160247
## 6 weiblich a3 13.5000000 2.5166115 1.2583057 8.500000 2.886751
## x2.se n
## 1 1.030776 4
## 2 1.030776 4
## 3 1.314978 4
## 4 1.547848 4
## 5 1.080123 4
## 6 1.443376 4
##
## $MANOVA
##
## Type III MANOVA Tests:
##
## Sum of squares and products for error:
## x1 x2
## x1 94.5 76.5
## x2 76.5 114.0
##
## ------------------------------------------
##
## Term: (Intercept)
##
## Sum of squares and products for the hypothesis:
## x1 x2
## x1 2281.5 2028.000
## x2 2028.0 1802.667
##
## Multivariate Tests: (Intercept)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.960659 207.5601 2 17 1.1381e-12 ***
## Wilks 1 0.039341 207.5601 2 17 1.1381e-12 ***
## Hotelling-Lawley 1 24.418839 207.5601 2 17 1.1381e-12 ***
## Roy 1 24.418839 207.5601 2 17 1.1381e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: C(Geschlecht, contr.helmert)
##
## Sum of squares and products for the hypothesis:
## x1 x2
## x1 0.6666667 0.6666667
## x2 0.6666667 0.6666667
##
## Multivariate Tests: C(Geschlecht, contr.helmert)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 1 0.0074631 0.06391302 2 17 0.93831
## Wilks 1 0.9925369 0.06391302 2 17 0.93831
## Hotelling-Lawley 1 0.0075192 0.06391302 2 17 0.93831
## Roy 1 0.0075192 0.06391302 2 17 0.93831
##
## ------------------------------------------
##
## Term: C(Medikament, contr.helmert)
##
## Sum of squares and products for the hypothesis:
## x1 x2
## x1 301.0 97.50000
## x2 97.5 36.33333
##
## Multivariate Tests: C(Medikament, contr.helmert)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 2 0.880378 7.07686 4 36 0.00026018 ***
## Wilks 2 0.168630 12.19913 4 34 2.9586e-06 ***
## Hotelling-Lawley 2 4.639537 18.55815 4 32 5.5968e-08 ***
## Roy 2 4.576027 41.18424 2 18 1.9190e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## ------------------------------------------
##
## Term: C(Geschlecht, contr.helmert):C(Medikament, contr.helmert)
##
## Sum of squares and products for the hypothesis:
## x1 x2
## x1 14.33333 21.33333
## x2 21.33333 32.33333
##
## Multivariate Tests: C(Geschlecht, contr.helmert):C(Medikament, contr.helmert)
## Df test stat approx F num Df den Df Pr(>F)
## Pillai 2 0.2269491 1.151992 4 36 0.34808
## Wilks 2 0.7743623 1.159326 4 34 0.34590
## Hotelling-Lawley 2 0.2896916 1.158766 4 32 0.34728
## Roy 2 0.2837227 2.553504 2 18 0.10562
##
## $Effektstärke
## eta^2
## (Intercept) 0.960659099
## C(Geschlecht, contr.helmert) 0.007463063
## C(Medikament, contr.helmert) 0.440189051
## C(Geschlecht, contr.helmert):C(Medikament, contr.helmert) 0.113474526
##
## $ANOVAs
## $ANOVAs$x1
## Anova Table (Type III tests)
##
## Response: dv
## Sum Sq Df F value Pr(>F)
## (Intercept) 2281.50 1 434.5714 4.704e-14 ***
## Geschlecht 0.67 1 0.1270 0.7257
## Medikament 301.00 2 28.6667 2.538e-06 ***
## Geschlecht:Medikament 14.33 2 1.3651 0.2806
## Residuals 94.50 18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## $ANOVAs$x2
## Anova Table (Type III tests)
##
## Response: dv
## Sum Sq Df F value Pr(>F)
## (Intercept) 1802.67 1 284.6316 1.776e-12 ***
## Geschlecht 0.67 1 0.1053 0.74934
## Medikament 36.33 2 2.8684 0.08292 .
## Geschlecht:Medikament 32.33 2 2.5526 0.10569
## Residuals 114.00 18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Grafiken
# HE-Plots
if (length(AV)==2){
heplot(fit)
} else{
pairs(fit)
}
# Diskriminazfunktionsplots
library(candisc)
term <- "Geschlecht:Medikament"
model <- as.formula(paste("cbind(",paste(AV,collapse=","),")~",paste(Faktoren,collapse="*"),sep=""))
fit2 <- lm(model, data=data2)
can <- candisc(fit2, term=term)
plot(can)
