penguins-simpson

lm
bayes
regression
causal
Published

October 10, 2024

1 Aufgabe 1

Laden Sie den Datensatz penguins (Palmerpenguins). Tipp: Es gibt ein R-Paket, in dem diese Daten wohnen. Im Skript QM2 finden sich diese Daten auch. Oder im weiten Internet.

1.1 Modell 1

1.1.1 Modell berechnen

Berechnen Sie ein Modell, um den Zusammenhang zwischen Schnabeltiefe (UV) und Körpergewicht (AV) statistisch zu schätzen.

1.1.2 Visualisieren Sie das Modell

Tipp 1: estimate_relation() aus dem Paket easystats, s. QM1 oder QM2. Aber es gibt auch andere Wege.

Tipp 2: Fragen Sie ChatGPT.

1.1.3 Präzision der Koeffizienten

Geben Sie die Präzision der Regressionskoeffizienten an. Interpretieren Sie das Ergebnis.

1.2 Modell 2: Gewicht als Funktion von Schnabeltiefe und von Spezies

1.2.1 Modell berechnen

Berechnen Sie folgendes Modell: Gewicht als Funktion von Schnabeltiefe und von Spezies

1.2.2 Visualisieren Sie das Modell

wie oben

1.2.3 Präzision der Koeffizienten

Vergleichen Sie die Präzision der Regressionskoeffizienten mit dem Modell 1. Interpretieren Sie das Ergebnis.













2 Lösung zu Aufgabe 1

2.1 Modell 1

Daten importieren:

data("penguins", package = "palmerpenguins")

Achtung: Das Paket muss installiert sein.

Pakete starten:

library(tidyverse)
library(easystats)

Modell berechnen mit lm, d.h. “frequentistisch”:

m1_freq <- 
  lm(body_mass_g ~ bill_depth_mm, data = penguins)

Mit Bayes:

library(rstanarm)  # Paket muss installiert sein
m1_bayes <- 
  stan_glm(body_mass_g ~ bill_depth_mm, data = penguins,
           refresh = 0)

Modellparameter:

parameters(m1_freq)
Parameter     | Coefficient |     SE |             95% CI | t(340) |      p
---------------------------------------------------------------------------
(Intercept)   |     7488.65 | 335.22 | [6829.29, 8148.01] |  22.34 | < .001
bill depth mm |     -191.64 |  19.42 | [-229.84, -153.45] |  -9.87 | < .001

Pro Millimeter Schnabeltiefe sinkt das Gewicht um knapp 200g, im Schnitt, laut Modell.

Die Null ist NICHT im Schätzbereich enthalten, also können wir die Hypothese, dass der Zusammenhang zwischen Schnabeltiefe und Gewicht 0 ist, verwerfen.

Wir entscheiden uns also zu glauben, dass es einen Zusammenhang gibt. Wir können nicht ganz sicher sein, aber das Modell befürwortet diese Entscheidung.

Allerdings sind wir nicht sicher, ob das ein Scheinzusammenhang ist oder ein “echter”, d.h. kausaler Zusammenhang.

parameters(m1_bayes)
Parameter     |  Median |             95% CI |   pd |  Rhat |     ESS |                       Prior
---------------------------------------------------------------------------------------------------
(Intercept)   | 7483.18 | [6853.18, 8102.24] | 100% | 1.000 | 4014.00 | Normal (4201.75 +- 2004.89)
bill_depth_mm | -191.16 | [-227.38, -154.91] | 100% | 1.000 | 4001.00 |    Normal (0.00 +- 1015.24)

Ein typischer Befund: Frequentistische und Bayes-Ergebnisse sind - bei genügend großen Stichproben - sehr ähnlich, was die Zahlen betrifft. Sehr unterschiedlich ist aber die Interpretation.

Bayes-Interpretation:

“Mit 95% Wahrscheinlichkeit liegt der Effekt zwischen -230g und 150g pro Millimeter Schnabeltiefe, laut dem Modell.”

Frequentistische Interpretation:

“Würde man sehhhr viele Stichproben aus der zugrundeliegenden Population ziehen und für jede Stichprobe ein 95%-Konfindenzintervall berechnen würde, dann würde in 95% der Fälle das wahre Populationsmittel in diesem Intervall liegen. In unserer konkreten Stichprobe lagen die Grenzen bei ca. -230 bis -150. Ob der wahre Wert in diesem bestimmten Intervall liegt, können wir aber nicht sagen.”

Modell visualisieren:

estimate_relation(m1_freq) |> plot()

2.2 Modell 2

Modell berechnen mit lm, d.h. “frequentistisch”:

m2_freq <- 
  lm(body_mass_g ~ bill_depth_mm + species, data = penguins)

Mit Bayes:

library(rstanarm)  # Paket muss installiert sein
m2_bayes <- 
  stan_glm(body_mass_g ~ bill_depth_mm + species, data = penguins, refresh = 0)

Modellparameter:

parameters(m2_freq)
Parameter           | Coefficient |     SE |              95% CI | t(338) |      p
----------------------------------------------------------------------------------
(Intercept)         |    -1007.28 | 323.56 | [-1643.73, -370.83] |  -3.11 | 0.002 
bill depth mm       |      256.61 |  17.56 | [  222.07,  291.16] |  14.61 | < .001
species [Chinstrap] |       13.38 |  52.95 | [  -90.77,  117.52] |   0.25 | 0.801 
species [Gentoo]    |     2238.67 |  73.68 | [ 2093.74, 2383.60] |  30.38 | < .001
parameters(m2_bayes)
Parameter        |   Median |              95% CI |     pd |  Rhat |     ESS |                       Prior
----------------------------------------------------------------------------------------------------------
(Intercept)      | -1020.57 | [-1658.20, -384.09] | 99.90% | 1.000 | 2313.00 | Normal (4201.75 +- 2004.89)
bill_depth_mm    |   257.45 | [  222.63,  291.87] |   100% | 1.000 | 2306.00 |    Normal (0.00 +- 1015.24)
speciesChinstrap |    12.68 | [  -90.53,  117.67] | 59.45% | 1.002 | 3678.00 |    Normal (0.00 +- 5015.92)
speciesGentoo    |  2239.94 | [ 2092.14, 2385.42] |   100% | 1.000 | 2669.00 |    Normal (0.00 +- 4171.63)

Äh, Moment … Jetzt ist der Zusammenhang zwischen Schnabeltiefe und Gewicht nicht mehr negativ, sondern POSITIV?! Der Effekt geht in die entgegengesetzte Richtung? Kann das sein?!

Ein Bild zur Hilfe:

m2_freq |> estimate_relation() |> plot()

Tatsächlich! Jetzt ist der Zusammenhang innerhalb jeder Gruppe (Spezies) POSITIV.

Das bedeutet: Wenn wir die Spezies berücksichtigen, dann ist der Zusammenhang zwischen Schnabeltiefe und Gewicht positiv.

Diesen Vorzeichenwechsel nennt man “Simpson-Paradox”.

2.2.1 Fazit: Welches Modell ist jetzt richtig?

Da sich die Effekte komplett widersprechen (negativ vs. positiver Zusammenhang) stellt sich die Frage: Welchem Modell - Modell 1 oder Modell 2 - glauben wir jetzt?

Die Antwort ist ein klares: Kommt drauf an. Kommt drauf an, welcher Theorie zum kausalen Zusammenhang der betreffenden Variablen wir glauben.