Statistik und Data Science - Korrelation & Regression

1 Datenauswertung

In den vorherigen Vorlesungen und Übungen bestand das Ziel darin, Datenvisualisierung, Datentransformation und Kommunikation zu erlernen. In den kommenden Vorlesungen steht nun die Auswertung der Daten im Fokus.

1.1 Herangehensweise/Hilfestellungen

1. Zunächst sollten Sie sich einen Überblick über Ihre Daten verschaffen. Stellen Sie sich dabei stets die folgenden Fragen:
   • Wie viele Variablen und Beobachtungen wurden gemessen?
   • Welche Datentypen liegen bei jeder Variablen vor? Qualitative (Kategorien) oder quantitative Daten?
   • Was ist die zugrundeliegende Fragestellung? Welches Ergebnis soll “abgeleitet” werden?
2. Betrachten Sie die einzelnen Variablen. Was beschreiben/bedeuten Sie im Zusammenhang zur Datenerhebung?
   1. Wie sind quantitative Variablen verteilt (z. B. Histogramme, Boxplots)?
      • Liegen Ausreißer vor? Warum?
      • Liegen fehlende Werte vor? Warum?
      • Warum sind die Variablen so verteilt?
   2. Wie viele Beobachtungen liegen pro Kategorie bei qualitativen Variablen vor?
      • Liegen die Kategorien (größtenteils) im selben Ausmaß vor?
      • Liegen Kategorien sehr viel häufiger/seltener vor als andere? Warum?
        • Kann einen direkten Einfluss auf statistische Auswertungen haben
   3. Betrachten Sie hierfür auch statistische Maßzahlen wie z. B. Mittelwerte oder Standardabweichungen, um Ihre Daten zu beschreiben.
3. Liegen Zusammenhänge vor?
   • Quantitative Variablen nach Kategorien aufteilen (z. B. durch Gitter/Farbhervorhebung in Grafiken)
     • Unterscheiden sich die unterschiedlichen Kategorien bezüglich der quantitativen Variablen? Warum?
   • Zusammenhänge von zwei (oder mehr) quantitativen Variablen
     • Gibt es einen steigenden/fallenden Trend oder keinen Zusammenhang? Warum?
       • Streudiagramm (ggf. pro Kategorie)
       • Korrelationskoeffizient
   • Zusammenhänge von qualitativen Variablen
     • Liegen manche Kategorien gemeinsam häufiger/seltener vor als andere? Warum?
       • Häufigkeitstabelle/Häufigkeitsgrafik
       • Gruppierte Balkendiagramme
4. Suchen Sie unter den vorherigen Punkten untypische Daten/Aspekte und interpretieren Sie diese.
5. Finden Sie für Ihre Fragestellungen die passende Auswertungstechnik:
   • ANOVA bei dem Einfluss von qualitativen Variablen auf quantitative Variablen
   • Regression bei dem Einfluss von quantitativen Variablen auf eine andere quantitative Variable

In manchen Fällen liegen auch interessante Beobachtungen vor, die zunächst nichts mit der ursprünglichen Fragestellung zutun haben. Es kann sich lohnen, nicht nur zielstrebig das gewünschte Ergebnis auszuwerten, sondern auch beiläufige Variablen/Aspekte zu betrachten. Es gibt natürlich noch weitere statistische Auswertungsmöglichkeiten. Allerdings lassen sich die meisten Experimente anhand von Varianzanalyse (ANOVA) oder Regression auswerten und sind daher in dem Rahmen des Kurses am sinnvollsten.

1.2 Kontinuierliche Daten

Bei der Datenanalyse-und auswertung wird zwischen qualitativen und quantitativen Daten unterschieden. Diese Unterscheidung ist essentiell bei der Entscheidung wie mit den Daten umgegangen werden kann. In dem Rahmen des Kurses haben Sie bereits einige Grafikdarstellungen kennengelernt und nun sollen diese verwendet werden, um Daten und die zugrundeliegenden Resultate auszuwerten.

1.2.1 Verteilung und Interpretation

Betrachten wir einen Datensatz über Kirschbäume und dessen Maße:

library(tidyverse)
trees <- as_tibble(trees)
trees

Der Datensatz beschreibt die Höhe, den Durchmesser und das Volumen von 31 gefällten Kirschbäumen.

Wir erhalten also bereits durch das Einsehen der ersten Spalten, worum es sich um unseren Datensatz handelt und wie viele Beobachtungen und Variablen es gibt. Ebenfalls sehen wir, dass alle Variablen quantitativ sind, da sie kontinuierliche Messwerte beschreiben.

Durch den Befehl ?trees sehen wir auch, dass der Durchmesser in Zoll, die Höhe in Fuß und das Volumen in Kubikfuß angegeben sind. Bereiten wir den Datensatz also zunächst so auf, dass ein Nichtamerikaner die (metrischen) Einheiten lesen kann:

trees <- trees %>%
    mutate(Girth = Girth * 2.54,
                 Height = Height * 30.48 / 100,
                 Volume = Volume * 35.315)
trees

Nun liegt der Durchmeser in Zentimeter, die Höhe in Meter und das Volumen in Kubikmeter vor.

Schauen wir uns nun die üblichen statistischen Maßzahlen

summary(trees)

     Girth           Height          Volume      
 Min.   :21.08   Min.   :19.20   Min.   : 360.2  
 1st Qu.:28.07   1st Qu.:21.95   1st Qu.: 685.1  
 Median :32.77   Median :23.16   Median : 854.6  
 Mean   :33.65   Mean   :23.16   Mean   :1065.5  
 3rd Qu.:38.73   3rd Qu.:24.38   3rd Qu.:1317.2  
 Max.   :52.32   Max.   :26.52   Max.   :2719.3  

und die Verteilungen der Variablen an

library(patchwork)
histo <- ggplot(trees) +
    geom_histogram(aes(x = Girth), binwidth = 3, fill = "steelblue", color = "black") +
    labs(x = NULL)
boxpl <- ggplot(trees) +
    geom_boxplot(aes(x = Girth), fill = "steelblue", color = "black") + 
    labs(x = "Durchmesser (cm)") +
    scale_y_continuous(breaks = NULL)
histo + boxpl + plot_layout(1, 2)

histo <- ggplot(trees) +
    geom_histogram(aes(x = Height), binwidth = 1, fill = "steelblue", color = "black") +
    labs(x = NULL)
boxpl <- ggplot(trees) +
    geom_boxplot(aes(x = Height), fill = "steelblue", color = "black") + 
    labs(x = "Höhe (m)") +
    scale_y_continuous(breaks = NULL)
histo + boxpl + plot_layout(1, 2) # funktioniert nur mit der patchwork Bibliothek

Hinweis: Bei der Klassenbreite (oder Anzahl der Klassen) sollten ein paar verschiedene Parameter ausprobiert werden (bins oder binwidth), um zu sehen, welcher Wert am sinnvollsten ist.

• Je nach Fragestellung kann die generelle Information über die Höhe und den Durchmesser von Kirschbäumen bereits interessant sein, beispielsweise wenn man herausfinden will, was man erwarten kann wie hoch ein Kirschbaum (einer bestimmten Gattung) wächst.
• Ohne jegliche weitere Information über die Kirschbäume zu haben, sehen wir bereits, dass diese Bäume sehr groß zu sein scheinen (Mittelwert 23,16m) und wohl nicht in dem Garten eines Kleingärtners stehen würden.
• Die Höhe scheint normalverteilt zu sein (Histogramm ähnelt Glockenkurve)
• Die Verteilung des Durchmessers ist nur schwer zu deuten
• Beide Variablen weisen scheinbar keine Ausreißer auf

Es gibt zusätzlich oder ergänzend zum Histogramm noch die Möglichkeit die Verteilung durch Linien zu verdeutlichen. Für die absolute Häufigkeit lässt sich mit geom_freqpoly eine grobe Linie einzeichnen:

ggplot(trees, aes(x = Girth)) +
    geom_histogram(binwidth = 5, fill = "steelblue") +
    geom_freqpoly(binwidth = 5) +
    scale_y_continuous(breaks = seq(0, 10, 2)) # sonst stehen ungerade Werte auf der y-Achse

• bins oder binwidth muss bei beiden Funktionen gleich sein
• Wir sehen eine etwas verdeutlichte Auf-und Abbewegung
• Je nach Klassenbreite ist dieser Trend weniger gut einsichtig

Um eine glattere Kurve zu betrachten, werden häufig Dichtelinien in ein Dichtehistogramm eingezeichnet. Bei dem Dichtehistogramm stehen hierbei die relativen Häufigkeiten geteilt durch Klassenbreite auf der y-Achse. Dichte lässt sich durch ..density.. als y-Wert in die Aesthetics einsetzen. Anstelle von geom_freqpoly verwenden wir für die Linie nun geom_density

ggplot(trees, aes(x = Girth, y = ..density..)) +
    geom_histogram(bins = 8, fill = "steelblue", color = "black") +
    geom_density(alpha = 0.25, fill = "grey", color = "grey")

• wohlmöglich sind die Durchmesser doch normalverteilt (oder Gamma oder Chi-Quadrat …)?

1.2.2 Daten auf Normalverteilung testen

Unsere Intuition mag uns in manchen Fällen helfen abzuschätzen, welcher Verteilung unsere Daten folgen. Um sicher zu gehen, gibt es hierfür statistische Tests. Bei der Normalverteilung wird beispielsweise der Shapiro-Wilk-Test verwendet.

Nun stellt sich generell die Frage, warum man überhaupt Verteilungen betrachten sollte. Natürliche Beobachtungen und Prozesse folgen häufig Normalverteilungen, d. h. es gibt einen “Mittelpunkt”, um den sich die Beobachtungen verteilen und extreme Werte weit ober-und unterhalb sind eher unwahrscheinlich und werden selten beobachtet. Normalverteilungen sind für viele statistische Tests eine Voraussetzung. Wenn diese verletzt wird, dann verliert der Test an Aussagekraft.

Zudem kann es interessant sein, andere Verteilungen aufzufinden wie z. B. die bimodale Verteilung und zu erkunden, warum es mehrere “Peaks” gibt:

In R gibt es die Funktion shapiro.test(), um den Shapiro-Wilk-Test für die Normalverteilung durchzuführen:

• Nullhypothese (H0): Die Daten sind normalverteilt

• Alternativhypothese (H1): Die Daten sind nicht normalverteilt

• Der Shapiro-Wilk-Test, wie auch jeder statistische Test, liefert eine Statistik sowie einen p-Wert.

Üblicherweise wird ein Signifikanzniveau ⍺ (= Schwellenwert) gesetzt, mit dem wir entscheiden, welche Hypothese angenommen werden soll. Häufig wird ein Niveau von 0.05 verwendet. Das p steht hierbei für “probability”, d. h. dieser Wert gibt uns an, wie wahrscheinlich es ist, dass ein extremerer Wert als die Teststatistik betrachtet wird.

In der Praxis bedeutet das, dass wir aus 20 durchgeführten Tests erwarten, dass ein Test ein falsches Ergebnis liefert, wenn wir ein Signifikanzniveau von 0.05 auswählen. Strengere Schwellenwerte führen also dazu, dass wir mehr Sicherheit haben, dass das Ergebnis tatsächlich wahr ist. Jedoch verlieren wir bei sehr strengen Schwellenwerte (z. B. 0.00001) wohlmöglich wahre Ergebnisse. Daher sind p-Werte ein sehr umstrittenes Thema.

Achten Sie immer darauf:

• Gilt p >= ⍺ => wir lehnen H0 nicht ab

• Gilt p < ⍺ => wir lehnen H0 ab und nehmen H1 an

Führen wir den Shapiro-Wilk-Test nun durch und wählen ⍺ = 0.05:

# Spalte aus Dataframe entnehmen:
trees$Girth

 [1] 21.082 21.844 22.352 26.670 27.178 27.432 27.940 27.940 28.194 28.448 28.702 28.956 28.956 29.718 30.480 32.766 32.766 33.782 34.798 35.052 35.560 36.068 36.830 40.640
[25] 41.402 43.942 44.450 45.466 45.720 45.720 52.324

# Shapiro-Test mit dieser Spalte:
shapiro.test(trees$Girth)

    Shapiro-Wilk normality test

data:  trees$Girth
W = 0.94117, p-value = 0.08893

Der Test nimmt einen Vektor (Variable) entgegen, d. h. wir müssen aus trees die Variable Girth mit dem $-Operator entnehmen.

• Es gilt (p-Wert =) 0.089 > 0.05 (= ⍺) => H0 wird nicht abgelehnt und wir nehmen an, dass die Daten normalverteilt sind.

• Dasselbe könnten wir nun für die anderen Variablen durchführen

Um das Testen mit tidyverse etwas zu vereinfachen, wollen wir das Paket rstatix verwenden:

install.packages("rstatix")

Mit rstatix können wir Tests in R durchführen und gleichzeitig weiter mit Dataframes arbeiten. Für den Shapiro-Wilk-Test gibt es hier die Funktion shapiro_test(), dem wir einfach die Variablen übergeben können:

library(rstatix)

Attache Paket: ‘rstatix’

Das folgende Objekt ist maskiert ‘package:stats’:

    filter

trees %>%
    shapiro_test(Girth, Height, Volume)

Der Vorteil hierbei ist, dass wir weiter mit Dataframes arbeiten, d. h. wir können die Daten auch direkt filtern (z. B. nach ⍺):

nicht_sig <- trees %>% 
    shapiro_test(Girth, Height, Volume) %>%
    filter(p >= 0.05)
nicht_sig

• Wir sehen direkt, dass Durchmesser und Höhe normalverteilt sind und Volumen scheinbar nicht.

• Insbesondere lässt sich der Prozess durch beispielsweise Gruppierung erweitern (dazu später mehr).

1.2.3 Zusammenhänge

Um Zusammenhänge von kontinuierlichen Variablen zu betrachten, werden Streudiagramme oder ähnliche Grafiken verwendet.
Falls verschiedene Kategorien/Gruppen vorhanden sind, kann dies ebenfalls durch Farbhervorhebung oder durch Gitter ergänzt werden.

Bei unserem Datensatz über Kirschbäume ist intuitiv klar, wie sich Durchmesser und Höhe verhalten. Oder?

ggplot(trees) +
    geom_point(aes(Girth, Height)) +
    labs(x = "Durchmesser (cm)",
             y = "Höhe (m)", 
             title = "Kirschbäume mit höherem Durchmesser sind tendenziell größer",
             caption = "Daten aus dem Buch 'The Minitab Student Handbook'")

• Die Tendenz ist eindeutig: steigt der Durchmesser, steigt ebenfalls auch die Höhe.
• Diese Aussage mag nicht für jeden einzelnen Punkt zutreffen. Wichtig ist hierbei, dass eine Tendenz zu sehen ist.
• Dieser Zusammenhang könnte linear sein: legen wir eine passende Gerade durch die Punkte, werden die meisten Punkte mit wenig Abweichung gut getroffen.

Zuvor hatten wir durch Streudiagramme mit geom_smooth eine geglättete Linie gelegt:

ggplot(trees, aes(Girth, Height)) +
    geom_point() +
    geom_smooth() +
    labs(x = "Durchmesser (cm)",
             y = "Höhe (m)", 
             title = "Kirschbäume mit höherem Durchmesser sind tendenziell größer",
             caption = "Daten aus dem Buch 'The Minitab Student Handbook'")

• Im Hintergrund von geom_smooth wird ein Modell erstellt, welches die Datenpunkte durch diese geglättete Linie annähert.
• Je nach Methode (Modellart) können das glatte komplexe Polynome sein oder (einfache) Geraden.
• Hierzu muss der method Parameter in geom_smooth angepasst werden (?geom_smooth).

Ein Gerade lässt sich beispielsweise durch "lm" einzeichnen:

ggplot(trees, aes(Girth, Height)) +
    geom_point() +
    geom_smooth(method = "lm", color = "red") +
    labs(x = "Durchmesser (cm)",
             y = "Höhe (m)", 
             title = "Kirschbäume mit höherem Durchmesser sind tendenziell größer",
             caption = "Daten aus dem Buch 'The Minitab Student Handbook'")

• Die Gerade scheint nicht (viel) schlechter zu sein als die Kurve von oben, aber scheint den Zusammenhang vielleicht nicht perfekt darzustellen
• Das ist ein typisches Phänomen. Später sehen wir, wie die Güte von diesem Modell zu bewerten ist.

1.2.3.1 Korrelationsmatrix & Korrelationsgrafik

Zusammenhänge sind die Tendenz, dass das Steigen oder Fallen einer Variable ebenfalls mit dem Steigen oder Fallen einer anderen Variable einhergeht. Um die Stärke eines Zusammenhangs zu bestimmen, wird ein Korrelationskoffizient verwendet:

• Liegt zwischen -1 (perfekt negativ) und +1 (perfekt positiv)

• Ein Wert von 0 deutet darauf, dass die Variablen keinen Zusammenhang aufweisen

• Es gibt verschiedene Korrelationskoeffiziente (Pearson, Spearman, Kendall)

  • Pearson: falls die Variablen (relativ) linear zusammenhängen und normalverteilt sind

  • Ansonsten: Kendall (ordinale Daten) oder Spearman verwenden

Wir hatten oben bereits getestet, dass Durchmesser und Höhe normalverteilt sind. Ebenfalls könnte ein linearer Zusammenhang möglich sein. Daher verwenden wir Pearsons Korrelationskoeffizient. Hierzu wird standardmäßig die cor() Funktion verwendet:

cor(trees$Girth, trees$Height)

[1] 0.5192801

• Der Korrelationskoeffizient weist eine Wert von 0,52 auf und deutet somit einen moderaten Zusammenhang an. Das entspricht dem, was vor zuvor betrachtet haben.

• Die Punkte liegen nicht exakt auf einer Geraden, sondern weichen etwas ab. Es liegt also ein tatsächlicher Zusammenhang vor.

• Die Höhe von einem Kirschbaum hängt scheinbar nicht ausschließlich von dem Durchmesser ab, denn sonst wäre der Zusammenhang perfekt.

Hinweis: den Koeffizienten können wir mit method = X ändern:

cor(trees$Girth, trees$Height, method = "pearson") # Default/Standard

[1] 0.5192801

cor(trees$Girth, trees$Height, method = "spearman")

[1] 0.4408387

cor(trees$Girth, trees$Height, method = "kendall")

[1] 0.3168641

• Die anderen Korrelationskoeffizienten weisen sogar noch etwas geringere Werte auf.

• Zu jedem Korrelationskoeffizienten gibt es auch einen zugehörigen Test, ob der Koeffizient sich signifikant von 0 unterscheidet.

• Anders formuliert: der Test gibt an, ob es einen Zusammenhang gibt. Die Testhypothesen lauten für alle Tests:

  • H0: Der Koeffizient ist gleich 0

  • H1: Der Koeffizient ist ungleich 0

Genauso wie cor() lässt sich ergänzend für diesen Test cor.test() verwenden:

cor.test(trees$Height, trees$Girth)

    Pearson's product-moment correlation

data:  trees$Height and trees$Girth
t = 3.2722, df = 29, p-value = 0.002758
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.2021327 0.7378538
sample estimates:
      cor 
0.5192801 

• Hier sehen Sie ebenfalls den Korrelationskoeffizient ganz unten
• Der p-Wert des Tests beträgt 0,0028 und ist somit kleiner als 0,05 (= typisches ⍺)
• Es folgt, dass der Zusammenhang signifikant bzw. vorhanden ist (unterschiedlich von 0)

In einem Bericht könnte man beispielsweise ein Streudiagramm mit der geglätteten Kurve darstellen und zusätzlich das Ergebnis des Korrelationskoeffizienten berichten, sofern dieser Zusammenhang interessant ist.

Es gibt für Korrelationskoeffizienten einen “alles-in-einem” Trick, um schnell Zusammenhänge darzustellen und zu erkennen, insbesondere wenn viele Variablen vorliegen. Hierzu wird eine Matrix mit den p-Werten/Korrelationskoeffizienten für jedes Paar der Variablen erstellt und zusätzlich als Grafik veranschaulicht.

Ebenfalls mit dem cor() Befehl lässt sich eine Matrix für alle numerischen Variablen berechnen. Hierbei muss der Datensatz auf die numerischen Variablen reduziert werden. Nehmen wir als Beispiel den mpg Datensatz. Es gibt hierbei mehrere Möglichkeiten:

# Variablen per Hand auswählen:
numerisch <- mpg %>%
    select(displ, year, cyl, cty, hwy)
numerisch

# Variablen nach Bedingung filtern:
numerisch <- mpg %>%
    select_if(is.numeric)
numerisch

• Bei manchen Variablen muss man sich überlegen, ob es Sinn ergibt, diese als Kategorie oder Zahlenwert zu betrachten (z. B. das Jahr)

Nun kann durch cor() jeder paarweiser Korrelationskoeffizient berechnet werden:

cor(numerisch)

           displ         year        cyl         cty          hwy
displ  1.0000000  0.147842816  0.9302271 -0.79852397 -0.766020021
year   0.1478428  1.000000000  0.1222453 -0.03723229  0.002157643
cyl    0.9302271  0.122245347  1.0000000 -0.80577141 -0.761912354
cty   -0.7985240 -0.037232291 -0.8057714  1.00000000  0.955915914
hwy   -0.7660200  0.002157643 -0.7619124  0.95591591  1.000000000

• Die Matrix-Darstellung ist sehr hilfreich, um direkt konkrete Werte zu vergleichen
• Beispielsweise scheint der Kraftstoffverbrauch in der Stadt (cty) sehr stark mit dem Kraftstoffverbrauch auf der Autobahn (hwy) zu korrelieren, was auch sinnvoll ist
• Anschaulich ist diese Darstellung jedoch nicht

Bei unserem Kirschbaum-Datensatz war das Filtern hingegen nicht notwendig, da bereits alle Variablen numerisch sind:

cor(trees)

           Girth    Height    Volume
Girth  1.0000000 0.5192801 0.9671194
Height 0.5192801 1.0000000 0.5982497
Volume 0.9671194 0.5982497 1.0000000

• Wir sehen für jedes Paar an Variablen einen Koeffizienten (als hätten wir für alle cor() separat aufgerufen)
• Die Diagonale ist immer 1, da eine Variable stets mit sich selbst perfekt korreliert ist

Alternativ können aus dem rstatix Paket die Funktionen cor_test(), cor_mat() und cor_get_pval() verwendet werden, um die Koeffizienten bzw. die p-Werte der zugehörigen Tests darzustellen. Hierbei wird Filterung und Berechnung zusammengefasst in eine Funktion:

# Ausführlicher Korrelationstest:
trees %>%
    cor_test()

# Korrelationsmatrix
trees %>%
    cor_mat()

# p-Werte-Matrix:
trees %>%
    cor_mat() %>%
    cor_get_pval()

# konkrete Benennung:
trees %>%
    cor_mat(Girth, Height, Volume)

Um die Matrix etwas detailierte anzusehen, können wir entweder cor_test() verwenden oder cor_mat() und dann cor_gather:

trees %>%
    cor_mat() %>%
    cor_gather()

Da Grafiken eine schnellere Orientierungshilfe sind, wollen wir diese Matrix direkt darstellen. Hierzu gibt es in rstatix die cor_plot() Funktion:

mpg %>%
    cor_mat(year, displ, cyl, cty, hwy) %>%
    cor_plot()

Es lassen sich direkt einige Dinge ablesen:

• Eine negative Korrelation wird in diesem Beispiel rot markiert, eine positive Korrelation durch blau

• Je dunkler der Farbton desto stärker ist die Korrelation

• Weiß stellt hierbei Korrelationskoeffizienten bei 0 dar

• In diesem Beispiel werden zusätzlich die Koeffizienten unterschiedlich groß dargestellt (je nach Wert)

• Sofern ein p-Wert oberhalb von 0.05 liegt, wird dieser ebenfalls durch ein Kreuz markiert (= kein Zusammenhang)

Einige dieser Optionen lassen sich anpassen (?cor_plot()):

korr <- mpg %>%
    cor_mat(year, displ, cyl, cty, hwy)
# Signifikanzlevel:
korr %>%
    cor_plot(significant.level = 0.1, insignificant = "blank")

# Form
korr %>%
    cor_plot(method = "square")

# Dreieck vs Quadrat
korr %>%
    cor_plot(type = "lower")

# Farbe und Label
korr %>%
    cor_plot(label = TRUE, palette = ggpubr::get_palette(c("tomato3", "white", "seagreen"), 200))

# Umsortieren:
korr %>%
    cor_reorder() %>%
    cor_plot()

# Kombination aus allem
korr %>%
    cor_reorder() %>%
    cor_plot(method = "ellipse", significant.level = 0.1, insignificant = "blank", type = "upper",
                     label = TRUE, palette = ggpubr::get_palette(c("darkred", "cornsilk", "darkcyan"), 200)
                     )

Für weitere Optionen bei der Arbeit mit Korrelationsgrafiken bietet sich ebenfalls das ggcorrplot Paket an. Hier lassen sich noch zusätzlich Einstellungen ändern mit dessen ggcorrplot() Funktion. Siehe hier für mehr Information zu ggcorrplot.

1.2.3.2 Regressionsmodelle

Bei der Korrelationsanalyse betrachten wir Zusammenhänge. Zusammenhänge müssen jedoch nicht kausal sein, d. h. es mag sein, dass eine Variable tendenziell mit einer anderen Variable steigt, jedoch nicht durch den Einfluss der anderen Variable.

Ein solcher Sachverhalt kann komplex sein, aber häufig steckt dahinter eine dritte (oder mehrere) Variable, die beide Variablen zugleich beinflusst. An dieser Stelle wollen wir das übliche Beispiel betrachten, um diesen Unterschied zu erläutern:

Mit einer erhöhten Anzahl an Störchen in einem beobachteten Gebiet steigt auch die Geburtenrate (tatsächliches Phänomen):

• Die Korrelation liegt hierbei vor (signifikant)

• Kausal ist dieser Zusammenhang jedoch nicht

• Eine dritte Variable (die Industrialisierung) scheint beide Variablen zu beeinflussen: auf dem Land siedeln sicher eher Familien an und in Stadtgebieten gibt es wenige Störche.

• Demnach ist der tatsächliche kausale Zusammenhang bei Industrialisierung und Geburtenrate/Anzahl Störche

Kleiner Einschub: betrachten wir einen Beispieldatensatz zu Störchen und der Geburtenrate. Wir sehen, dass die “Storchenrate” und die Geburtenrate stark positiv korreliert ist (und Anzahl von Störchen und Anzahl von Babies):

library(TeachingDemos)
stork <- as_tibble(stork)
stork

stork %>%
    cor_mat() %>%
    cor_plot()

Mit Hilfe der partiellen Korrelation können wir mögliche Störfaktoren aus der Korrelation “entfernen”:

library(ppcor)

Lade nötiges Paket: MASS

Attache Paket: ‘MASS’

Das folgende Objekt ist maskiert ‘package:rstatix’:

    select

Das folgende Objekt ist maskiert ‘package:patchwork’:

    area

Das folgende Objekt ist maskiert ‘package:dplyr’:

    select

stork %>%
    pcor() %>% # partielle Korrelation
    magrittr::extract2(1) %>% # ersten Eintrag entnehmen
    round(digits = 3) %>% # Dezimalstellen runden auf 3
    cor_plot()

• Anzahl der Störche und Anzahl der Babies sind nun moderat negativ korreliert.

• Wenn wir die Storchrate und die Geburtenrate betrachten, so liegt hier nur noch eine sehr schwache Korrelation vor.

Kehren wir zurück zur Regression.

Wichtig bei der Regressionsanalyse ist, dass wir einen kausalen Zusammenhang unterstellen. Es sollte, anders als bei dem Storchenbeispiel, ein echter kausaler Zusammenhang vorliegen.
Bei dem Datensatz der Kirschbäume kann es gut sein, dass der Durchmesser des Baums eine wichtige Rolle spielt, sodass der Baum höher wachsen kann.

In R werden lineare Regressionsmodelle typischerweise mit der lm() (= “Lineares Modell”) Funktion durchgeführt. Im einfachsten Fall, so wie wir es bereits mehrfach gesehen haben, wird hiermit eine “passende” Gerade durch unsere Punkte gezogen. Erstellen wir ein solches Modell für den Durchmesser und die Höhe der Kirschbäume:

lin_modell <- lm(Height ~ Girth, data = trees)
summary(lin_modell)

Call:
lm(formula = Height ~ Girth, data = trees)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.8349 -0.8439  0.0964  0.7537  3.0314 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 18.90714    1.33603  14.152 1.49e-14 ***
Girth        0.12652    0.03867   3.272  0.00276 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.688 on 29 degrees of freedom
Multiple R-squared:  0.2697,    Adjusted R-squared:  0.2445 
F-statistic: 10.71 on 1 and 29 DF,  p-value: 0.002758

• Hier sehen wir erneut die Verwendung von dem Tilde “~” Symbol
• Links von der Tilde steht die abhängige Variable, rechts von der Tilde steht die unabhängige Variable
• Der Effekt der unabhängigen Variablen wird auf die abhängige Variable betrachtet, z. B. welchen Einfluss der Durchmesser auf die Baumhöhe hat
• Der p-Wert ganz unten gibt an, ob das Modell signifikant ist, d. h. ob der Durchmesser einen Einfluss auf die Baumhöhe hat:

  • H0: Die unabhängigen Variablen haben keinen Einfluss

  • H1: Es gibt mindestens eine unabhängige Variable mit einem Einfluss

• Da der p-Wert relativ klein ist, gehen wir davon aus, dass der Durchmesser einen signifikanten Einfluss auf die Baumhöhe hat
• Multiple R-squared (R-Quadrat) gibt die Güte des Modells an:
  • Werte zwischen 0 (schlecht) und 1 (perfekt)
  • “Wie viel Varianz kann das Modell beschreiben” = “Wie viel weichen die Werte von der Modellkurve/Gerade ab?”
  • Da R-Quadrat = 0.27 ist, sehen/schlussfolgern wir erneut, dass der Durchmesser Einfluss auf die Baumhöhe hat, aber wohlmöglich nicht der einzige Faktor ist, der einen Einfluss hat (denn sonst wäre R-Quadrat sehr nahe bei 1)
• Unter Coefficients sehen wir die Modellgerade: Height = 18,9 + 0.127 * Girth

Wollen wir die Modellgerade in eine Grafik einzeichnen, gibt es mehrere Möglichkeiten. Wir können predict() auf das Modell anwenden, um die Vorhersagen, d. h. y-Werte, des Modells zu erhalten und diese als Linien einzeichnen. Alternativ können wir einfach das Modell in geom_smooth() erstellen lassen:

# Die Vorhersagen des Modells verwenden:
ggplot(trees, aes(Girth, Height)) +
    geom_point() +
    geom_line(aes(y = predict(lin_modell)), color = "red")

# geom_smooth verwenden:
ggplot(trees, aes(Girth, Height)) +
    geom_point() +
    geom_smooth(method = "lm", formula = y ~ x, color = "red", se = FALSE)

Betrachten wir abschließend noch den Durchmesser, die Höhe und das Volumen.

ggplot(trees, aes(x = Girth, y = Volume)) +
    geom_point() +
    geom_smooth()

ggplot(trees, aes(x = Height, y = Volume)) +
    geom_point() +
    geom_smooth()

• Wenn wir Durchmesser und Volumen betrachten, könnte der Zusammenhang wohlmöglich (leicht) quadratisch (Kurve) sein.

• Bei Höhe und Volumen ist das nicht ganz eindeutig. Hier könnte es sich aber vermutlich um einen linearen Zusammenhang handeln.

Es können auch alle drei kontinuierlichen Variablen gleichzeitig betrachtet werden. Hierzu werden Gitter erstellt, wobei eine Variable in Intervalle mit einer cut_*() Funktion aufgeteilt wird. Um beispielsweise den Zusammenhang von Volumen, Durchmesser und Höhe zu betrachten, wird hier ein normales Streudiagramm mit Durchmesser und Volumen erstellt, wobei Gitter basierend auf den Intervallen der Höhe eingeteilt werden:

ggplot(trees, aes(x = Girth, y = Volume)) +
    geom_point() +
    geom_smooth(method = "lm", se = FALSE) +
    facet_wrap(~ cut_width(Height, 3))

Abschließend noch ein Ausblick, was sich mit Regressionsmodellen noch machen lässt: Ebenfalls können wir zu Volumen und Höhe/Durchmesser Regressionsmodelle erstellen:

lin_modell2 <- lm(Volume ~ Girth, data = trees)
summary(lin_modell2)

Call:
lm(formula = Volume ~ Girth, data = trees)

Residuals:
    Min      1Q  Median      3Q     Max 
-284.83 -109.71    5.37  123.42  338.56 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -1304.658    118.840  -10.98 7.62e-12 ***
Girth          70.433      3.439   20.48  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 150.2 on 29 degrees of freedom
Multiple R-squared:  0.9353,    Adjusted R-squared:  0.9331 
F-statistic: 419.4 on 1 and 29 DF,  p-value: < 2.2e-16

lin_modell3 <- lm(Volume ~ Height, data = trees)
summary(lin_modell3)

Call:
lm(formula = Volume ~ Height, data = trees)

Residuals:
   Min     1Q Median     3Q    Max 
-751.3 -349.4 -102.2  426.2 1054.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -3076.77    1033.78  -2.976 0.005835 ** 
Height        178.82      44.48   4.021 0.000378 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 473.1 on 29 degrees of freedom
Multiple R-squared:  0.3579,    Adjusted R-squared:  0.3358 
F-statistic: 16.16 on 1 and 29 DF,  p-value: 0.0003784

Der Durchmesser scheint relativ gut zu sein, um das Volumen zu bestimmen (R-Quadrat = 0.935). Allerdings wissen wir auch, dass sich das Volumen im Allgemeinen (von zylinder-artigen Objekten) berechnen lässt durch die Formel:

pi * Radius^2 * Höhe

Hierzu muss für den Durchmesser I(X ** 2) geschrieben werden, damit der Durchmesser auch quadriert in die Formel eingeht. Statt der Multiplikation * (~ Addition und Multiplikation) wird in Modellen : geschrieben.

Betrachten wir zunächst ein Modell mit dem Durchmesser zum Quadrat:

lin_modell4 <- lm(Volume ~ I(Girth ** 2), data = trees)
summary(lin_modell4)

Call:
lm(formula = Volume ~ I(Girth^2), data = trees)

Residuals:
     Min       1Q   Median       3Q      Max 
-220.630  -90.894    9.025   75.367  247.420 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -118.4868    50.0379  -2.368   0.0248 *  
I(Girth^2)     0.9917     0.0379  26.169   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 119 on 29 degrees of freedom
Multiple R-squared:  0.9594,    Adjusted R-squared:  0.958 
F-statistic: 684.8 on 1 and 29 DF,  p-value: < 2.2e-16

ggplot(trees, aes(Girth, Volume)) +
    geom_point() +
    geom_smooth(formula = y ~ I(x ** 2), method = "lm")

Tatsächlich scheint eine quadratische Kurve gut zu passen. Darauf basierend können wir unser Modell weiter verbessern. Hierzu multiplizieren wir die Höhe mit dem Durchmesser im Quadrat:

lin_modell5 <- lm(Volume ~ Height:I(Girth ** 2), data = trees)
summary(lin_modell5)

Call:
lm(formula = Volume ~ Height:I(Girth^2), data = trees)

Residuals:
     Min       1Q   Median       3Q      Max 
-163.136  -38.853   -5.848   61.629  148.237 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -10.512549  34.027959  -0.309     0.76    
Height:I(Girth^2)   0.038151   0.001068  35.711   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 88.04 on 29 degrees of freedom
Multiple R-squared:  0.9778,    Adjusted R-squared:  0.977 
F-statistic:  1275 on 1 and 29 DF,  p-value: < 2.2e-16

Mit einem R-Quadrat von 0.978 ist das Modell sogar noch etwas besser als zuvor und scheint nun beinahe perfekt das Volumen der Kirschbäume anhand von Durchmesser und Höhe zu bestimmen, auch wenn der Zusammenhang bereits bekannt war.

Was könnte der Grund sein, dass typischerweise (trotz eindeutiger funktionaler Beziehung) kein perfekter R-Quadrat von 1 erhalten werden kann?

2 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)
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

2.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, 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.