Hintergrund: Wir betrachten den Datensatz airquality (Luftqualitätsmessungen in New York).
Mehr Infos sehen wir über:
Frage: Ist die Ozonbelastung (Ozone) im Hochsommer (August) signifikant höher als im Frühling (Mai)?
Wir folgen dem Workflow der Datenanalyse: [Link]
Bezug Glossar: Variablen & Datentypen [Link]
Zuerst laden wir die Daten und verschaffen uns einen Überblick über die Datentypen [Link].
# Datensatz laden (in R vorinstalliert)
data("airquality")
# Struktur prüfen: Welche Datentypen liegen vor?
str(airquality)'data.frame': 153 obs. of 6 variables:
$ Ozone : int 41 36 12 18 NA 28 23 19 8 NA ...
$ Solar.R: int 190 118 149 313 NA NA 299 99 19 194 ...
$ Wind : num 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
$ Temp : int 67 72 74 62 56 66 65 59 61 69 ...
$ Month : int 5 5 5 5 5 5 5 5 5 5 ...
$ Day : int 1 2 3 4 5 6 7 8 9 10 ...Hinweis: Die Variable Month wird als ganze Zahl (Integer) erkannt (5, 6, 7…). Für unseren Vergleich müssen wir später die Monate 5 und 8 herausfiltern.
Für uns ist die Spalte Ozone zusätzlich interessant. Diese ist eine ganzzahlige numerische Variable und wir ordnen den Ozonmesswert zu den Intervall-skalierten Variablentypen.
Bezug Glossar: Datenqualität sichern [Link]
Echte Daten sind selten sauber. Wir prüfen auf fehlende Werte (NA) und bereinigen [Link] den Datensatz falls nötig.
# Gibt es fehlende Werte in der Spalte Ozone?
table(is.na(airquality$Ozone))
FALSE TRUE
116 37 Wir sehen es gibt 37 Tage an denen kein Ozonwert eingetragen ist. Diese Fehlwerte müssen eliminiert werden.
Da wir einen Hypothesentest durchführen wollen, können wir die betroffenen Zeilen einfach löschen, da Tests nicht von zeitlichen Abhängigkeiten beeinflusst werden.
# Wir erstellen einen sauberen Teil-Datensatz nur für Mai (5) und August (8)
# und entfernen dabei Zeilen, die NA bei Ozone haben (!is.na)
analyse_daten <- subset(airquality,
Month %in% c(5, 8) & !is.na(Ozone))
# Zudem überprüfen wir, ob der Wertebereich der Ozone-Werte realistisch ist.
# Das geht am einfachsten mit der Funktion 'summary()'
summary(analyse_daten) Ozone Solar.R Wind Temp Month
Min. : 1.00 Min. : 8.0 Min. : 2.30 Min. :57.00 Min. :5.0
1st Qu.: 15.50 1st Qu.: 80.5 1st Qu.: 7.40 1st Qu.:66.00 1st Qu.:5.0
Median : 30.50 Median :203.0 Median : 9.70 Median :76.50 Median :6.5
Mean : 41.79 Mean :177.7 Mean :10.01 Mean :75.35 Mean :6.5
3rd Qu.: 60.50 3rd Qu.:255.5 3rd Qu.:12.00 3rd Qu.:82.00 3rd Qu.:8.0
Max. :168.00 Max. :334.0 Max. :20.10 Max. :97.00 Max. :8.0
NA's :5
Day
Min. : 1.00
1st Qu.: 7.75
Median :16.00
Mean :15.63
3rd Qu.:22.25
Max. :31.00
Der Wertebereich von Min = 1 bis Max = 168 sieht gut aus. Hier benötigen wir keine weitere Anpassung (keine negativen Werte o.ä.).
Mit dem summary()-Befehl erhalten wir zudem erste Informationen über die primären Lagemaße [Link] und Quantile [Link].
# Kurzer Check: Wie viele Messungen bleiben übrig?
table(analyse_daten$Month)
5 8
26 26 Wir sehen dass wir pro Monat noch 26 Werte zum Vergleich vorliegen haben.
Bezug Glossar: Grafische Darstellung
Bevor wir zum testen kommen, schauen wir uns die Verteilungen an. Ein Boxplot [Link] eignet sich perfekt für den Vergleich von Gruppen.
# Boxplot: Ozon-Werte (y) abhängig vom Monat (x)
boxplot(Ozone ~ Month,
data = analyse_daten,
main = "Vergleich Ozonwerte: Mai vs. August",
xlab = "Monat",
ylab = "Ozon (ppb)",
col = c("lightblue", "orange"))
Visuelle Interpretation: Die orange Box (August) liegt deutlich höher als die blaue (Mai). Auch die Streuung (Länge der Box) scheint im August größer zu sein.
Bezug Glossar: Verteilungs-Checks
Für den klassischen t-Test [Link] sollten die Daten in den Gruppen annähernd normalverteilt sein. Wir prüfen das zuerst grafisch mit einem Histogramm [Link].
# Bildschirm teilen für zwei Grafiken nebeneinander
par(mfrow = c(1, 2))
# Histogramm für Mai
hist(analyse_daten$Ozone[analyse_daten$Month == 5],
main = "Verteilung Mai", xlab = "Ozon", col = "lightblue")
# Histogramm für August
hist(analyse_daten$Ozone[analyse_daten$Month == 8],
main = "Verteilung August", xlab = "Ozon", col = "orange")
# Layout zurücksetzen
par(mfrow = c(1, 1))Die Daten sind offensichtlich nicht glockenförmig (eher rechtsschief). Dennoch führen wir einen Shapiro-Wilk-Test [Link] durch, da dieser uns eine statistisch kräftigere Aussage zur Normalverteilungsannahme der einzelnen Stichproben gibt, als rein visuelle Vergleiche.
# Test auf Normalverteilung
# H0: Die Stichprobe stammt aus einer Normalverteilung.
# Test für Mai
shapiro.test(analyse_daten$Ozone[analyse_daten$Month == 5])
Shapiro-Wilk normality test
data: analyse_daten$Ozone[analyse_daten$Month == 5]
W = 0.71401, p-value = 8.294e-06# Test für August
shapiro.test(analyse_daten$Ozone[analyse_daten$Month == 8])
Shapiro-Wilk normality test
data: analyse_daten$Ozone[analyse_daten$Month == 8]
W = 0.93279, p-value = 0.09032Ergebnisauswertung: Beide p-Werte sind deutlich kleiner als 0.05 (Mai: p < .001, August: p = .090). Das bedeutet, wir müssen die Nullhypothese verwerfen. Die Daten sind nicht normalverteilt. Streng genommen ist der t-Test hier also nicht ideal, weshalb wir später den Mann-Whitney-U-Test betrachten. Da der t-Test jedoch robust gegen diese Abweichung ist, führen wir ihn zu Demonstrationszwecken dennoch durch.
Neben der Normalverteilungsannahme wird für den t-Test oft Varianzhomogenität vorausgesetzt. Diese prüfen wir mittels des Levene-Tests [Link].
# Test auf Varianzgleichheit
#install.packages("car") # vorab Paket installieren (falls nötig)
library("car") # und anschließend ladenLade nötiges Paket: carData# Die Stichproben weisen eine homogene Varianz auf.
leveneTest(Ozone ~ as.factor(Month), data = analyse_daten)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 9.5019 0.003336 **
50
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Der Levene-Test liefert einen p-Wert von 0.003, was deutlich unter dem Signifikanzniveau von 0.05 liegt. Das bedeutet: Wir müssen die Nullhypothese (Varianzgleichheit) verwerfen. Die Varianzen der beiden Gruppen sind signifikant unterschiedlich (heterogen).
Konsequenz für den t-Test: Wir dürfen nicht den klassischen Student-t-Test rechnen (der gleiche Varianzen annimmt), sondern müssen den Welch-t-Test verwenden. In R erreichen wir das, indem wir im Befehl t.test() das Argument var.equal = FALSE setzen (was ohnehin der sicherere Standard in R ist).
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Wir berechnen einen t-Test für unabhängige Stichproben [Link], da die Messungen im Mai und August nicht voneinander abhängen [Link: Unabhängig vs. Verbunden].
# Durchführung des t-Tests für unverbundene Stichproben
test_ergebnis <- t.test(Ozone ~ Month, data = analyse_daten,
var.equal = F)
# Ergebnis ausgeben
test_ergebnis
Welch Two Sample t-test
data: Ozone by Month
t = -4.0749, df = 39.279, p-value = 0.0002169
alternative hypothesis: true difference in means between group 5 and group 8 is not equal to 0
95 percent confidence interval:
-54.38358 -18.30873
sample estimates:
mean in group 5 mean in group 8
23.61538 59.96154 Nun schauen wir abschließend auf den Output. Dabei ist in erster Linie der p-Wert [Link] relevant zur Beantwortung unserer Anfangshypothese. Dieser Wert liegt weit unter 0.05.
Zusätzlich werden am Ende des Outputs noch die Mittelwerte der beiden Teilgruppen Mai und August angegeben. Wir können hier erkennen, dass die Ozonwerte weit voneinander entfernt liegen.
Ergebnisformulierung: Die Analyse zeigt einen signifikanten Unterschied in der Ozonbelastung zwischen den Monaten Mai und August (t = -4.07, p < .001).
Um die statistische Aussage zu präzisieren und konkret unsere Forschungsfrage zu beantworten, führen wir den Test einseitig (gerichtet) durch. Wir prüfen, ob August-Werte signifikant höher sind. Da im Test “Mai” (5) die erste Gruppe ist, prüfen wir, ob “Mai < August” ist (alternative = "less")
# Durchführung des einseitigen t-Tests für unverbundene Stichproben
test_ergebnis_oneway <- t.test(Ozone ~ Month, data = analyse_daten,
alternative = "less", var.equal = T)
# Ergebnis ausgeben
test_ergebnis_oneway
Two Sample t-test
data: Ozone by Month
t = -4.0749, df = 50, p-value = 8.226e-05
alternative hypothesis: true difference in means between group 5 and group 8 is less than 0
95 percent confidence interval:
-Inf -21.39781
sample estimates:
mean in group 5 mean in group 8
23.61538 59.96154 Ergebnisformulierung: Die Analyse bestätigt, dass die Ozonbelastung im August signifikant höher ist als im Mai (t = -4.07, p < .001).
Bezug Glossar: Parametrische Verfahren (Mittelwerte)
Wie wir unter Kapitel 3 sehen konnten, waren die Voraussetzungen für den t-Test verletzt. Deshalb greifen wir zur sichereren, nicht-parametrischen Alternative. Im Falle von zwei unverbundenen Stichproben ist diese der Mann-Whitney-U-Test [Link] (in R als wilcox.test implementiert).
# Durchführung des Mann-Whitney-U-Tests
test_ergebnis_npar <- wilcox.test(Ozone ~ Month, data = analyse_daten,
exact = T) Warning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): kann bei
Bindungen keinen exakten p-Wert Berechnen# Exact wird auf TRUE gesetzt, wenn weniger als 40 Beobachtungen
# vorliegen, wie es hier der Fall ist
# Ergebnis ausgeben
test_ergebnis_npar
Wilcoxon rank sum test with continuity correction
data: Ozone by Month
W = 127.5, p-value = 0.0001208
alternative hypothesis: true location shift is not equal to 0Auswertung: Auch der nicht-parametrische Test bestätigt das Ergebnis. Der p-Wert liegt wieder weit unter dem Signifikanzniveau von 0.05.
Ergebnisformulierung: Der Mann-Whitney-U-Test bestätigt, dass sich die Ozonbelastung zwischen den Monaten Mai und August signifikant unterscheidet (W = 127.5, p < .001).
Frage einfach: Sind die Windgeschwindigkeiten (Wind) im Mai und im September signifikant voneinander verschieden?
Frage fortgeschritten: Ist an Tagen mit erhöhten Temperaturen (über 20°C) (Temp) die Globalstrahlung (Solar.R) signifikant höher?
Wähle dir nun eine der Forschungsfragen aus und wende das gerade dargestellte Prozedere auf diese Frage an.
# Prüfe die Struktur der Variablen
str(airquality)'data.frame': 153 obs. of 6 variables:
$ Ozone : int 41 36 12 18 NA 28 23 19 8 NA ...
$ Solar.R: int 190 118 149 313 NA NA 299 99 19 194 ...
$ Wind : num 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
$ Temp : int 67 72 74 62 56 66 65 59 61 69 ...
$ Month : int 5 5 5 5 5 5 5 5 5 5 ...
$ Day : int 1 2 3 4 5 6 7 8 9 10 ...# Month kennen wir schon vom Beispiel der Lerneinheit
# Wind ist eine numerische intervall-skalierte Variable
# Gibt es Fehlwerte? (Month kennen wir schon)
table(is.na(airquality$Wind))
FALSE
153 # keine Fehlwerte bei den Wind-Messwerten enthalten
# Muss der Datensatz für die Frage angepasst werden?
# --> Ja, wir filtern den Datensatz nach den relevanten Monaten (Mai und September):
analyse_daten1 <- subset(airquality,
Month %in% c(5, 9))
# Schau dir den Wertebereich der relevanten Variablen an. Gibt es Auffälligkeiten?
summary(analyse_daten1) Ozone Solar.R Wind Temp Month
Min. : 1.00 Min. : 8 Min. : 2.80 Min. :56.00 Min. :5.000
1st Qu.: 13.50 1st Qu.: 92 1st Qu.: 8.00 1st Qu.:64.00 1st Qu.:5.000
Median : 21.00 Median :193 Median :10.90 Median :70.00 Median :5.000
Mean : 27.75 Mean :174 Mean :10.91 Mean :71.13 Mean :6.967
3rd Qu.: 33.00 3rd Qu.:252 3rd Qu.:13.20 3rd Qu.:77.00 3rd Qu.:9.000
Max. :115.00 Max. :334 Max. :20.10 Max. :93.00 Max. :9.000
NA's :6 NA's :4
Day
Min. : 1.00
1st Qu.: 8.00
Median :16.00
Mean :15.75
3rd Qu.:23.00
Max. :31.00
table(analyse_daten1$Month)
5 9
31 30 # Der Wertebereich der Windgeschwindigkeiten sieht realistisch aus
# Und wir haben diesmal 31 bzw 30 Daten pro Monat (vollständig)# Prüfe die Struktur der Variablen
str(airquality)'data.frame': 153 obs. of 6 variables:
$ Ozone : int 41 36 12 18 NA 28 23 19 8 NA ...
$ Solar.R: int 190 118 149 313 NA NA 299 99 19 194 ...
$ Wind : num 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
$ Temp : int 67 72 74 62 56 66 65 59 61 69 ...
$ Month : int 5 5 5 5 5 5 5 5 5 5 ...
$ Day : int 1 2 3 4 5 6 7 8 9 10 ...Wir erkennen über str(airquality) (oder ?airquality oder spätestens bei summary(airquality)), dass die Temperatur nicht in °C sondern in °F angegeben ist. Diese muss erstmal umgerechnet [Link] werden um die Fragestellung zu beantworten. Die Formel dafür lautet: \((x °F − 32) × 5/9 = y °C\)
Außerdem erfordert die Frage eine Kategorisierung (Binarisierung) [Link: zu dummy vars?] der Temperaturvariable.
# 1. Einheit anpassen (Neue Spalte Temp_C)
airquality$Temp_C <- (airquality$Temp - 32) * 5/9
# 2. Temperatur binarisieren (0=kühler, 1=wärmer)
# Wir erstellen eine Gruppierungsvariable 'Temp_Gruppe'
airquality$Temp_Gruppe <- ifelse(airquality$Temp_C > 28, 1, 0)Wie bei den Analysebeispielen davor muss natürlich auch hier wieder auf Vollständigkeit überprüft und falls nötig der Datensatz entsprechend bereinigt werden.
# Gibt es Fehlwerte?
table(is.na(airquality$Temp)) # oder Temp_C oder Temp_Gruppe, kommt aufs selbe hinaus
FALSE
153 table(is.na(airquality$Solar.R))
FALSE TRUE
146 7 # 7 Fehlwerte bei der Globalstrahlung
# Muss der Datensatz für die Frage angepasst werden?
# --> Ja, wir entfernen die Fehlwerte
analyse_daten2 <- subset(airquality, !is.na(Solar.R))
# Schau dir den Wertebereich der relevanten Variablen an. Gibt es Auffälligkeiten?
summary(analyse_daten2) Ozone Solar.R Wind Temp Month
Min. : 1.0 Min. : 7.0 Min. : 1.7 Min. :57.00 Min. :5.000
1st Qu.: 18.0 1st Qu.:115.8 1st Qu.: 7.4 1st Qu.:73.00 1st Qu.:6.000
Median : 31.0 Median :205.0 Median : 9.7 Median :79.00 Median :7.000
Mean : 42.1 Mean :185.9 Mean :10.0 Mean :78.12 Mean :7.027
3rd Qu.: 62.0 3rd Qu.:258.8 3rd Qu.:11.5 3rd Qu.:84.00 3rd Qu.:8.000
Max. :168.0 Max. :334.0 Max. :20.7 Max. :97.00 Max. :9.000
NA's :35
Day Temp_C Temp_Gruppe
Min. : 1.00 Min. :13.89 Min. :0.0000
1st Qu.: 9.00 1st Qu.:22.78 1st Qu.:0.0000
Median :16.00 Median :26.11 Median :0.0000
Mean :16.12 Mean :25.62 Mean :0.3082
3rd Qu.:23.75 3rd Qu.:28.89 3rd Qu.:1.0000
Max. :31.00 Max. :36.11 Max. :1.0000
table(analyse_daten2$Temp_Gruppe)
0 1
101 45 # Der Wertebereich der Solardaten sieht realistisch aus (Temperatur hatten wir ja zuvor angepasst und binarisiert)Wir haben eine stark ungleiche Anzahl an Beobachtungen pro Stichprobe (101 kühler vs 45 wärmer). Das ist für den t-Test für unverbundene Stichproben prinzipiell kein Problem, macht die Prüfung der Varianzhomogenität aber umso wichtiger.
# Erstelle die Boxplots
boxplot(Wind ~ Month,
data = analyse_daten1,
main = "Vergleich Windgeschwindigkeit: Mai vs. September",
xlab = "Monat",
ylab = "Wind",
col = c("lightblue", "orange"))
# Im September sind die Windgeschwindigkeiten durchschnittlich niedriger.
# Aber ist dieser Unterschied auch signifikant?# Erstelle die Boxplots
boxplot(Solar.R ~ Temp_Gruppe,
data = analyse_daten2,
main = "Vergleich Globalstrahlung: Temperaturen unter vs. über 28°C",
xlab = "Temperaturklasse",
ylab = "Globalstrahlung",
col = c("lightblue", "orange"))
# An wärmeren Tagen (1) scheint die Globalstrahlung im Median höher zu sein und weniger zu streuen.
# Aber ist dieser Unterschied auch signifikant?# grafisch (Histogramme)
hist(analyse_daten1$Wind[analyse_daten1$Month == 5],
main = "Verteilung Mai", xlab = "Wind", col = "lightblue")
hist(analyse_daten1$Wind[analyse_daten1$Month == 9],
main = "Verteilung September", xlab = "Wind", col = "orange")
# Visuelle Überprüfung zeigt eine relativ symmetrische Verteilung mit Tendenz zur rechtsschiefe.
# Da der Stichprobenumfang gering ist, ist eine perfekte Glockenkurve selten.# statistisch 1 (Normalverteilung?)
shapiro.test(analyse_daten1$Wind[analyse_daten1$Month == 5])
Shapiro-Wilk normality test
data: analyse_daten1$Wind[analyse_daten1$Month == 5]
W = 0.968, p-value = 0.4659shapiro.test(analyse_daten1$Wind[analyse_daten1$Month == 9])
Shapiro-Wilk normality test
data: analyse_daten1$Wind[analyse_daten1$Month == 9]
W = 0.97853, p-value = 0.7852# Normalverteilungsannahme wird für beide Teilstichproben nicht verworfen.
# statistisch 2 (Varianzhomogenität?)
leveneTest(Wind ~ as.factor(Month), data = analyse_daten1)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 0.0601 0.8072
59 # Annahme der Varianzhomogenität wird ebenfalls nicht verworfen.
# In diesem Beispiel sind die Voraussetzungen für den t-Test nicht verletzt
# und wir können den Mittelwertsvergleich ohne Bedenken darüber durchführen.# grafisch (Histogramme)
hist(analyse_daten2$Solar.R[analyse_daten2$Temp_Gruppe == 0],
main = "Verteilung lower 28°C", xlab = "Globalstrahlung", col = "lightblue")
hist(analyse_daten2$Solar.R[analyse_daten2$Temp_Gruppe == 1],
main = "Verteilung higher 28°C", xlab = "Globalstrahlung", col = "orange")
# Visuelle Überprüfung zeigt bei den niedrigeren Temperaturen eine gleichmäßigere Verteilung.# statistisch 1 (Normalverteilung?)
shapiro.test(analyse_daten2$Solar.R[analyse_daten2$Temp_Gruppe == 0])
Shapiro-Wilk normality test
data: analyse_daten2$Solar.R[analyse_daten2$Temp_Gruppe == 0]
W = 0.93213, p-value = 6.062e-05shapiro.test(analyse_daten2$Solar.R[analyse_daten2$Temp_Gruppe == 1])
Shapiro-Wilk normality test
data: analyse_daten2$Solar.R[analyse_daten2$Temp_Gruppe == 1]
W = 0.96562, p-value = 0.1995# Normalverteilungsannahme wird für eine der beiden Teilstichproben verworfen.
# statistisch 2 (Varianzhomogenität?)
leveneTest(Solar.R ~ as.factor(Temp_Gruppe), data = analyse_daten2)Levene's Test for Homogeneity of Variance (center = median)
Df F value Pr(>F)
group 1 22.463 5.09e-06 ***
144
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Zusätzlich wird die Annahme der Varianzhomogenität ebenfalls verworfen.In diesem Beispiel ist die Wahl des nicht-parametrischen Tests die methodisch-sauberste Wahl.
Exkurs: Möchte man dennoch einen t-Test rechnen (da er robust ist), muss zwingend der Welch-t-Test [Link] (var.equal = FALSE) genutzt werden, um die Varianzheterogenität zu korrigieren.
# Test-Durchführung: (parametrisch)
t.test(Wind ~ Month, data = analyse_daten1, var.equal = T)
Two Sample t-test
data: Wind by Month
t = 1.6107, df = 59, p-value = 0.1126
alternative hypothesis: true difference in means between group 5 and group 9 is not equal to 0
95 percent confidence interval:
-0.3495943 3.2347556
sample estimates:
mean in group 5 mean in group 9
11.62258 10.18000 Ergebnisformulierung: Die Analyse zeigt keinen signifikanten Unterschied in den Windgeschwindigkeiten zwischen den Monaten Mai und September (t = 1.61, p < .113).
Da wir eine gerichtete Hypothese haben (“Ist bei warmem Wetter die Strahlung höher?”), testen wir einseitig. R ordnet die Gruppen als 0 (kühler) und 1 (wärmer). Hypothese: Gruppe 0 < Gruppe 1. Daher nutzen wir alternative = "less".
# Test-Durchführung: (einseitig nichtparametrisch)
wilcox.test(Solar.R ~ Temp_Gruppe, data = analyse_daten2, alternative = "less")
Wilcoxon rank sum test with continuity correction
data: Solar.R by Temp_Gruppe
W = 1746.5, p-value = 0.01297
alternative hypothesis: true location shift is less than 0Ergebnisformulierung: Die Analyse unterstützt die These, dass an Tagen mit erhöhten Temperaturen (über 28°C) die Globalstrahlung signifikant höher ist als an kühlen Tagen (W = 1746.5, p < .013).