Einfaktorielle ANOVA in R
Quantis
24 5 2022
1 Datensätze in R
In R gibt es einige Standarddatensätze (die bei dem Programm “mit dabei” sind), wie z.B. den Datensatz “ToothGrwoth”: Ein Datensatz der zu einem Experiment gehört, in dem der Einfluss von Vitamin C auf das Zahnwachstum bei Meerschweinchen untersucht wurde).
Eine Übersicht über die Standarddatensätze gibt es unter folgendem Link http://www.sthda.com/english/wiki/r-built-in-data-sets.
2 Zahnwachstumsrate bei Meerschweinchen
Wenn man genauer wissen will, was genau das Szenario hinter den Daten ist, kann man einfach “?NAMEDATENSATZ” als Code eingeben, hier also z.B. “?ToothGrowth”
#?ToothGrowthIm Hilfefenster würde man dann folgende Information sehen:
The response is the length of odontoblasts (cells responsible for tooth growth) in 60 guinea pigs. Each animal received one of three dose levels of vitamin C (0.5, 1, and 2 mg/day) by one of two delivery methods, orange juice or ascorbic acid (a form of vitamin C and coded as VC).
Heißt: Es gab drei Dosis-Bedingungen. Außerdem wurde noch variiert, in welcher Weise das Vitamin C verabreicht wurde (Orangensaft vs. Askorbinsäure [wahrscheinlich aufgelöst in Wasser]). Gemessen (Kriterium bzw. AV) wurde die Länge der Odontoblasten (die Zellen, die für Zahnlänge verantwortlich sind).
3 Einfaktorielle ANOVA: Der Einfluss der Vitmain-C-Dosis
3.1 Aufgaben
Der Zweite Faktor (die Art der Verabreichung) soll nachfolgend zunächst ignoriert werden.
Testen Sie die folgenden Hypothesen:
Hypothese 1: Die Odontoblastenlänge hängt von der Dosierung ab (einfaktorielle ANOVA).
Hypothese 2: Die Odontoblastenlänge nimmt mit der Dosierung zu (Post-hoc Tests).
Beginnen Sie die statistische Auswertung mit einer grafischen Veranschaulichung der Daten. Erstellen Sie:
Eine Abbildung, die die Gruppenmittelwerte und Fehlerbalken zeigt,
Verteilungsabbildungen: (1) Boxplot, (2) Strip-Chart und (3) Density-Plot.
Für den Strip-Chart: Diese Art Abbildung hatten wir bisher nicht erstellt. Benutzen Sie deshalb zur Erstllung die Hilfe in R. Tippen Sie dazu ?stripchart in die Console, um das Übersichtsmenü aufzurufen. Nutzen Sie zur Erstellung der Grafik den Funktionsparameter “jitter” (method = “jitter”). Finden Sie einen Jitter-Wert, der eine angemessene Darstellung der Datenpunkte ermöglicht. Färben Sie die Datenpunkte blau. Bennen Sie auch die Achsen der Grafik.
Welchen Eindruck vermittelt Ihnen der Strip-Chart bzgl. des Kriteriums “Varianzhomogenität”?
Prüfen Sie dieses Kriterium optisch auch anhand des Density-Plots.
Testen Sie das Kriterium danach mit einem Test (s. Formelsammlung).
Prüfen Sie dann die aufgestellten Hypothesen (ANOVA und post-hoc-Tets)
Zusatz:
Erstellen Sie abschließend ein lineares Modell mit Kriterium Odontoblastenlänge und Prädiktr Dosierung. Vergleichen Sie den Output mit dem der Varianzanalyse und dem der Post-hoc Tests für Mittelwertsunterschiede.
3.2 Datensatz laden
data("ToothGrowth") # "Einlesen" der Daten (hier besser gesagt "aktivieren" der Daten)
head(ToothGrowth) # Anzeige der ersten 10 Einträge## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.53.3 Faktor erzeugen
die experimentelle Variable “dose” (also Dosierung von Vitamin C) ist hier numerisch. Aber eigentlich ist es ja ein experimenteller Faktor. Also zunächst diese Variable in einen Faktor konvertieren.
ToothGrowth <- within(ToothGrowth, {
dose <- factor(dose, labels=c('low','medium','large'))
})3.4 Grafiken
3.4.1 Plot für Mittelwerte
library(Rcmdr)with(ToothGrowth, plotMeans(len, dose, error.bars="conf.int", level=0.95, xlab="Dosierung",
ylab="Odontoblastenlänge", main="Mittlere Odontoblastenlänge (Fehler = 95% KI)", connect=FALSE))
3.4.2 Boxplot
Boxplot(len~dose, data=ToothGrowth, id=list(method="y"))
3.4.3 Strip-Chart
stripchart(len ~ dose, vertical=TRUE, method="jitter", jitter = 0.1, ylab="Odontoblastenlänge",
xlab = "Dosierung", col = "blue", data=ToothGrowth)
3.4.4 Density Plot
densityPlot(len~dose, data=ToothGrowth, bw=bw.SJ, adjust=1, kernel=dnorm, method="adaptive",
xlab="Dosierung", ylab="Odontoblastenlänge", main = "Verteilung der Odontoblastenlängen")
3.5 Test auf Varianzhomogenität
Der nötige Test heißt “Levene-Test”. ACHTUNG: der R-Commander lädt dazu einige Pakete.
library(mvtnorm, pos=17)
library(survival, pos=17)
library(MASS, pos=17)
library(TH.data, pos=17)##
## Attaching package: 'TH.data'## The following object is masked from 'package:MASS':
##
## geyserlibrary(multcomp, pos=17)
library(abind, pos=22)
Tapply(len ~ dose, var, na.action=na.omit, data=ToothGrowth) # variances by group## low medium large
## 20.24787 19.49608 14.24421leveneTest(len ~ dose, data=ToothGrowth, center="mean") # Test auf Var-Homog.## Levene's Test for Homogeneity of Variance (center = "mean")
## Df F value Pr(>F)
## group 2 0.7328 0.485
## 57Der Levene-Test ist nicht signifikant. Das Kriterium der Varianzhomogenität kann als erfüllt betrachtet werden.
3.6 Statistische Hypothesentests
3.6.1 Einfaktorielle ANOVA
Zur Prüfung der Hypothese, dass die Dosierung einen Einfluss hat.
Im R-Commander-Menü dazu kann bereits der paarweise Vergleich der Mittelwerte aktiviert werden (Post-hoc-Tests). Hier aber zunächst nur der Output der Code und Output der ANOVA.
AnovaModel.1 <- aov(len ~ dose, data=ToothGrowth) # Erstellt das ANOVA-Modell
summary(AnovaModel.1) # Gibt den Output aus## Df Sum Sq Mean Sq F value Pr(>F)
## dose 2 2426 1213 67.42 9.53e-16 ***
## Residuals 57 1026 18
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Die einfaktorielle ANOVA ist signifikant, \(F(2,57)~= 67.42\), \(p~< .001\). Die Hypothese, dass die Dosierung mit einer veränderten Odontoblastenlänge einhergeht, kann als bestätigt betrachtet werden.
Die vom R-Commander genutzte Funktion der Einfaktoriellen ANOVA liefert keine Effektgrößen. Ein alternatives Paket, das das tut, ist das Paket “afex()”. Das müssen Sie ggf. installieren. Um damit eine ANOVA durchzuführen muss allerdings eine ID-Spalte in den Daten sein, die die einzelnen Beobachtungen kennzeichnet. In unserem Fall müssen wir diese zunächst den Daten hinzufügen.
id <- c(1:length(ToothGrowth$len))
ToothGrowth$id <- idlibrary(afex)## Loading required package: lme4## Loading required package: Matrix## ************
## Welcome to afex. For support visit: http://afex.singmann.science/## - Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
## - Methods for calculating p-values with mixed(): 'KR', 'S', 'LRT', and 'PB'
## - 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
## - NEWS: library('emmeans') now needs to be called explicitly!
## - Get and set global package options with: afex_options()
## - Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
## - For example analyses see: browseVignettes("afex")
## ************##
## Attaching package: 'afex'## The following object is masked from 'package:lme4':
##
## lmera1 <- aov_car(len ~ dose + Error(id), ToothGrowth, anova_table = list(es = "pes"))## Contrasts set to contr.sum for the following variables: dosea1## Anova Table (Type 3 tests)
##
## Response: len
## Effect df MSE F pes p.value
## 1 dose 2, 57 18.00 67.42 *** .703 <.001
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1# Die Effektstärke "pes" ist das partielle Eta-Quadrate (partial eta squared), welches in G-Power verwendet werden kann.3.6.2 Post-hoc Tests
Der folgende Code stammt aus dem R-Commander und kann erzeugt werden, wenn man bei der einfaktoriellen ANOVA die Option “paarweiser Vergleich” der Mittelwerte wählt.
Der Code erzeugt zwei Ergebnisse:
Die Mittelwerte sowie °>Tests für die Vergleiche aller möglichen Mittelwertsunterschiede (Hier wird automatisch für multiple Vgl. kontrolliert; Methode = TUKEY).
Eine grafische Zusammenfassung der Mittelwertsdifferenzen (für alle möglichen Vergleiche).
with(ToothGrowth, numSummary(len, groups=dose, statistics=c("mean", "sd")))## mean sd data:n
## low 10.605 4.499763 20
## medium 19.735 4.415436 20
## large 26.100 3.774150 20local({
.Pairs <- glht(AnovaModel.1, linfct = mcp(dose = "Tukey"))
print(summary(.Pairs)) # pairwise tests
print(confint(.Pairs, level=0.95)) # confidence intervals
print(cld(.Pairs, level=0.05)) # compact letter display
old.oma <- par(oma=c(0, 5, 0, 0))
plot(confint(.Pairs))
par(old.oma)
})##
## Simultaneous Tests for General Linear Hypotheses
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = len ~ dose, data = ToothGrowth)
##
## Linear Hypotheses:
## Estimate Std. Error t value Pr(>|t|)
## medium - low == 0 9.130 1.341 6.806 < 1e-05 ***
## large - low == 0 15.495 1.341 11.551 < 1e-05 ***
## large - medium == 0 6.365 1.341 4.745 3.94e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## (Adjusted p values reported -- single-step method)
##
##
## Simultaneous Confidence Intervals
##
## Multiple Comparisons of Means: Tukey Contrasts
##
##
## Fit: aov(formula = len ~ dose, data = ToothGrowth)
##
## Quantile = 2.4061
## 95% family-wise confidence level
##
##
## Linear Hypotheses:
## Estimate lwr upr
## medium - low == 0 9.1300 5.9022 12.3578
## large - low == 0 15.4950 12.2672 18.7228
## large - medium == 0 6.3650 3.1372 9.5928
##
## low medium large
## "a" "b" "c"
Die Post-hoc-Tests bestätigen die zweite Hypothese. Höhere Dosierung geht mit längeren Odontoblasten einher.
3.7 Zusatz
Anpassung eines linearen Modells (Regressionsanalyse)
LinearModel.2 <- lm(len ~ dose, data=ToothGrowth)
summary(LinearModel.2)##
## Call:
## lm(formula = len ~ dose, data = ToothGrowth)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.6000 -3.2350 -0.6025 3.3250 10.8950
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.6050 0.9486 11.180 5.39e-16 ***
## dosemedium 9.1300 1.3415 6.806 6.70e-09 ***
## doselarge 15.4950 1.3415 11.551 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.242 on 57 degrees of freedom
## Multiple R-squared: 0.7029, Adjusted R-squared: 0.6924
## F-statistic: 67.42 on 2 and 57 DF, p-value: 9.533e-16Der F-Test für das Gesamtmodell entspricht dem F-Test der ANOVA. Das multiple R-Quadrat entspricht dem partiellen Eta-Quadrat. Die Estimates korrespondieren mit den Mittelwertsvergleichen der post-hoc-Tests:
- Intercept = Mittelwert der Gruppe mit niedriger Dosierung
- dose[T.medium] = Differenz medium und low dose
- dose[T.large] = Differenz large und low dose
Es fehlt hier also der Vgl. von large und medium dose.
(ACHTUNG: Die beschriebenen Übereinstimmungen gelten für die einfaktorielle ANOVA).
3.7.1 Mit dummy-kodierten Daten
Diese Schritte sind nur zur Veranschaulichung. R ist bereits selbst in der Lage, für eine Regressionsanalyse mit einem mehrfachgestuften Faktor eine Dummykodierung vorzunehmen.
# Erstellen von Dummy-Spalten:
ToothGrowth$med[ToothGrowth$dose == "low"] <- 0
ToothGrowth$med[ToothGrowth$dose == "medium"] <- 1
ToothGrowth$med[is.na(ToothGrowth$med)] <- 0
ToothGrowth$high[ToothGrowth$dose == "low"] <- 0
ToothGrowth$high[ToothGrowth$dose == "medium"] <- 0
ToothGrowth$high[is.na(ToothGrowth$high)] <- 1RegModel.3 <- lm(len~med+high, data=ToothGrowth)
summary(RegModel.3)##
## Call:
## lm(formula = len ~ med + high, data = ToothGrowth)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.6000 -3.2350 -0.6025 3.3250 10.8950
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.6050 0.9486 11.180 5.39e-16 ***
## med 9.1300 1.3415 6.806 6.70e-09 ***
## high 15.4950 1.3415 11.551 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.242 on 57 degrees of freedom
## Multiple R-squared: 0.7029, Adjusted R-squared: 0.6924
## F-statistic: 67.42 on 2 and 57 DF, p-value: 9.533e-16Analyse liefert denselben Output.