Regressionsmodelle visualisieren in R: Mit Interaktionseffekten, 3D (ggplot2, plotly)

Regressionsmodell: 3D-Darstellung, Ebene im Raum statt Regressionsgerade

Regressionsmodelle sind nach wie vor sehr populär in der Statistik, dem Data Mining, Data Science und Machine Learning – das belegen aktuelle Zahlen, die KDNuggets kürzlich via Twitter präsentierte:

Heute geht es um Möglichkeiten, solche Modelle mit der frei erhältlichen Software R / RStudio zu visualisieren. Wir nutzen den weit verbreiteten Datensatz mtcars, der in R integriert ist.

Modell 1: Einfache lineare Regression

Zunächst eine einfache lineare Regression. Zur Darstellung benötigen wir nicht mal ein Modell – ggplot2 übernimmt das für uns. Modelliert wird der Verbrauch von einigen alten US-Automodellen in Abhängigkeit von der PS-Zahl des Motors. Anders als in Deutschland üblich, wird der Verbrauch in Gallonen pro Meile angegeben, d. h. je höher der Wert, desto sparsamer das Auto (weil es eine größere Entfernung mit der gleichen Spritmenge zurücklegt).

Einfache lineare Regression (ggplot2)
Einfache lineare Regression (R, ggplot2)

Hier der Code dazu:

library(ggplot2)

ggplot(mtcars, aes(x = hp, y = mpg)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE, size = 0.5) +
  labs(x = "hp (PS, horsepower)", y = "mpg - Verbrauch in miles per gallon\n(Je höher, desto sparsamer)",
       title = "lm(mpg ~ hp, data = mtcars)")

Mit geom_smooth() wird die Regressionsgerade in das Streudiagramm eingefügt. „lm“ steht für lineares Modell.

Modell 2: Zwei parallele Regressionsgeraden

Nun fügen wir eine kategoriale Variable mit zwei Ausprägungen hinzu: Schaltgetriebe vs. Automatik. Wir möchten den gleichen Zusammenhang wie eben darstellen, aber separat für die beiden Autotypen.

Parallele Regressionsgeraden
Parallele Regressionsgeraden (R / ggplot2, broom)

Autos mit Schaltgetrieben sind laut dieser Darstellung sparsamer (sie schaffen mehr Meilen pro Gallone). Englische Modellbezeichnung: parallel slopes model.

Eine elegante Möglichkeit, Modellvorhersagen für Grafiken zu nutzen, bietet das broom-Paket von David Robinson, das sich bestens in Hadley Wickhams tidyverse einfügt. Man kann damit Modellergebnisse in „saubere“ (tidy) Datensätze umwandeln und einfach weiterverarbeiten, auch für Diagramme.

Hier der Code:

library(dplyr)
mtcars$am <- factor(mtcars$am)
mtcars$am <- recode(mtcars$am, "0" = "Automatik", "1" = "Schaltgetriebe")

mod2 <- lm(mpg ~ disp + am, data = mtcars)

library(broom)
ggplot(augment(mod2), aes(x = disp, y = mpg, color = am)) +
 geom_point() +
 geom_line(aes(y = .fitted), size = 1) +
 labs(x = "disp (Verdrängung / Hubraum in cubic inch)", y = "mpg - Verbrauch in miles per gallon\n(Je höher, desto sparsamer)",
 title = "lm(mpg ~ disp + am, data = mtcars)")

Der „Trick“ ist die augment-Funktion, die die Modellvorhersagen in den Datensatz aufnimmt. So können wir die Linien nach den Modellwerten einzeichnen (geom_line, .fitted). Durch die Farb-Angabe in der ersten ggplot-Zeile werden automatisch zwei Linien erstellt für die beiden Kategorien von am.

Wichtig für die Interpretation: Dass die Linien parallel verlaufen, ist eine Modellannahme und damit kein empirisches Ergebnis! Wir wollen nun prüfen, ob diese Modellannahme gerechtfertigt ist.

Verlaufen die Geraden wirklich parallel?

Dazu zeichnen wir die Regressionsgeraden separat nach den beiden Autotypen ein, wie oben mit geom_smooth:

Regressionsmodell mit zwei sich schneidenden Geraden
Regressionsmodell mit zwei sich schneidenden Geraden

Diese Überprüfung hat sich gelohnt: Die Geraden verlaufen offenbar nicht parallel. Mit zunehmendem Hubraum fällt bei Autos mit Schaltgetriebe die Reichweite schneller als bei Automatik-Autos. (In anderen Worten: Der Verbrauch steigt bei Autos mit Schaltgetriebe schneller an.)

Hier wieder der Code:

ggplot(mtcars, aes(x = disp, y = mpg, color = am)) +
 geom_point() +
 geom_smooth(method = "lm", se = FALSE) +
 labs(x = "disp (Verdrängung / Hubraum)", y = "mpg - Verbrauch in miles per gallon\n(Je höher, desto sparsamer)",
 title = "lm(mpg ~ disp * am, data = mtcars)")

Welches Regressionsmodell kann diesen Zusammenhang abbilden? Sich schneidende bzw. nicht parallele Regressionsgeraden verweisen auf Interaktionseffekte bzw. Moderatoreffekte. Die Getriebeart moderiert den Zusammenhang zwischen Hubraum und Verbrauch.

Modell 3: Regressionsmodell mit Interaktionseffekt

In R kann man Interaktionseffekte sehr einfach modellieren, indem man die betroffenen Variablen direkt in der Modellformel multipliziert, hier: disp * am. R bildet dann ein Modell, das automatisch die beiden Haupteffekte und den Interaktionseffekt enthält. (Mit disp:am könnte man nur den Interaktionseffekt abbilden.)

Ist dieser Interaktionseffekt statistisch signifikant?

mod3 <- lm(mpg ~ disp * am, data = mtcars)
summary(mod3)
Regressionsmodell mit Interaktionseffekt
Regressionsmodell mit Interaktionseffekt

Ja, da ist er: p = 0,01 (disp:amSchaltgetriebe).

Haben wir dieses Modell mit der obigen Darstellung korrekt wiedergegeben? Zur Kontrolle verwenden wir einen Code, der nicht die lm-Funktion des ggplot2-Befehls nutzt, sondern die Modellwerte einsetzt. Ähnlich zu oben greifen wir wieder auf die augment-Funktion des broom-Pakets zurück:

ggplot(augment(mod3), aes(x = disp, y = mpg, color = am)) +
 geom_point() +
 geom_line(aes(y = .fitted), size = 1) +
 labs(x = "disp (Verdrängung / Hubraum in cubic inch)",
 y = "mpg - Verbrauch in miles per gallon\n(Je höher, desto sparsamer)",
 title = "lm(mpg ~ disp * am, data = mtcars)")

Tatsächlich erhalten wir das gleiche Diagramm.

Seit dem Umstieg auf R verzichte ich gern auf Excel-Tools, um Interaktionseffekte zu visualisieren.

Die dritte Dimension: Zwei metrische Prädiktoren – die Gerade wird zur Ebene

Was passiert, wenn wir zwei metrische Prädiktoren verwenden, hier z. B. hp (PS) und disp (Hubraum)? Dann begeben wir uns in die dritte Dimension, aus der Regressionsgeraden wird eine Ebene, eine Fläche im Raum. Das ist schwierig darzustellen, aber zum Beispiel mit dem plotly-Paket möglich. Hier als statisches Bild:

Regressionsmodell: 3D-Darstellung, Ebene im Raum statt Regressionsgerade
Regressionsmodell: 3D-Darstellung, Ebene im Raum statt Regressionsgerade (R, plotly)
lm(mpg ~ hp + disp, data = mtcars)
(Klicken für größere Darstellung)

Die Erstellung ist etwas aufwändiger, da man eine Matrix mit Vorhersagewerten berechnen muss, die dann die Ebene darstellt. Hier der Code fürs Diagramm:

library(plotly)
p <- plot_ly(data = mtcars, z = ~mpg, x = ~disp, y = ~hp, opacity = 0.6) %>% add_markers()
p %>% add_surface(z = ~plane, x = ~disp, y = ~hp, showscale = FALSE) %>% layout(showlegend = FALSE)

Im Browser kann man solche Diagramme sogar interaktiv darstellen, d. h. man kann es drehen und die Datenpunkte aus verschiedenen Blickwinkeln sehen. Ich bin etwas skeptisch, was die Lesbarkeit solcher Darstellungen betrifft: Dreidimensionale Grafiken auf zweidimensionalen Oberflächen (Bildschirm, Papier) stellen einen Kompromiß dar mit der Gefahr der Fehl-Interpretation. Nützlich finde ich die Darstellung, um verständlicher zu machen, was in multiplen Regressionsmodellen passiert (ohne dass man aus dem Diagramm bestimmte Messwerte genau ablesen muss).

Diagnostische Plots / Regressions-Diagnostik

An dieser Stelle kann sich der Forscher wie ein Arzt fühlen: Es gilt, das erstellte Modell zu diagnostizieren. In Base R geht das nahezu unschlagbar einfach. plot(mod3) genügt – ich habe lediglich zwei Zeilen hinzugefügt, um die vier Diagramme gemeinsam darzustellen.

par(mfrow = c(2, 2))
plot(mod3)
par(mfrow = c(1, 1))

Ergebnis:

Regressions-Diagnostik: Base R
Regressions-Diagnostik: Base R

Eleganter ist es, auch hier auf ggplot2 zurückzugreifen. Dabei unterstützt uns das ggfortify-Paket von Masaaki Horikoshi und Yuan Tang und macht uns die Arbeit sehr leicht:

library(ggfortify)
autoplot(mod3)

Ergebnis:

Regressionsdiagnostik mit ggplot2 / ggfortify
Regressionsdiagnostik mit ggplot2 / ggfortify

Natürlich sind noch weitere Diagramme möglich, z. B. vorhergesagte Werte vs. tatsächliche Werte.

R-Schulungen

Visualisierung statistischer Daten

Buchempfehlungen:


R for Data Science

2 Gedanken zu „Regressionsmodelle visualisieren in R: Mit Interaktionseffekten, 3D (ggplot2, plotly)“

Freue mich über Kommentare!