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 Meilen pro Gallone 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:

mod3 <- lm(mpg ~ hp + disp, data = mtcars)

hp <- mtcars$hp
disp <- mtcars$disp

grid <- expand.grid(hp, disp)
d <- setNames(data.frame(grid), c("hp", "disp"))
vals <- predict(mod3, newdata = d)

mpg <- matrix(vals, nrow = length(d$hp), ncol = length(d$disp))
plane <- mpg

rm(d, grid, vals)

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.

Für eine ausführlichere Diskussion über die Auswahl der geeigneten statistischen Methode siehe den Beitrag Methodenberatung: Welcher statistische Test passt zu meiner Fragestellung und meinen Daten?

R-Schulungen

Buchempfehlungen (Hinweis: bezahlte Links):

ggplot2: Elegant Graphics for Data Analysis (Use R!)

R for Data Science

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

  1. Hallo Herr Riepl,

    vielen Dank für den Beitrag, war mir eine große Hilfe! Jetzt weiß ich, dass ich eine signifikante Interaktion habe.

    Was ich nur noch nicht verstehe ist, ob mein Moderator den Zusammenhang zwischen UV und AV verstärkt und schwächt. Hierzu habe ich lediglich nur die Werte über R mit dem Befehl aus ihrem Beispiel:

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

    Muss ich da das Vorzeichen beim Interaktionsterm in der letzten Zeile bei Estimate überprüfen? bei Minus = schwächt es, bei plus = verstärkt es? Oder muss ich da die Vorzeichen mit anderen Koeffizienten vergleichen??

    Wäre Ihnen dankbar, wenn Sie mir hierbei helfen könnten. Durch die grafische Visualisierung kann ich es leider nicht feststellen.

    LG Zara

    1. Hallo Zara,
      mit den Vorzeichen tue ich mich auch schwer. Es geht nicht um das Vorzeichen des Moderators allein, sondern um die Kombination der drei Vorzeichen (im Beispiel: disp, am, Interaktionsterm). „Verstärken“ trifft es nicht immer. Im mtcars-Beispiel hat am nur zwei Ausprägungen. Bei am = 1 (manuelle Schaltung) fällt mpg mit steigendem disp stärker als bei am = 0 (Automatik). Ich würde also die beiden Geraden vergleichend beschreiben anstatt zu überlegen, ob ein Effekt verstärkt wird oder nicht.

      1. Vielen Dank für Ihre schnelle Antwort!

        Leider ist mir die grafische Darstellung nicht möglich, da mir ständig ein Fehler angezeigt wird. Meine Hypothesen lautet: Je höher der Moderator, desto schwacher der negative Zusammenhang von UV und AV.
        UV und Moderator haben beide in dem o.g. Befehl minus als Vorzeichen und mein Interaktionsterm plus. Kann ich es anhand dessen sagen?

        Oder kann man es anhand des Korrelations- oder linearen Regressionskoeffizienten von UV und AV herausfinden?

        Danke und LG Zara

        1. Ich kann es so nicht sagen.
          Klingt nach kontinuierlichem (metrischem) Moderator. Interpretationshilfe könnte sein: Moderator in zwei Gruppen teilen (hoch und niedrig) und für die beiden Gruppen separate Modelle aufstellen, ohne Moderator, nur mit UV und AV, und Koeffizienten vergleichen. Vielleicht klappt es dann auch mit der Visualisierung – das finde ich immer hilfreich.

          1. Sie sind toll, vielen vielen lieben Dank!!! Mit Ihrem Tipp hat bei mir jetzt endlich die grafische Darstellung funktioniert

  2. Hallo Herr Riepl,

    vielen Dank für Ihre Beiträge und die Mühe!!

    Ich schreibe aktuelle meine Bachelorarbeit und untersuche u.a. mit einer Hypothese, inwieweit Zukunftsorienterung (Moderationsvariable) den Zusammenhang zwischen Prokrastination (UV) und Finanzverhalten (AV). Hierbei habe ich in R den Interaktionseffekt getestet und bekomme folgenden Werte raus:

    Coefficients:
    Estimate Std. Error t value Pr(>|t|)
    (Intercept) 27.332766 3.461101 7.897 7.37e-14 ***
    pr -0.580451 0.139500 -4.161 4.27e-05 ***
    zo -0.264282 0.215033 -1.229 0.220139
    pr:zo 0.032562 0.008929 3.647 0.000319 ***

    Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Residual standard error: 3.014 on 268 degrees of freedom
    Multiple R-squared: 0.2847, Adjusted R-squared: 0.2767
    F-statistic: 35.55 on 3 and 268 DF, p-value: < 2.2e-16

    Wie kann ich das interpretieren? Ich will ja testen wie stark meine Moderationsvariable "Zukunftsorientierung" den Zusammenhang zwischen "Prokrastination" und "Finanzverhalten" beeinflusst. Aber oben bekomme ich Werte für pr:zo. Wie kann ich diese auf mein Finanzverhalten beziehen? Da rein logisch das eine sinkt und das andere steigt bei einem verantwortungsbewussten Finanzverhalten. (verantw. Finanzverhalten sinkt bei Prokrastination und steigt bei Zukunftsorientierung.)

    Ich hoffe, dass Sie mir weiterhelfen können und bedanke mich vorab für Ihre Bemühungen!

    Liebe Grüße
    Dersim

    1. Es gibt eine signifikante Wechselwirkung (Interaktion). Für eine inhaltliche Interpretation empfehle ich Visualisierung.

  3. Klasse Beitrag! Wären Sie so freundlich, auch den Code für die als Objekt ,,plane“ gespeicherte Matrix bereitzustellen! Das wäre eine enorme Hilfe, um ein solches Interaktionsdiagramm mit einer Ebene zu erstellen. Meinen besten Dank im Voraus!

  4. Guten Tag Herr Riepl
    Vielen Dank für Ihre ausgezeichneten & hilfreichen Artikel. Im Abschnitt „Modell 1: Einfache lineare Regression“ sollte wahrscheinlich der Satz „Anders als in Deutschland üblich, wird der Verbrauch in Gallonen pro Meile angegeben, d. h. …“ anders lauten, nämlich: „Anders als in Deutschland üblich, wird der Verbrauch in Meilen pro Gallonen angegeben, d. h. …“. Zumindest lässt die dazugehörende Grafik bzw. deren vertikale Achse (mpg) darauf schliessen.
    Freundliche Grüsse
    Stefan Etterlin, St. Gallen

    1. Guten Tag Herr Etterlin,
      vielen Dank für den Hinweis! Sie haben natürlich Recht. Habe den Satz korrigiert.
      Viele Grüße nach St. Gallen!
      Wolf Riepl

Freue mich über Kommentare!

Wir benutzen Cookies um die Nutzerfreundlichkeit der Webseite zu verbessen. Durch Deinen Besuch stimmst Du dem zu.