Übungszettel 8

1 Aufgaben

In dieser Übung wollen wir uns mit kontinuierlichen Daten beschäftigen. Hierzu betrachten wir den Datensatz airquality, welcher verschiedene Messungen bezüglich Luftqualität zwischen Mai und September im Jahre 1973 umfasst.
Für mehr Information siehe ?airquality.

library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
── Attaching packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──✔ ggplot2 3.4.0      ✔ purrr   0.3.5 
✔ tibble  3.1.8      ✔ dplyr   1.0.10
✔ tidyr   1.2.1      ✔ stringr 1.4.1 
✔ readr   2.1.3      ✔ forcats 0.5.2 ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()

airquality <- as_tibble(airquality)
# Faktoren aus Monat und Tag machen:
airquality <- airquality %>%
  mutate_at(c("Month", "Day"), factor)
# Alternativ:
# airquality %>%
#   mutate(Month = factor(Month), Day = factor(Day))
airquality

1.1 Übungsaufgaben

1. Schaffen Sie sich zunächst eine (statistische) Übersicht über die Daten. Wie viele Beobachtungen und Variablen wurden erfasst?
2. Versichern Sie sich zunächst, dass die Daten einigermaßen kontinuierlich innerhalb des Zeitraums erfasst wurden.
3. Veranschaulichen Sie die Verteilungen der kontinuierlichen Variablen von airquality. Welche Aspekte fallen ins Auge und könnten interessant sein?
4. Finden Sie Zusammenhänge zwischen den verschiedenen Variablen.
   1. Visualisieren Sie für alle sinnvollen Variablen die paarweisen Beziehungen. Sind diese alle linear?
   2. Berechnen Sie für alle sinnvollen Variablen eine Korrelationsmatrix. Bietet sich Pearsons Korrelationskoeffizient an?
   3. Bestimmen Sie durch Visualisierung, wo Zusammenhänge vorliegen.
5. Erstellen Sie für relevante (kausale) Zusammenhänge Regressionsmodelle.

1.2 Lösungen

1.2.1 Schaffen Sie sich zunächst eine (statistische) Übersicht über die Daten. Wie viele Beobachtungen und Variablen wurden erfasst?

Anhang des tibbles sehen wir direkt die Anzahl der Beobachtungen und Variablen. Alternativ auch oben rechts unter Environment oder mit dim()

airquality

dim(airquality)

[1] 153   6

summary(airquality)

     Ozone           Solar.R           Wind             Temp       Month       Day     
 Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00   5:31   1      :  5  
 1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00   6:30   2      :  5  
 Median : 31.50   Median :205.0   Median : 9.700   Median :79.00   7:31   3      :  5  
 Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88   8:31   4      :  5  
 3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00   9:30   5      :  5  
 Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00          6      :  5  
 NA's   :37       NA's   :7                                               (Other):123  

1.2.2 Versichern Sie sich zunächst, dass die Daten einigermaßen kontinuierlich innerhalb des Zeitraums erfasst wurden.

Hier gibt es mehrere Möglichkeiten. Naiv könnten wir die Tabelle durchschauen und jeden einzelnen Tag und Monat nachschlagen, ob dieser vorhanden ist. Hier kann aber leicht das Gesamtbild aus den Augen verloren gehen, daher wollen wir das kurz veranschaulichen:

ggplot(airquality) +
  geom_bar(aes(x = Day))

ggplot(airquality) +
  geom_bar(aes(x = Month))

Wenn wir Balkendiagramme für Tag und Monat separat zeichnen, deutet sich an, dass diese scheinbar kontinuierlich erfasst wurden.
Seien Sie hierbei vorsichtig, denn es wäre prinzipiell auch möglich, mehrere Messungen an demselben Tag durchzuführen, um das Ergebnis zu verändern.

Für eine deutlichere Einsicht können wir uns das gemeinsame Balkendiagramm anzeigen:

ggplot(airquality) +
  geom_bar(aes(x = Day, fill = Month))

ggplot(airquality) +
  geom_bar(aes(x = Month, fill = Day)) +
  facet_wrap(~ Day)

Scheinbar liegen die kleinen Unterschiede lediglich an dem 31. Tag im Monat. Die Daten wurden somit kontinuierlich erfasst.

1.2.3 Veranschaulichen Sie die Verteilungen der kontinuierlichen Variablen von airquality. Welche Aspekte fallen ins Auge und könnten interessant sein?

Wir könnten uns hier jede einzelne Variable im Detail anschauen, aber verbleiben beispielhaft bei der Variablen Ozone.

Schauen wir uns zunächst die allgemeine Verteilung anhand von Boxplot und Histogramm an:

ggplot(airquality) +
  geom_boxplot(aes(y = Ozone), fill = "steelblue")

ggplot(airquality, aes(x = Ozone)) +
  geom_histogram(fill = "steelblue", bins = 15) +
  geom_freqpoly(bins = 15)

Viele Werte scheinen sich zwischen ~10 und ~30 aufzuhalten. Das erste Quartil wird durch 18 und das dritte Quartial durch 63 markiert.

Aus der Auskunft (?airquality) sehen wir, dass es sich hierbei um durchschnittliche Teilchen pro Milliarde zwischen 13 und 15 Uhr, gemessen auf Roosevelt Island, handelt. Scheinbar liegen auch einige fehlende Werte und Ausreißer vor, welche wir zunächst betrachten wollen.

Ausreißer können wir erhalten, in dem wir (in diesem Fall) nach hohen Werte von Ozone filtern:

airquality %>%
  filter(Ozone > 125)

Wir sehen, dass am 1.07.1973 und am 25.08.1973 Ausreißer bei den Ozonwerten vorlagen. Um die Ausreißer der Boxplots direkt zu erhalten (ohne selbstständiges Filtern), können wir aus rstatix die Funktion identify_outliers() verwenden:

library(rstatix)

Attache Paket: ‘rstatix’

Das folgende Objekt ist maskiert ‘package:stats’:

    filter

airquality %>%
  identify_outliers(Ozone)

Wir sehen, dass dieselben zwei Werte in der Spalte is.outlier als Ausreißer (TRUE) markiert wurden, jedoch nicht als extreme Ausreißer (FALSE). Ausreißer weichen um einen IQR von mehr als 1.5 von dem 3. Quartil ab, extreme Ausreißer ab dem dreifachen IQR.

Um die fehlenden Einträge zu betrachten, können wir nach is.na() filtern:

airquality %>%
  filter(is.na(Ozone))

Scheinbar wurden an 2/3 der Tage des Junis keine Daten bezüglich der Ozonwerte erhoben. Wohlmöglich gab es Störungen der Messgeräte?

Nun stellt sich die Frage, wieso so viele Messwerte fehlen. Das lässt sich ohne weitere Information des Datenerhebers nicht einfach klären, allerdings können wir nach Mustern suchen:

airquality %>%
  group_by(Month) %>%
  summarise(n = sum(is.na(Ozone)))

Schauen wir, ob es an zusammenhängenden Zeitpunkten Probleme bei der Datenerhebung gab. Hierzu fügen wir zunächst Monat und Tag zu einem Objekt für (Kalender-)Daten zusammen:

airquality <- airquality %>%
  mutate(date = as.Date(paste("1973", Month, Day, sep = "-")))
airquality

ggplot(airquality) + 
  geom_line(aes(x = date, y = Ozone), color = "tomato1")

In den meisten Fällen scheinen die fehlenden Einträge nicht mehrere aufeinanderfolgende Tage zu beschreiben. Im Juni scheinen gänzlich am Anfang und Ende des Monats die Einträge zu fehlen.

Es könnte auch interessant sein, den Verlauf der Ozonwerte an den verschiedenen Monatstagen zu vergleichen oder zwischen den Monat (oder beides zugleich):

geom_histogram(aes(x = Ozone, fill = Day)) +

facet_wrap(~ Day)

ggplot(airquality) +
  geom_boxplot(aes(Day, Ozone, fill = Day))

ggplot(airquality) +
  geom_histogram(aes(x = Ozone, fill = Day)) +
  facet_wrap(~ Day)

ggplot(airquality) +
  geom_boxplot(aes(Month, Ozone, fill = Month))

ggplot(airquality, aes(x = Ozone, fill = Month)) +
  geom_histogram(aes(y = after_stat(density))) +
  geom_density(alpha = 0.2) +
  facet_wrap(~ Month) +
    theme(legend.position = NULL)

Bei der Betrachtung der Monatstage sehen wir tendenziell am Anfang und Ende des Monats mehr Variabilität (wohlmöglich wegen der fehlenden Werte im Juni). In der Mitte des Monats scheinen die Werte weniger auseinanderzugehen, auch wenn es Ausreißer gibt.

Bei der Betrachtung der Monate sehen wir, dass es durchschnittlich vermutlich mehr Ozon im Juli und August gab, wobei hier auch die Varianz höher war. Besonders im September gab es viele Ausreißer. In dem Histogramm vom Juni sehen wir insbesondere die fehlenden Werte. Auch hier werden die Muster der Ozonwerte etwas deutlicher, dass im Mai und September (und ggf. Juni) tendenziell (gehäuft) geringere Werte vorliegen.

Der Vorteil an rstatix liegt auch darin, dass sich verschiedene Funktionen auf Gruppen berechnen lassen. So können wir auch separat für jeden Tag/Monat die Ausreißer von Ozone bestimmen:

Wir hatten auch bereits gelernt, wie wir für Gruppen/Kategorien statistische Zusammenfassungen machen können, um die Sachverhalte auch durch Zahlenwerte betrachten zu können. Hierbei muss auf die NAs geachtet werden:

Meist ist es hilfreich, bei interessanten Aspekten zusätzlich zu einer Grafik auch markante statistische Maßzahlen zu benennen. Beispielsweise sind die Standardabweichungen und Mediane/Mittelwerte von Juli und August hoch (hohe Varianz). Somit kann auch die entsprechende Maßzahl angegeben werden.

Betrachten wir letztlich noch, ob die Variable Ozone normalverteilt ist:

airquality %>%
  shapiro_test(Ozone)

• Da der p-Wert wesentlich kleiner ist als übliche Signifikanzniveaus ⍺, lehnen wir die Nullhypothese ab, d. h. wir gehen davon aus, dass die Daten nicht normalverteilt sind.
• Es stellt sich dennoch die Frage, ob die Daten der einzelnen Monate normalverteilt sind:

airquality %>%
  group_by(Month) %>%
  shapiro_test(Ozone)

• Scheinbar gilt p < ⍺ für Mai und September => keine Normalverteilung
• Für Juni, Juli und August sind die Daten (vermutlich) normalverteilt, da der p-Wert kleiner als ⍺ ist.
• Wenn wir die anderen Variablen in die Funktion shapiro_test() eingeben würden, könnten wir direkt alle gemeinsam auf Normalverteilung testen und müssten dies nicht einzeln machen.

Entsprechende Überlegungen könnten auch für die anderen Variablen in Betracht gezogen werden.

1.2.4 Finden Sie Zusammenhänge zwischen den verschiedenen Variablen.

1.2.4.1 Visualisieren Sie für alle sinnvollen Variablen die paarweisen Beziehungen. Sind diese alle linear?

Eine Herangehensweise ist, für jedes Paar an Variablen ein neues Streudiagramm zu zeichnen:

ggplot(airquality, aes(x = Ozone, y = Temp)) +
    geom_point(aes(color = Month)) +
    geom_smooth(color = "red")

• Hier sehen wir beispielsweise einen scheinbar quadratischen Zusammenhang (Parabel)
• Das ist jedoch aufwendig, wenn kein spezifisches Interesse an jedem Variablenpaar besteht

Eine Standardfunktion in R, die einen Überblick über einen gesamten Dataframe liefert, ist die plot() Funktion:

plot(airquality)

• Die Grafik ist jedoch nicht sehr ordentlich
• Weitere Alternativen sind pairs() oder ggpairs() aus dem Paket GGally:

pairs(airquality, panel = panel.smooth) # panel = panel.smooth zeichnet die Kurven ein

library(GGally)
airquality %>%
    select(-c(Day, date)) %>%
ggpairs(progress = FALSE)

• Beachten Sie, dass für jede Grafik die NAs entfernt wurden.
• Diese Grafik ist in vielerlei Hinsicht hilfreich:
  • Zunächst sehen wir auf der Diagonalen die Verteilung der Variablen
  • Oberhalb der Diagonalen sehen wir die Korrelationskoeffizienten für die Variablenpaare oder Boxplots nach Gruppen arrangiert
  • Unterhalb der Diagonalen sehen wir die Streudiagramme für die Variablenpaare oder Histogramme nach der Gruppe arrangiert
  • Die Grafiken sind im Gegensatz zu den obigen Grafiken nicht “doppelt” (x- und y-Achse)
  • Die Grafik lässt sich noch weiter verfeinern
  • Durch die Vielzahl an Grafiken und dem hohen Informationsgehalt ist die Grafik gut geeignet, um sich eine Übersicht zu schaffen
  • Die Grafik mag in dieser Ausführlichkeit jedoch nicht unbedingt die erste Wahl sein, wenn Daten/Ergebnisse vorgestellt werden sollen (es sei denn, es gibt Platzmangel …)

Fokussieren wir uns nur auf die kontinuierlichen Daten und erweitern die Grafik etwas:

# In columns lassen sich ebenfalls die Variablen auswählen:
airquality %>%
    ggpairs(columns = 1:4, aes(color = Month, alpha = 0.01), progress = FALSE)

• Durch Aesthetics in ggpairs() können wir ebenfalls alle Grafiken gruppieren (hier nach dem Monat)
• Zusätzlich erhalten wir für jede Gruppe separat den Korrelationskoeffizienten
• Weitere Information finden Sie hier.

1.2.4.2 Berechnen Sie für alle sinnvollen Variablen eine Korrelationsmatrix. Bietet sich Pearsons Korrelationskoeffizient an?

Wir hatten zuvor die Korrelationskoeffizienten schon gesehen, wollen das an dieser Stelle noch einmal gesamteinheitlich machen:

airquality %>%
    select(-c(Day, Month, date)) %>% # Faktoren/Gruppen nicht beachten
    cor_mat()

• Temp und Ozone sind stark positiv korreliert
• Wind und Ozone sowie Temp und Wind sind moderat negativ korreliert

1.2.4.3 Bestimmen Sie durch Visualisierung, wo Zusammenhänge vorliegen.

airquality %>%
    select(-c(Day, Month, date)) %>% # Faktoren/Gruppen nicht beachten
    cor_mat() %>%
    cor_reorder() %>%
    cor_plot(type = "upper", label = TRUE,
                     palette = ggpubr::get_palette(c("darkred", "cornsilk", "darkcyan"), 200)
                     )

Die Zusammenhänge wurden oben bereits erläutert. Zusätzlich haben Wind und Solar.R keinen Zusammenhang und die restlichen Paare nur schwache Zusammenhänge.

1.2.5 Erstellen Sie für relevante (kausale) Zusammenhänge Regressionsmodelle.

Inbesondere Wind und Ozone sowie Temp und Ozone könnten interessante Zusammenhänge sein:

• Wind kann das das Ozongas wegtransportieren, sodass die Anzahl der Ozonmoleküle in der Luft sinkt

• Sonneneinstrahlung und hohe Temperaturen begünstigen die Reaktion von Stickstoff und Sauerstoff zu Ozon

ggplot(airquality, aes(Wind, Ozone)) +
    geom_point() +
    geom_smooth()

• Die Trendkurve zeigt, dass der Ozongehalt in der Luft mit steigendem Wind sinkt
• Der Zusammenhang könnte quadratisch sein, wohlmöglich aber auch linear

Erstellen wir das zugehörige Modell:

# Linearer Zusammenhang:
lin_modell <- lm(Ozone ~ Wind, data = airquality)
summary(lin_modell)
# Quadratischer Zusammenhang:
quad_modell <- lm(Ozone ~ Wind * I(Wind ** 2), data = airquality)
summary(quad_modell)

• Da der R-Quadrat-Wert von dem quadratischen Zusammenhang deutlich besser ist, scheint dieser zu bevorzugen

Dasselbe können wir nun auch für Temp und Ozone machen:

ggplot(airquality, aes(Temp, Ozone)) +
    geom_point() +
    geom_smooth()

• Es liegt ein steigender Trend vor
• Auch hier ist nicht ganz eindeutig, ob dieser linear ist

Erstellen wir wieder Regressionsmodelle:

# Linearer Zusammenhang:
lin_modell <- lm(Ozone ~ Temp, data = airquality)
summary(lin_modell)
# Quadratischer Zusammenhang:
quad_modell <- lm(Ozone ~ Temp * I(Temp ** 2), data = airquality)
summary(quad_modell)

• Das Modell mit dem quadratischen Koeffizienten ist nicht deutlich besser als das Modell mit dem linearen Koeffizienten
• Dennoch sehen wir auch hier, dass sich ein großer Teil der Varianz durch die Temperatur erklären kann
• Anschaulich bedeutet dies, dass wir einigermaßen gut vorhersagen können, welcher Ozongehalt vorliegt, wenn wir lediglich die Temperatur oder die Windgeschwindigkeit kennen

1.2.5.1 Zusätzliches/Ausschau

• Angenommen, Sie sagen, es liegen 80 Fahrenheit (~26,6 Grad Celsius) vor, dann könnte ich Ihnen einigermaßen akurat sagen, wie hoch der Ozongehalt ungefähr sein muss (in Roosevelt Island):

predict(quad_modell, list(Temp = 80))

Vergleichen wir diesen Wert (42) mit den tatsächlich gemessenen Werten:

# Werte nahe 80 Fahrenheit:
nahe_80f <- airquality %>%
    filter(between(Temp, 75, 85) & !is.na(Ozone)) %>%
    select(Ozone, Temp)
nahe_80f %>%
    summarise(mittelwert = mean(Ozone))
ggplot(nahe_80f, aes(Temp, Ozone)) +
    geom_jitter() +
    geom_smooth(se = FALSE) +
    geom_hline(yintercept = 42, color = "red")

• Wir sehen, dass die Punkte im Mittelwert nahe bei unserer Vorhersage 42 liegen
• Die individuellen Punkte liegen allerdings gestreut um die Linie

Letztendlich können wir auch ein etwas komplexeres Modell erstellen, um den Wert von Ozone zu berechnen:

modell <- lm(Ozone ~ Wind + I(Wind ** 2) + I(Temp ** 2), data = airquality)
summary(modell)

• Das Modell wird durch jeden weiteren Parameter besser in seiner Vorhersage
  • Das kann jedoch dazu führen, dass die Vorhersage nur noch für den eigenen Datensatz relevant wird (Overfitting)
  • Testen wir einfach alle möglichen Kombinationen und der R-Quadrat-Wert wird sehr hoch sein:

modell <- lm(Ozone ~ Wind * I(Wind ** 2) * Temp * Solar.R * Month, data = airquality)
summary(modell)

• Manchmal ergibt es Sinn, einfach ein paar Zusammenhänge zu betrachten und Vermutungen anzustellen, welche Zusammenhänge vorliegen
• Innerhalb des Modells sehen wir an dem t-Test von jedem Koeffiziente, welcher Koeffizient einen signifikanten Einfluss hat

Um ein Modell visuell darzustellen, kann das Paket GGally mit den Funktionen ggnostic() und ggcoef_model() helfen:

Um sich die Koeffizienten schnell anzuschauen (anstelle der Ansicht in der R-Konsole), kann ggcoef_model() verwendet werden. Hierzu wird zunächst broom.helpers benötigt (install.packages("broom.helpers"):

library(broom.helpers)
modell <- lm(Ozone ~ Temp + Wind + Solar.R + Month, data = airquality)
summary(modell)
ggcoef_model(modell)

• Beta zeigt hierbei den standardisierten Koeffizienten und dessen Unsicherheit an

• Temp und Wind haben in dem Modell einen signifikanten Einfluss (kleiner p-Wert).

• Solar.R ist ebenfalls noch signifikant.

• Month hingegen scheint nur im September einen signifikanten Einfluss auf Ozone zu haben und weist hohe Schwankungen auf.

Um weitere Aspekte des Modells zu betrachten und um seine Güte zu beurteilen, kann danach ggnostic() verwendet werden:

ggnostic(modell, aes(color = Month), progress = FALSE)

• Die verschiedenen Bereiche auf der y-Achse stellen das Modell auf unterschiedliche Arten dar (Residueen, Diagonale der Hat-Matrix, Schätzung der SD von Residueen, wenn die Beobachtung aus dem Modell genommen wurde, Cooks-Distanz)

• Der wichtigste Aspekt liegt darin, dass Ausreißer anhand der Trennlinien erkennbar werden

• Auf der x-Achse sind die Variablen des Modells

Filtern wir die zwei Ausreißer von Ozone, um ggf. das Modell zu verbessern:

# Entfernen von hohen Ozonwerten
airquality_gefiltert <- airquality %>%
    filter(Ozone < 125)
modell <- lm(Ozone ~ Temp + Wind + Month, data = airquality_gefiltert)
summary(modell)
ggnostic(modell, aes(color = Month), progress = FALSE)

• Das Modell scheint (ein wenig) besser geworden zu sein