12  Fallbeispiele

Schlüsselwörter

Statistik, Bayes, Wahrscheinlichkeit, Inferenz, R, Kausalität

12.1 Lernsteuerung

Abbildung 1.1 gibt einen Überblick zum aktuellen Standort im Modulverlauf.

Nach Absolvieren des jeweiligen Kapitels sollen folgende Lernziele erreicht sein.

Sie können…

• typische, deskriptive Forschungsfragen spezifizieren als Regression
• Forschungsfragen in Regressionsterme übersetzen
• typische Forschungsfragen auswerten

Der Stoff dieses Kapitels orientiert sich an McElreath (2020), Kap. 4.4 sowie Gelman et al. (2021), Kap. 7 und 10.

Frischen Sie Ihr Wissen in den Grundlagen der einfachen und multiplen Regression (inkl. Interaktionseffekte) auf. Dazu sind z.B. folgende Literaturstellen geeignet.

• Statistik1, Kap. “Geradenmodelle 1”¹
• Statistik1, Kap. “Geradenmodelle 2”²

In diesem Kapitel werden die üblichen R-Pakete benötigt.

library(rstanarm)   # Bayes-Modelle
library(tidyverse)  # Datenjudo
library(easystats)  # Statistik-Komfort

Wir benötigen in diesem Kapitel folgende Datensätze: kidiq, penguins. Den Datensatz kidiq importieren Sie am einfachsten aus dem R-Paket rstanarm, das Sie schon installiert haben. Alternativ können Sie die Daten hier herunterladen.³

data("kidiq", package = "rstanarm")

Außerdem benötigen wir den Datensatz penguins. Sie können den Datensatz penguins entweder über den Pfad importieren oder über das zugehörige R-Paket. Beide Möglichkeiten sind okay.

penguins_url <- paste0(
  "https://vincentarelbundock.github.io/Rdatasets/",
  "csv/palmerpenguins/penguins.csv"
)

penguins <- read.csv(penguins_url)

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

Beispiel 12.1 (Was waren noch mal die Skalenniveaus?) Um Forschungsfragen zu klassifizieren, müssen Sie wissen, was die Skalenniveaus der beteiligten AV und der UV(s) sind.⁴ \(\square\)

Beispiel 12.2 (Was war noch einmal die Interaktion?) Erklären Sie die Grundkonzepte der Interaktion (hier synonym: Moderation) im Rahmen einer Regressionsanalyse!⁵ \(\square\)

Wenn Sie die Skalenniveaus kennen, können Sie die Forschungsfrage korrekt auswerten, also das korrekte (Regressions-)Modell spezifizieren. Wir werden hier viele der typischen Forschungsfragen (aus psychologischen und ähnlichen Fragestellungen) mit Hilfe von Regressionsmodellen beantworten. Das hat den Vorteil, dass sie nicht viele verschiedene Auswertungsmethoden (t-Test, Varianzanalyse, …) lernen müssen. Außerdem ist die Regressionsanalyse (für viele Situationen) die beste Herangehensweise, da sie viele Möglichkeiten für Erweiterungen bietet. Entsprechend ist das Thema dieses Kapitels, gängige Forschungsfragen mit Hilfe der Regressionsanalyse zu untersuchen. Wenn Sie die Grundkonzepte der Regression schon kennen, wird Ihnen vieles sehr bekannt vorkommen. Natürlich würzen wir das Ganze mit einer ordentlichen Portion Post-Verteilungen aus der Bayes-Küche. Allerdings kommt auch dabei nichts Wesentliches mehr hinzu, abgesehen von ein paar Erweiterungen.

Hinweis

Die Folien zu diesem Kapitel finden Sie hier.

12.2 Taxonomie von Forschungsfragen

Wir konzentrieren uns im Folgenden auf Forschungsfragen auf Basis von Regressionsmodellen mit metrischer AV. Andere Skalenniveaus bei der AV klammern wir aus. Im Folgenden sind für die UV(s) nominale sowie metrische Skalenniveaus erlaubt. Modelle mit mehreren UV (und mehreren Stufen an UV) sind ebenfalls erlaubt.

Wir untersuchen in diesem Kapitel häufig verwendete Arten von Forschungsfragen mittels Regressionsanalysen. Für jede Variante ist zumeist ein Beispiel, die Modellformel, der Kausalgraph⁶, die Forschungsfrage sowie die Grundlagen der Auswertung dargestellt.

Dabei wird folgende Nomenklatur verwendet, um die Skalenniveaus der beteiligten Variablen einer Forschungsfrage zu benennen:

• y: metrische AV
• g: Gruppierungsvariable; nominal skalierte UV (querschnittlich)
• b: binäre UV
• x: metrische UV
• u: ungemessene (unbekannte) Variable

Übungsaufgabe 12.1 (Einstieg: Was ist Ihre Lieblingsforschungsfrage?) Welche Forschungsfrage würden Sie selber gerne untersuchen – etwa im Kontext eines entsprechenden Moduls in einem Folgesemester? Auch in Ihrer Abschlussarbeit haben Sie Gelegenheit zu eigener Forschung.

Überlegen Sie sich (in Kleingruppen) eine oder mehrere Forschungsfragen, die Sie interessieren würden! Posten Sie Ihre Antwort in das von der Lehrkraft bereitgestellte Forum. \(\square\)

12.3 y ~ b

Hintergrund:

Eine Psychologin, die im öffentlichen Dienst als Schulpsychologin arbeitet, versucht herauszufinden, warum einige Kinder intelligenter sind als andere. Dazu wurden in einer aufwändigen Studie die Intelligenz vieler Kinder gemessen. Zusätzliche wurden verschiedene Korrelate der Intelligenz erhoben, in der Hoffnung, “Risikofaktoren” für geringere Intelligenz zu entdecken.

Forschungsfrage:

Ist der mittlere IQ-Wert (kid_score) von Kindern, deren jeweilige Mutter über einen Schulabschluss (mom_hs, \(x=1\)) verfügt, höher als bei Kindern, deren jeweilige Mutter nicht über einen Schulabschluss verfügt (\(x=0\))? (ceteris paribus)⁷.

Die Modellformel zur Forschungsfrage lautet allgemein y ~ b bzw. konkret für den vorliegenden Fall kid_iq ~ mom_hs.

Theorem 12.1 drückt die Forschungsfrage formaler und als Hypothese (Behauptung) aus.

Theorem 12.1 (Hypothese für ungleiche Mittelwerte) \[H_A: \mu_{x=1} > \mu_{x=0}\quad \square\]

In Worten: “Der mittlere IQ-Wert für Kinder, deren Mütter über einen Schulabschluss verfügen ist höher als in der Gruppe von Kindern, deren Mütter über keinen Schulabschluss verfügen”. Zu beachten ist, dass sich eine Hypothese immer auf Parameterwerte bezieht, also auf die Population, nicht auf die Statistiken der Stichprobe.

Die zugehörige Nullhypothese, \(H_0\) lautet:

\[H_0: \mu_{x=1} = \mu_{x=0}\quad \square\]

Für die Praktisch-Null-Hypothese (ROPE-Hypothese) wählen wir 5 IQ-Punkte als Grenzwert für gerade noch unbedeutend.

\[H_{ROPE}: \mu_{x=1} > \mu_{x=0} + 5\quad \square\]

Die Regressionsformel zur Forschungsfrage lautet: y ~ b bzw. kid_iq ~ mom_hs.

Der Kausalgraph zur Modellformel ist in Abbildung 12.1 dargestellt. Y hat, laut unserem Modell, zwei Ursachen:

1. mom_hs (b)
2. u, das steht für “unbekannt”⁸

Abbildung 12.1: DAG für kid_iq ~ mom_hs

Auf dieser Basis berechnen wir unser Regressionsmodell mit Stan: stan_glm(kid_score ~ mom_hs, data = kidiq)

m10.1 <- stan_glm(
  kid_score ~ mom_hs, 
  data = kidiq)

Der Einfachheit halber übernehmen wir die Prioriwerte von Stan, s. Listing 12.1.

Listing 12.1: Standard-Prioriwerte für m10.1, von Stan vergeben

prior_summary(m10.1)
## Priors for model 'm10.1' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 87, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 87, scale = 51)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = 0, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 0, scale = 124)
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.049)
## ------
## See help('prior_summary.stanreg') for more details

Die komplette Modellspezifikation ist in Gleichung 12.1 aufgeführt. Wie man sieht, sind die Priori-Verteilungen sehr breit – zu breit vielleicht. Das sollten wir beim nächsten Mal verbessern.

\[\begin{aligned} \text{kid score}_i &\sim \operatorname{Normal}(\mu_i, \sigma) && \text{Likelihood} \\ \mu_i &= \beta_0 + \beta_1 \cdot \text{mom hs}_i && \text{Lineares Modell} \\ \beta_0 &\sim \operatorname{Normal}(87, 51) && \text{Prior Achsenabschnitt} \\ \beta_1 &\sim \operatorname{Normal}(0, 124) && \text{Prior Regressionsgewicht} \\ \sigma &\sim \operatorname{Exp}(0.049) && \text{Prior Vorhersagegüte} \end{aligned} \tag{12.1}\]

Mit parameters(m10.1) bekommt man die Parameter des Modells, s. Tabelle 12.1.

Tabelle 12.1: Parameter des Modells m10.1 (sigma ist nicht dargestellt, da meistens nicht von hohem Interesse)

Parameter	Median	95% CI	pd	Rhat	ESS (tail)	Prior
(Intercept)	77.56	(73.40, 81.67)	100%	1.000	2667	Normal (86.80 +- 51.03)
mom_hs	11.75	(7.28, 16.37)	100%	1.000	2781	Normal (0.00 +- 124.21)

m10.1: kid_score = 78 + 12*mom_hs + error

Der Achsenabschnitt (intercept, \(\beta_0\) oder auch mit \(\alpha\) bezeichnet) ist der mittlere vorhergesagte IQ-Wert von Kindern, deren Mütter über keinen Schulabschluss (mom_hs = 0) verfügen:

kid_score = 78 + 0*12 + error

Das Regressionsgewicht (slope, \(\beta_1\), \(\beta\)) ist der Unterschied im IQ-Wert von Kindern mit Müttern mit Schulabschluss (im Vergleich zum IQ-Wert von Kindern mit Müttern ohne Schulabschluss). Dieser Unterschied entspricht der Steigung der Regressionsgeraden.

kid_score = 78 + 1*12 + error = 90 + error

Der Wert von error (\(e\)) zeigt, wie genau die Schätzung (Vorhersage) ist bzw. wie stark UV und AV zusammenhängen. error entspricht dem Vorhersagefehler, also dem Unterschied vom tatsächlichen IQ-Wert des Kindes (\(y\)) zum vom Modell vorhergesagten Wert (\(\hat{y}\)).

Ein lineares Modell der Art y ~ g kann man als Berechnung des Unterschieds im Mittelwert von y zwischen beiden Gruppen (g0 vs. g1) verstehen.

👨🏫 Hey Stan! Nimm den Datensatz kidiq, gruppiere nach mom_hs und fasse zusammen anhand des Mittelwerts. Die resultierende Zahl soll heißen kid_score_avg. An die Arbeit!

🤖 Schon mal was von “bitte” gehört?

R-Code

kidiq |> 
  group_by(mom_hs) |> 
  summarise(kid_score_avg = 
              mean(kid_score))

Ausgabe

mom_hs	kid_score_avg
0	77.55
1	89.32

Der mittlere (average, avg) IQ-Wert unterscheidet sich um ca. 12 Punkte (89.4-77.6), zugunsten der Kinder von Müttern mit Abschluss.

In Abbildung 12.2 ist der Unterschied im IQ der Kinder als Funktion des Schulabschlusses der Mutter dargestellt, auf Basis des Datensatzes kidiq.

estimate_relation(m10.1) |> plot()

Abbildung 12.2: Kinder, deren Mütter über einen Schulabschluss verfügen (mom_hs=1), haben im Mittel einen höheren Intelligenztestwert (als Kinder, deren Mütter keinen Schulabschluss haben, mom_hs=0), laut dem vorliegenden Modell. Die Regressionsgerade ist als durchgezogene Linie dargestellt. Die Mittelwerte pro Gruppe als Punkte und als gestrichelte, horizontale Linie.

Prüfen wir mit rope(m10.1), ob der Effekt der UV (Unterschied zwischen den Gruppen) “praktisch Null” ist; dazu nutzen wir das ROPE-Verfahren.

rope(m10.1, range = c(-5, 5))

Das Ergebnis zeigt uns, dass es 0% Überlappung vom Rope und dem 95%-HDI (der Posterior-Verteilung) gibt.

Fazit: Wir verwerfen die Praktisch-Null-Hypothese. Adios! Abbildung 12.3 visualisiert die Erstreckung der Posteriori-Verteilung (und des 95%-HDI) sowie des Rope.

(a) Diagramm mit rope(m10.1) |> plot()

(b) Diagramm mit parameters(m10.1) |> plot()

Abbildung 12.3: Der Bereich plausibler Werte (der Post-Verteilung von m10.1) liegt komplett außerhalb des Bereichs zu vernachlässigend kleiner Werte (ROPE-Bereich). Wir verwerfen die Praktisch-Null-Hypothese (ROPE-Hypothese).

In der frequentistischen Statistik (die mehrheitlich unterrichtet wird) untersucht man diese Datensituation – Mittelwertsdifferenz zwischen zwei Gruppen – mit einem sog. t-Test. Der t-Test (für zwei Gruppen) ist ein inferenzstatistisches Verfahren, das prüft, ob die Mittelwertsdifferenz zwischen zwei Gruppen (in der Population) \(\mu_d\) Null ist: \(\mu_d = 0\).⁹ In der Bayes-Statistik betrachtet man dazu stattdessen die Posteriori-Verteilung (z.B. mit 95%PI).

Alternativ zum t-Test kann man – unabhängig, ob man Frequentistisch oder Bayesianisch unterwegs ist – mit einer Regression vom Typ y ~ b das in etwa gleiche Ergebnis erreichen.¹⁰

Zur Beantwortung der Forschungsfrage betrachten wir die Ergebnisse von m10.1.

R-Code

m10.1_post <-
  m10.1 |> 
  as_tibble() 

names(m10.1_post) <- c("Achsenabschnitt", "momhs", "sigma")  # schönere Namen

Ausgabe

Hier sind die ersten paar Zeilen, s. Tabelle 12.2.

Tabelle 12.2: m10.1, Stichprobe aus der Post-Verteilung, ersten paar Zeilen

Achsenabschnitt	momhs	sigma
80.5	9.5	19.7
77.9	12.5	19.7
75.9	12.9	20.0
77.8	9.5	20.7
78.7	11.9	20.0

Berechnen wir zur Übung ein 95%-PI von Hand; komfortabler geht es mit eti(m10.1), s. Tabelle 12.1.

pi_mom_hs <-
  m10.1_post |>
  reframe(pi_95 = quantile(momhs, c(.025, .975)))

pi_mom_hs

Mit 95% Wahrscheinlichkeit liegt der Unterschied im mittleren IQ-Wert zwischen Kindern von Müttern mit bzw. ohne Schulabschluss im Bereich von 7 bis 14 IQ-Punkten, laut unserem Modell: \(95\%PI: [7,16]\). Die Hypothese, dass es keinen Unterschied oder einen Unterschied in die andere Richtung geben sollte, ist vor diesem Hintergrund als unwahrscheinlich abzulehnen. Mit plot(eti(m10.1)) kann man das 95%-ETI der Post-Verteilung visualisieren.

Wichtig

Zur Erinnerung: Korrelation ungleich Kausation. Von einem “Effekt” zu sprechen, lässt in den meisten Köpfen wohl die Assoziation zu einem kausalen Effekt entstehen. Ein Kausaleffekt ist eine starke (und sehr interessante und wichtige) Behauptung, die mehr Fundierung bedarf als eine einfache Korrelation bzw. ein einfacher Zusammenhang. Für eine Kausalaussage braucht man ein Argument, etwa einen Verweis auf bestehende Studien oder eine Theorie. \(\square\)

12.4 Vertiefung: Toleranzbereich

🏎️VERTIEFUNG, nicht prüfungsrelevant🏎️

Berechnet man ein Regressionsmodell mit stan_glm (🤖😁), dann zieht man dabei Zufallszahlen 🎲. Der Hintergrund ist, dass Stan eine Stichproben-Post-Verteilung erstellt, und das Ziehen der Stichproben erfolgt zufällig. Das erklärt, warum Ihre Ergebnisse einer Regressionsanalyse mittels stan_glm von denen in diesem Buch abweichen können.

Um zu prüfen, ob Ihre Ergebnisse “ähnlich genug” oder “innerhalb eines Toleranzbereichs” sind, kann man die Funktion is_in_tolerance() aus dem R-Paket prada nutzen.

HinweisGröße des Toleranzbereichs

Die Größe des relativen Toleranzbereichs ist in is_in_tolerance() auf 5% festgelegt. Das heißt, ein Unterschied von 5% zwischen einem Referenzwert (dem “wahren” Wert) und Ihrem Wert ist okay, also im Toleranzbereich. Außerdem gibt es noch einen absoluten Toleranzbereich, der auf 10% der SD der AV festgelegt ist (bei Regressionsmodellen). Der größere der beiden Werte gilt. \(\square\)

Wenn Sie diese Funktion nutzen wollen, müssen Sie zunächst das Paket installieren (von Github, nicht vom Standard-R-App-Store CRAN) und dann wie gewohnt starten.

pak::pak("sebastiansauer/prada")  # `pak` können Sie mit `install.packages("pak")` installieren
library(prada)

Dann testen Sie, ob Ihr Modellparameter, z.B. \(\beta_1\) innerhalb eines Toleranzbereichs liegt.

Sagen wir der “richtige” oder “wahre” Wert (oder schlicht der Wert einer Musterlösung) für \(\beta_0\) ist 77. Unser Wert sei 77.56. Liegt dieser Wert noch innerhalb eines Toleranzbereichs?

is_in_tolerance(asis = 77.56,  # Ihr Wert
                tobe = 77,   # Referenzwert
                tol_rel = .05,   # relative Toleranz
                tol_abs = .10 * sd(kidiq$kid_score)  # absolute Toleranz
                )
## [1] TRUE

Ja, unser Wert ist innerhalb des Toleranzbereichs. ✅

Übungsaufgabe 12.2 (Peer-Instruction: Wählen Sie plausible Priori-Verteilungen für Modell m10.1!) Die Apriori-Verteilungen, die Stan als Voreinstellung gewählt hat, sind recht weit. Das ist nicht unbedingt falsch, aber z.B. weniger Streuung in den Verteilungen wäre wohl plausibler.

Bearbeiten Sie dazu folgende Aufgaben in Kleingruppen:

1. Welche Parameter gibt es in m10.1?
2. Geben Sie passende Parameter für diese Apriori-Verteilungen an in mathematischer Notation.
3. Geben Sie den passenden R-Befehl für diese Apriori-Werte an.

Bereiten Sie Ihre Lösungen schriftlich elektronisch vor. Die Lehrkraft wird Sie ggf. bitten, Ihre Aufgaben im Plenum zu präsentieren. \(\square\)

12.5 y ~ x + b

Forschungsfrage:

Wie stark ist der statistische Effekt von jeweils dem Schulabschluss der Mutter (mom_hs) und dem IQ der Mutter (mom_iq) auf den IQ des Kindes (kid_score)?

Die Modellformel zur Forschungsfrage lautet: y ~ x + b bzw. kid_score ~ mom_iq + mom_hs.

Die Hypothesen lauten:

1. Der Schulabschluss der Mutter hat einen positiven Effekt auf den IQ des Kindes: \(\beta_{momhs} > 0\).
2. Der IQ der Mutter hat einen positiven Effekt auf den IQ des Kindes: \(\beta_{momiq} > 0\).

Der Kausalgraph¹¹ zur Modellformel ist in Abbildung 12.4 dargestellt. Laut unserem Modell ist y also eine Funktion zweier (kausaler) Einflüsse, b und u, wobei u für “unbekannt” steht, also für alle sonstigen Einflüsse.¹²

Abbildung 12.4: DAG für y ~ b

Deskriptive Statistiken zum Datensatz sind in Tabelle 12.3 dargestellt.

describe_distribution(kidiq)

Tabelle 12.3: Variablen und ihre Verteilung im Datensatz kidiq

Variable	Mean	SD	IQR	Min	Max	n
kid_score	86.80	20.41	28	20	144	434
mom_hs	0.79	0.41	0	0	1	434
mom_iq	100.00	15.00	22	71	139	434
mom_age	22.79	2.70	4	17	29	434

Berechnen wir als Erstes folgendes einfaches Modell mit nur einer UV: kid_score ~ mom_iq (m10.2), s. Tab. Tabelle 12.4.

m10.2 <-
  stan_glm(kid_score ~ mom_iq, data = kidiq)

m10.2 |> 
  parameters()

Tabelle 12.4: Parameter des Modells m10.2

Parameter	Median	95% CI	pd	Rhat	ESS (tail)	Prior
(Intercept)	25.78	(14.06, 37.44)	100%	1.000	2422	Normal (86.80 +- 51.03)
mom_iq	0.61	(0.49, 0.72)	100%	1.000	2391	Normal (0.00 +- 3.40)

Setzt man die (Punktschätzer der) Parameter in die allgemeine Regressionsgleichung von m10.2 ein, so erhält man: kid_score = 26 + 0.6 * mom_iq + error. Visualisieren wir uns noch das Modell m10.2, s. Abbildung 12.5.

mit ggplot2

kidiq |> 
  ggplot(aes(x = mom_iq, y = kid_score)) +
  geom_point(alpha = .7) +
  geom_abline(slope = coef(m10.2)[2],
              intercept = coef(m10.2)[1],
              color = "blue")

Abbildung 12.5: Die Intelligenz eines Kindes als Funktion der Intelligenz der Mutter (m10.2)

Mit easystats

Alternativ kann man sich das Modell (m10.2) mit plot(estimate_relation) visualisieren, mit Hilfe des R-Pakets easystats, s. Abbildung 12.6.

plot(estimate_relation(m10.2)) +
  geom_point(data = kidiq, aes(x = mom_iq, y = kid_score))

Abbildung 12.6: Die geschätzten Erwartungswerte von m10.2 visualisiert

Die Linie zeigt die vorhergesagten IQ-Werte der Kinder für verschiedene IQ-Werte der Mütter. Vergleicht man Teilpopulationen von Müttern mit einem mittleren Unterschied von einem IQ-Punkt, so findet man 0.6 IQ-Punkte Unterschied bei ihren Kindern im Durchschnitt, laut dem Modell m10.2. Der Achsenabschnitt hilft uns nicht weiter, da es keine Menschen mit einem IQ von 0 gibt.

Berechnen wir als nächstes ein Modell mit beiden Prädiktoren: kid_score ~ mom_hs + mom_iq, s. Tabelle 12.5; Modell m10.3.

m10.3 <- 
  stan_glm(
    kid_score ~ mom_iq + mom_hs, 
    refresh = 0,
    data = kidiq)

Tabelle 12.5: Parameter des Modells m10.3 (ohne sigma; ETI-Intervalle)

Parameter	Median	95% CI	pd	Rhat	ESS (tail)	Prior
(Intercept)	25.86	(14.14, 37.04)	100%	1.000	3018	Normal (86.80 +- 51.03)
mom_iq	0.56	(0.44, 0.69)	100%	1.000	2786	Normal (0.00 +- 3.40)
mom_hs	5.97	(1.72, 10.32)	99.60%	1.000	2822	Normal (0.00 +- 124.21)

Will man nur schnell die Koeffizienten des Modells (d.h. Punktschätzer der Modellparameter, in diesem Fall den Median) wissen, so kann man anstelle von parameters(mein_modell) auch coef(mein_modell) schreiben. Aber natürlich ist es möglich (und einfacher) anstelle von coef den Befehl parameters zu verwenden.

coef(m10.3)
## (Intercept)      mom_iq      mom_hs 
##       25.86        0.56        5.97

m10.3: kid_score = 26 + 0.6*mom_iq + 6*mom_hs + error

Möchte man nur z.B. den 3. Wert aus diesem Vektor, so kann man schreiben:

coef(m10.3)[3]
## mom_hs 
##      6

Abbildung 12.7 visualisiert das Modell m10.3.

kidiq2 <-
  kidiq |> 
  mutate(mom_hs = as.factor(mom_hs))

m10.3a <- 
  stan_glm(
    kid_score ~ mom_iq + mom_hs, 
    refresh = 0,
    data = kidiq2)

pred <- estimate_relation(m10.3a)

plot(pred) +
  geom_point(data = kidiq2, aes(x = mom_iq, y = kid_score, color = mom_hs)) +
  scale_color_okabeito()

Abbildung 12.7: Der Effekt von sowohl mütterlicher Intelligenz als auch mütterlichem Schulabschluss.

• Achsenabschnitt: Hat das Kind eine Mutter mit einem IQ von 0 und ohne Schulabschluss, dann schätzt das Modell den IQ-Wert des Kindes auf 26.
• Koeffizient zum mütterlichen Schulabschluss: Vergleicht man Kinder von Müttern gleicher Intelligenz, aber mit Unterschied im Schulabschluss, so sagt das Modell einen Unterschied von 6 Punkten im IQ voraus.
• Koeffizient zum mütterlichen IQ: Vergleicht man Kinder von Müttern mit gleichem Wert im Schulabschluss, aber mit 1 IQ-Punkt Unterschied, so sagt das Modell einen Unterschied von 0.6 IQ-Punkten bei den Kindern voraus.

Mit 95% Wahrscheinlichkeit liegt der Unterschied im mittleren IQ-Wert zwischen Kindern von Müttern mit bzw. ohne Schulabschluss im Bereich von 1.6 bis 10.1 IQ-Punkten, laut unserem Modell. Der Effekt des mütterlichen IQs wird auf 0.5 bis 0.7 geschätzt (95%-ETI). Da für beide UV die Null nicht im Intervall plausibler Werte liegt, kann ein Null-Effekt (die exakte Nullhypothese) abgelehnt werden.

12.6 y ~ x + b + x:b

Forschungsfrage:

Gibt es einen Interaktionseffekt zwischen mütterlichem Schulabschluss und mütterlichem IQ (auf den IQ-Wert des Kindes)?

Außerdem ist man vermutlich auch an den Effekten der beiden UV auf die AV interessiert; diese Fragen haben wir im letzten Abschnitt untersucht (und greifen sie daher nicht noch mal ausführlich auf).

Die Modellformel zur Forschungsfrage lautet: y ~ x + b + x:b. Der Einfachheit halber übernehmen wir wieder die Prioris wie vom R-Paket rstanarm bereitgestellt.

Der DAG zur Modellformel ist in Abbildung 12.8 dargestellt (momiqmom_hs bezeichnet den Interaktionseffekt der Intelligenz der Mutter mit ihrem Schlulabschluss auf die Intelligenz des Kindes).

Abbildung 12.8: DAG für y ~ x + b + x:b

In m10.3 hat das Modell die Regressionsgeraden gezwungen, parallel zu sein. Betrachtet man das Streudiagramm, so sieht man, das nicht-parallele Geraden besser passen. Sind die Regressionsgeraden nicht parallel, so spricht man von einer Interaktion (synonym: Interaktionseffekt, Moderation). Fügen wir also im nächsten Modell, m10.4 einen Interaktionseffekt hinzu.

Wichtig

Liegt eine Interaktion vor, so unterscheidet sich die Steigung der Geraden in den Gruppen. Liegt keine Interaktion vor, so sind die Geraden parallel. \(\square\)

Wir berechnen mit m10.4 das Modell mit folgender Modellformel: kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, s. Listing 12.2, Abbildung 12.9 und Tabelle 12.6.

Listing 12.2: Die Modelldefition von m10.4 mit stanglm

m10.4 <- 
  stan_glm(kid_score ~ mom_iq + mom_hs + mom_hs:mom_iq, 
           data = kidiq, 
           refresh = 0)

In der Regressionsformel sieht man, dass ein zusätzlicher Parameter, eben der Interaktionseffekt, in das Modell aufgenommen wurde.

Tabelle 12.6: Parameter von m10.4

## Parameter     |          95% CI
## -------------------------------
## (Intercept)   | [-35.42, 16.94]
## mom_iq        | [  0.67,  1.23]
## mom_hs        | [ 20.41, 78.60]
## mom_iq:mom_hs | [ -0.77, -0.15]

Mit estimate_relation(m10.4) |> plot() kann man sich das Modell visualisieren, s. Abbildung 12.9.

Abbildung 12.9: m10.4: Wie m10.3, aber mit Interaktionseffekt. Es ist gut zu erkennen, dass der Achsenabschnitt für diese Daten kaum zu interpretieren ist.

Der Achsenabschnitt gibt die IQ-Schätzwerte für Kinder mit Müttern ohne Abschluss und mit einem IQ von 0 an. Kaum zu interpretieren. mom_hs: Dieser Koeffizient zeigt den Unterschied der IQ-Schätzwerte zwischen Kindern mit Mutter ohne bzw. mit Schulabschluss und jeweils mit einem IQ von 0. Puh. mom_iq: Unterschied der IQ-Schätzwerte zwischen Kindern mit Müttern, die sich um einen IQ-Punkt unterscheiden, aber jeweils ohne Schulabschluss. Interaktion (mom_hs:mom_iq): Der Unterschied in den Steigungen der Regressionsgeraden, also der Unterschied des Koeffizienten für mom_iq zwischen Müttern mit bzw. ohne Schulabschluss. Für beide Gruppen, mom_hs=0 und mom_hs=1, gilt folgende Regressionsformel, s. Gleichung 12.2.

\[ \text{kid score} = \beta_0 + \beta_1 \cdot \text{mom hs} + \beta_2 \cdot \text{mom iq} + \beta_3 \cdot \text{mom hs} \cdot \text{mom iq} \tag{12.2}\]

\(\beta_3\) gibt die Stärke des Interaktionseffekts an.

Auf Errisch schreibt man Gleichung 12.2 so (s. Listing 12.2):

kid_score ~ mom_iq + mom_hs + mom_hs:mom_iq.

Der Doppelpunkt zwischen mom_hs und mom_iq steht für den Interaktionseffekt der beiden Variablen.

Trägt man die Punkteschätzer der Koeffizienten (\(\beta_0, \beta_1, \beta_2, \beta_3\)) ein, so erhält man Gleichung 12.3.

\[ \text{kid score} = -10 + 49.1 \cdot \text{mom hs} + 1 \cdot \text{mom iq} + -0.5 \cdot \text{mom hs} \cdot \text{mom iq} \tag{12.3}\]

Teilen wir die Regressionsformel einmal auf die beiden Gruppen (mom_hs=0 bzw. mom_hs=1) auf:

mom_hs=0:

kid_score = -10 + 49.1*0 + 1*mom_iq  - -0.5*0*mom_iq
          = -10 + 1.1*mom_iq

mom_hs=1:

kid_score = -10 + 49.1*mom_hs + 1*mom_iq  - -0.5*mom_hs*mom_iq
          = -10 + 1.1*mom_iq

Nach der Interpretation von 20 unzentrierten Koeffizienten …

via GIPHY

Wir müssen dringend die unzentrierten Prädiktoren loswerden …

Wie in Tabelle 12.6 ersichtlich, kann für alle drei Effekte (mütterliche IQ, mütterlicher Schulabschluss und Interaktion von mütterlichem IQ mit mütterlichem Schulabschluss) ein Nulleffekt ausgeschlossen werden. Ob die Effekte stärker als “praktisch Null” sind, kann mittels des ROPE-Verfahren untersucht werden.

12.7 y ~ x_c + b + x_c:b

Unter Zentrieren (to center) versteht man das Bilden der Differenz eines Messwerts zu seinem Mittelwert.¹³ Zentrierte Werte geben also an, wie weit ein Messwert vom mittleren (typischen) Messwert entfernt ist. Mit zentrierten Werten ist eine Regression einfacher zu interpretieren. Hier zentrieren wir (nur) mom_iq; die zentrierte Variable kennzeichnen wir durch den Suffix _c, also mom_iq_c. Man könnte auch mom_hs zentrieren, aber für eine einfache Interpretation ist es meist nützlich, nur metrische Prädiktoren zu zentrieren.

kidiq <-
  kidiq |> 
  mutate(mom_iq_c = mom_iq - mean(mom_iq))

m10.5 <- stan_glm(kid_score ~ mom_hs + mom_iq_c + mom_hs:mom_iq_c, 
                  data = kidiq, 
                  refresh = 0)

coef(m10.5)  # nur die Punktschätzer für die Koeffizienten ausgeben

Tabelle 12.7 zeigt die Punktschätzer der Koeffizienten von m10.5.

Tabelle 12.7: Punktschätzer von m10.5 (zentrierte UV)

Fixed Effects

Parameter	Median	CI
(Intercept)	85.32	(81.07, 89.47)
mom_hs	2.92	(-1.83, 7.52)
mom_iq_c	0.96	(0.69, 1.24)
mom_hs:mom_iq_c	-0.48	(-0.79, -0.17)

Zur Interpretation von m10.5:

• Der Achsenabschnitt (Intercept) gibt den geschätzten mittleren IQ des Kindes an, wenn man eine Mutter mittlerer Intelligenz und ohne Schulabschluss betrachtet.
• mom_hs gibt den Unterschied im geschätzten mittleren IQ des Kindes an, wenn man Mütter mittlerer Intelligenz aber mit bzw. ohne Schulabschluss vergleicht.
• mom_iq_c gibt den Unterschied im geschätzten mittleren IQ des Kindes an, wenn man Mütter ohne Schulabschluss aber mit einem IQ-Punkt Unterschied vergleicht.
• mom_hs:mom_iq_c gibt den mittleren Unterschied in den Koeffizienten für mom_iq_c an zwischen den beiden Grupen von mom_hs.

Mit estimate_relation(m10.5) |> plot() kann man sich das Modell visualisieren, s. Abbildung 12.10.

Abbildung 12.10: m10.5: Mit zentrierten Prädiktoren gibt der Achsenabschnitt den Y-Wert für eine Beobachtung mit mittleren X-Wert an; daher ist der Achsenabschnitt besser zu interpretieren als ohne Zentrierung.

Zentrieren ändert nichts an den Vorhersagen. Betrachten wir die Vorhersagen von m10.4. Mit predict kann man sich die Vorhersagen eines Modells ausgeben lassen.

new <- tibble(mom_hs = 0, mom_iq = mean(kidiq$mom_iq))
pred_new <- predict(m10.4, newdata = new)
mean(pred_new)
## [1] 85

Und vergleichen wir Vorhersagen von m10.4 mit denen von m10.5: Wir sehen, die Vorhersagen sind (bis auf Rundungsfehler) identisch.

new <- tibble(mom_hs = 0, mom_iq_c = 0)
pred_new <- predict(m10.5, newdata = new)
mean(pred_new)
## [1] 85

Auch die Streuungen der vorhergesagten Werte unterscheiden sich nicht (wirklich): \(\sigma_{m10.4}= 18\); \(\sigma_{m10.5}= 18\).

Das Zentrieren ändert auch nicht die Regressionskoeffizienten, da die Streuungen dieser Variablen nicht verändert wurden durch das Zentrieren.

Tabelle 12.8 zeigt die Punktschätzer der Parameter für m10.5 sowie ihre ETI. Nutzen Sie dafür parameters(m10.5), s. Tabelle 12.8. Highest Density (Posterior) Intervalle (HDI oder HDPI) kann man sich komfortabel ausgeben lassen mit hdi(m10.5) oder mit parameters(m10.5, ci_method = "hdi"), s. Tabelle 12.9. Im Falle symmetrischer Posteriori-Verteilungen (wie hier) kommen beide Arten von Intervallen zu gleichen Ergebnissen.

Tabelle 12.8: Parameter von m10.5 und ETIs

Fixed Effects

Parameter	Median	CI
(Intercept)	85.32	(81.07, 89.47)
mom_hs	2.92	(-1.83, 7.52)
mom_iq_c	0.96	(0.69, 1.24)
mom_hs:mom_iq_c	-0.48	(-0.79, -0.17)

parameters(m10.5, ci_method = "hdi") |> 
  select(1,2,4,5) |> 
  display()

Tabelle 12.9: Parameter von m10.5 und HDIs

Fixed Effects

Parameter	Median	CI
(Intercept)	85.32	(81.16, 89.53)
mom_hs	2.92	(-1.49, 7.79)
mom_iq_c	0.96	(0.68, 1.24)
mom_hs:mom_iq_c	-0.48	(-0.79, -0.17)

Das Model zeigt keine Belege, dass sich die mittlere Intelligenz von Kindern bei Müttern mit bzw. ohne Schulabschluss unterscheidet (95%PI: [-2.0, 7.8]); die Befundlage ist unklar. Hingegen fand sich ein Effekt der mütterlichen Intelligenz; pro Punkt Unterschied in mütterlichem IQ fand sich ein Unterschied von 0.7 bis 1.3 IQ-Punkte beim Kind (95%PI). Außerdem fand sich ein Beleg, dass der Zusammenhang des IQ zwischen Mutter und Kind durch den Schulabschluss moderiert wird: Bei Mütter mit Schulabschluss war der Regressionskoeffizient zwischen Mutter-IQ und Kind-IQ geringer (95%PI: [-0.80, -0.17]).

Wichtig

Das Modell hat mittels Abbildung 12.11 mutig Kausalaussagen postuliert. Das ist zwar schön, bedarf aber einer Begründung mit Rückgriff auf die Literatur (was hier nicht getan wurde). Ohne eine Begründung ist die Behauptung der Kausalaussage nicht zu verteidigen.

12.8 y ~ g

Hier untersuchen wir ein Modell mit einer nominalen UV mit mehreren Stufen.

Nach Ihrem Studium wurden Sie reich als Unternehmensberaterin; Ihre Kompetenz als Wirtschaftspsychologin war heiß begehrt. Von Statistik wollte niemand etwas wissen… Doch nach einiger Zeit kamen Sie in eine Sinnkrise. Sie warfen Ihren Job hin und beschlossen, in die Wissenschaft zu gehen. Kurz entschlossen bewarben Sie sich auf das erste Stellenangebot als Nachwuchswissenschaftler:in. Ihr Forschungsprojekt führte Sie in die Antarktis… Nun, das war zumindest ein Gegenentwurf zu Ihrem bisherigen Jet-Set-Leben. Ihre Aufgabe bestand nun darin, Pinguine zu untersuchen. Genauer gesagt ging es um Größenunterschiede zwischen drei Pinguinarten. Ja, stimmt, an so ein Forschungsprojekt hatten Sie vorher nie auch nur im Traum gedacht.

Forschungsfrage:

Unterscheiden sich die mittleren Körpergewichte der drei Pinguinarten?

Die allgemeine Modellformel zur Forschungsfrage lautet: y ~ g.

Der DAG zur Modellformel ist in Abbildung 12.11 dargestellt.

Abbildung 12.11: DAG für y ~ g

Prüfen wir die Nullhypothese, \(H_0\), dass alle Mittelwerte gleich sind, exakt gleich (?)

Formal: \(\mu_1 = \mu_2 = \ldots = \mu_k\) mit \(k\) verschiedenen Gruppen von Pinguinarten. Hypothesen, die keinen (Null) Unterschied zwischen Gruppen oder keinen Zusammenhang zwischen Variablen postulieren, kann man als Nullhypothesen bezeichnen. Moment. Dass sich alle Mittelwerte um 0,00000000 unterscheiden, ist wohl nicht zu vermuten. Daher ist die bessere Forschungsfrage:

Wie sehr unterscheiden sich mittlere Körpergewichte in Abhängigkeit von der Pinguinart?

Alternativ können wir die Hypothese prüfen, ob die Mittelwerte “praktisch” gleich sind, also sich “kaum” unterscheiden. Der Grenzwert für “praktisch gleich” bzw. “kaum unterschiedlich” ist subjektiv. Dazu in Kapitel 11.4 mehr.

Werfen wir einen Blick in den Datensatz penguins. Hier ist die Verteilung des Gewichts jeder Spezies im Datensatz, Tabelle 12.10.

Palmer Penguins

Datenquelle¹⁴, Beschreibung des Datensatzes¹⁵

penguins |> 
  select(body_mass_g, species) |> 
  group_by(species) |> 
  describe_distribution(range = FALSE, iqr = FALSE)

Tabelle 12.10: Die Verteilung des Körpergewichts pro Spezies der Pinguine

species	Mean	SD	IQR	Min	Max	n
Adelie	3701	459	650	2850	4775	151
Chinstrap	3733	384	488	2700	4800	68
Gentoo	5076	504	800	3950	6300	123

Was fällt Ihnen auf? Visualisieren wir als Nächstes die Verteilung der AV, um den Datensatz genauer zu verstehen. Hier kommen die Pinguine! Wie schwer sind die Tiere in unserer Stichprobe, s. Abbildung 12.12?

Abbildung 12.12: Verteilung des Körpergewichts dreier Arten von Pinguinen – Geom Violine

Berechnen wir nun das mittlere Gewicht pro Spezies (Gruppe) der Pinguine, s. m10.6 und Tabelle 12.11. Die Modellformel für m10.6 lautet also body_mass_g ~ species.

options(mc.cores = parallel::detectCores())  # Turbo einschalten

m10.6 <- stan_glm(body_mass_g ~ species, 
                  data = penguins, 
refresh = 0,  # unterdrückt Ausgabe der Posteriori-Stichproben
                  )

m10.6 |> parameters() 

Tabelle 12.11: Parameter des Modells m10.6; neben dem Achsenabschnitt sind die Effekte der Gruppe Adelie und Chinstrap ausgewiesen

Fixed Effects

Parameter	Median	CI_low	pd
(Intercept)	3702.05	3628.80	100%
speciesChinstrap	33.38	-105.46	67.60%
speciesGentoo	1375.11	1263.06	100%

Zur Interpretation von m10.6:

Die UV hat drei verschiedene Stufen (Werte, Ausprägungen; hier: Spezies), aber es werden in Tabelle 12.11 nur zwei Stufen angezeigt (also eine weniger) zusätzlich zum Achsenabschnitt. Die fehlende Stufe (Adelie, nicht ausgegeben) ist die Vergleichs- oder Referenzkategorie (baseline) und ist im Achsenabschnitt ausgedrückt (Intercept). Die Koeffizienten für species geben jeweils den (vorhergesagten) Unterschied zur Vergleichskategorie wieder. Pinguine der Spezies Adelie haben laut Modell ein mittleres Gewicht von ca. 3700g. Pinguine der Spezies Gentoo sind laut Modell im Mittel gut 1000g schwerer als Pinguine der Spezies Adelie, etc.

Der Unterschied im mittleren Gewicht von den Gruppen Chinstrap und Gentoo zur Referenzgruppe (Adelie) ist in Abbildung 12.13 verdeutlicht.

plot(hdi(m10.6)) + scale_fill_okabeito()

Abbildung 12.13: Effekt der UV: Unterschiede zur Referenzgruppe (95%-HDI)

Das Farbschema nach Okabe und Ito¹⁶ ist gut geeignet, um nominal skalierte Farben zu kodieren (s. Details hier).

Glauben wir jetzt an Gruppeneffekte?

Glauben wir jetzt, auf Basis der Modellparameter, an Unterschiede (hinsichtlich der AV) zwischen den Gruppen (UV)?

Es scheinen sich nicht alle Gruppen voneinander zu unterscheiden. So ist der Mittelwert der Gruppe Gentoo deutlich höher als der der beiden anderen Gruppen. Umgekehrt sind sich die Pinguinarten Adelie und Chinstrap in ihren Mittelwerten ziemlich ähnlich. Wie in Abbildung 12.13 ersichtlich, überlappt sich der Schätzbereich für den Parameter von Gentoo nicht mit der Null; hingegen überlappt sich der Schätzbereich des Parameters für Chinstrap deutlich mit der Nullinie. Auf Basis unseres Modells verwerfen wir also (mit hoher Sicherheit) die Hypothese, dass alle Mittelwerte exakt identisch sind. Ehrlicherweise ist sowieso zweifelhaft, dass die exakte Nullhypothese \(\mu_1 = \mu_2 = \ldots = \mu_k\) bis in die letzte Dezimale gilt. Anders gesagt: Die Wahrscheinlichkeit eines bestimmten Wertes einer stetigen Zufallsvariable ist praktisch Null. Aber: Viele Forscherinnen und Forscher prüfen gerne die Nullhypothese, daher diskutieren wir hier den Begriff der (exakten) Nullhypothese. Das Verfahren der Frequentistischen Statistik, um die Nullhypothese \(\mu_1 = \mu_2 = \ldots = \mu_k\) zu testen, nennt man Varianzanalyse (analysis of variance, kurz ANOVA). In der Bayes-Statistik nutzt man – wie immer – primär die Post-Verteilung, um Fragen der Inferenz (z.B. Gruppenunterschiede dieser Art) zu beurteilen.

Unser Modell m10.6 verwendet schwach informierte (weakly informative) Priors. Für Achsenabschnitt und die Regressionskoeffizienten trifft unser Golem Stan folgende Annahmen in der Voreinstellung:

• Achsenabschnitt und Regressionsgewichte werden als normalverteilt angenommen
• mit Mittelwert entsprechend den Stichprobendaten
• und einer Streuung des Mittelwerts, die der 2.5-fachen der Streuung in der Stichprobe entspricht
• für Sigma wird eine Exponentialverteilung mit Rate \(\lambda=1\) angenommen, skaliert mit der Streuung der AV.

Mehr Infos kann man sich mit prior_summary(modell) ausgeben lassen.

prior_summary(m10.6)
## Priors for model 'm10.6' 
## ------
## Intercept (after predictors centered)
##   Specified prior:
##     ~ normal(location = 4202, scale = 2.5)
##   Adjusted prior:
##     ~ normal(location = 4202, scale = 2005)
## 
## Coefficients
##   Specified prior:
##     ~ normal(location = [0,0], scale = [2.5,2.5])
##   Adjusted prior:
##     ~ normal(location = [0,0], scale = [5015.92,4171.63])
## 
## Auxiliary (sigma)
##   Specified prior:
##     ~ exponential(rate = 1)
##   Adjusted prior:
##     ~ exponential(rate = 0.0012)
## ------
## See help('prior_summary.stanreg') for more details

Wenn man über mehr inhaltliches Wissen verfügt, wird man die Prioris anpassen wollen. So könnte man z.B. auf Basis von Fachwissen über das Gewicht von Pinguinen postulieren, dass Adelie-Pinguine im Mittel 3000 g wiegen. Und dass die anderen zwei Pinguin-Arten im Mittel sich nicht unterscheiden vom Mittelwert der Adelie-Pinguine.

m10.6b <- stan_glm(
  body_mass_g ~ species, 
  data = penguins, 
  refresh = 0,
  prior = normal(location = c(0, 0),  # betas, Mittelwert
                 scale = c(500, 500)),  # betas, Streuung
  prior_intercept = normal(3000, 500),  # Achsenabschnitt, Mittelwert und Streuung
  prior_aux = exponential(0.001)
)
coef(m10.6b)
##      (Intercept) speciesChinstrap    speciesGentoo 
##             3707               24             1356

Anstelle von Rohwerten (hier Angabe von Gramm Gewicht) kann man die Streuung auch in z-Werten eingeben, das macht es etwas einfacher, s. m10.6c. Dazu gibt man bei dem oder den entsprechenden Parametern den Zusatz autoscale = TRUE an.

m10.6c <- stan_glm(
  body_mass_g ~ species, 
  data = penguins, 
  refresh = 0,
  prior = normal(location = c(0, 0),  # betas, Mittelwert
                 scale = c(2.5, 2.5),  # betas, Streuung
                 autoscale = TRUE),  # in z-Einheiten
  prior_intercept = normal(3000, 2.5,   # Achsenabschnitt, Mittelwert und Streuung
                           autoscale = TRUE), 
  prior_aux = exponential(1, autoscale = TRUE)
)
coef(m10.6c)
##      (Intercept) speciesChinstrap    speciesGentoo 
##             3701               32             1374

Den Parameter für die Streuung des Modells, \(\sigma\), kann man sich mit sigma(modell) ausgeben lassen.

sigma(m10.6b)
## [1] 463

Implizit bekommt man die Informationen zu \(\sigma\) mitgeteilt durch die Größe der Konfidenzintervalle.

Übrigens macht es meistens keinen Sinn, extrem weite Prioris zu definieren.¹⁷

12.9 Vertiefung: Wechsel der Referenzkategorie

species ist eine nominale Variable, da passt in R der Typ factor (Faktor) am besten. Aktuell ist der Typ noch character (Text):

penguins <- penguins |> 
  mutate(species = factor(species))

Im Standard sortiert R die Faktorstufen alphabetisch, aber man kann die Reihenfolge ändern.

levels(penguins$species)
## [1] "Adelie"    "Chinstrap" "Gentoo"

Setzen wir Gentoo als Referenzkategorie und lassen die restliche Reihenfolge, wie sie ist:

library(forcats)
penguins <- penguins |> 
  mutate(species = factor(species),
    species = fct_relevel(species, "Gentoo"))

Beachten Sie, dass dazu das Paket forcats verfügbar sein muss.

Jetzt haben wir die Referenzkategorie geändert:

levels(penguins$species)
## [1] "Gentoo"    "Adelie"    "Chinstrap"

Der Wechsel der Referenzkategorie ändert nichts Wesentliches am Modell, s. Tabelle 12.12.

m10.6a <- stan_glm(body_mass_g ~ species, data = penguins, refresh = 0)
hdi(m10.6a)

Tabelle 12.12: m10.6a mit geänderter Referenzkategorie; die Effekte der UVs bleiben gleich.

Highest Density Interval

Parameter	95% HDI
(Intercept)	[ 4995.54, 5160.03]
speciesAdelie	[-1481.34, -1266.61]
speciesChinstrap	[-1479.09, -1207.72]

12.10 y ~ x1 + x2

Hier untersuchen wir Forschungsfragen mit zwei metrischen UV (und einer metrischen AV).

Forschungsfrage:

Stehen sowohl der IQ der Mutter als auch, unabhängig davon, das Alter der Mutter im Zusammenhang mit dem IQ des Kindes?

Das ist eine deskriptive Forschungsfrage. Keine Kausalwirkung (etwa “IQ der Mutter ist die Ursache zum IQ des Kindes”) wird impliziert. Es geht in dieser Forschungsfrage rein darum, Zusammenhänge in den Daten – bzw. in der Population – aufzuzeigen. Viele Forschungsfragen gehen allerdings weiter und haben explizit Kausalwirkungen im Fokus. Für solche Fragen ist ein Kausalmodell nötig: Fachlich fundierte Annahmen über Kausalzusammenhänge zwischen UV und AV.

Was heißt, X hängt mit Y zusammen?

Der Begriff “Zusammenhang” ist nicht exakt. Häufig wird er (für metrische Variablen) verstanden als lineare Korrelation \(\rho\) bzw. \(r\) oder als lineare Regression \(\beta\) bzw. \(b\).

Der Regressionskoeffizient misst die Steigung der Regressionsgerade und zeigt, wie groß der vorhergesagte Unterschied in Y ist, wenn man zwei Personen (Beobachtungseinheiten) vergleicht, die sich um eine Einheit in X unterscheiden. Der Regressionskoeffizient wird manchmal mit dem “Effekt von X auf Y” übersetzt. Vorsicht: “Effekt” klingt nach Kausalzusammenhang. Eine Regression ist keine hinreichende Begründung für einen Kausalzusammenhang.

Der Korrelationskoeffizient misst eine Art der Stärke des linearen Zusammenhangs. Er zeigt, wie klein die Vorhersagefehler der zugehörigen Regression im Schnitt sind. Korrelation ist nicht (automatisch) Kausation.¹⁸

Es ist hilfreich, sich die Korrelationen zwischen den (metrischen) Variablen zu betrachten, bevor man ein (Regressions-)Modell aufstellt, s. Tabelle 12.13.

kidiq |> 
  correlation(bayesian = TRUE)

Tabelle 12.13: Korrelation der Variablen im Datensatz kidiq

Parameter1	Parameter2	rho	CI
kid_score	mom_hs	0.23	(0.14, 0.32)
kid_score	mom_iq	0.44	(0.37, 0.52)
kid_score	mom_age	0.09	(-4.59e-03, 0.18)
kid_score	mom_iq_c	0.44	(0.36, 0.51)
mom_hs	mom_iq	0.28	(0.19, 0.37)
mom_hs	mom_age	0.21	(0.13, 0.30)
mom_hs	mom_iq_c	0.28	(0.19, 0.37)
mom_iq	mom_age	0.09	(-2.50e-03, 0.18)
mom_iq	mom_iq_c	1.00	(1.00, 1.00)
mom_age	mom_iq_c	0.09	(-7.19e-03, 0.18)

Tabelle 12.14 zeigt die Korrelationsmatrix.

kidiq |> 
  correlation(bayesian = TRUE) |> 
  summary()

Tabelle 12.14: Die Korrelationen zwischen den Variablen der Tabelle kidiq. Die Sterne geben einen Bereich des p-Werts an (wir ignorieren die Sterne hier).

Correlation Matrix (pearson-method)

Parameter	mom_iq_c	mom_age	mom_iq	mom_hs
kid_score	0.45***	0.09	0.45***	0.24***
mom_hs	0.28***	0.21***	0.28***
mom_iq	1.00***	0.09
mom_age	0.09

p-value adjustment method: Holm (1979)

Nützlich ist auch die Visualisierung der Korrelationstabelle als Heatmap, Abbildung 12.14.

kidiq |> 
  correlation() |> 
  summary() |> 
  plot()

Abbildung 12.14: Visualisierung der Korrelationsmatrix als Heatmap

Wir berechnen zunächst jeweils eine univariate Regression, pro Prädiktor, also eine für mom_iq und eine für mom_age.

m10.7 <- stan_glm(kid_score ~ mom_iq, data = kidiq, refresh = 0)
m10.8 <- stan_glm(kid_score ~ mom_age, data = kidiq, refresh = 0)

Tabelle 12.15 zeigt die Ergebnisse für mom_iq.

Tabelle 12.15: Parameter für m10.7

Fixed Effects

Parameter	Median	CI
(Intercept)	25.68	(14.29, 37.97)
mom_iq	0.61	(0.49, 0.72)

Tabelle 12.16 zeigt die Ergebnisse für mom_age.

Tabelle 12.16: Parameter für m10.8

Fixed Effects

Parameter	Median	CI
(Intercept)	71.08	(54.69, 87.33)
mom_age	0.69	(-0.03, 1.40)

Visualisieren wir nun die univariaten Regressionen. In Abbildung 12.15 ist die univariate Regression mit jeweils einem der beiden Prädiktoren dargestellt.

m10.7: Die Steigung beträgt 0.6. m10.8: Die Steigung beträgt 0.7.

Abbildung 12.15: Zwei univariate Regressionen

Berechnen wir im Folgenden das multiple Modell (beide Prädiktoren), m10.9. m10.9 stellt das multiple Regressionsmodell dar; multipel bedeutet in diesem Fall, dass mehr als ein Prädiktor im Modell aufgenommen ist, s. Tabelle 12.17.

m10.9 <- stan_glm(kid_score ~ mom_iq + mom_age, 
                  data = kidiq, 
                  refresh = 0)
parameters(m10.9)

Tabelle 12.17: Parameter für m10.9

Fixed Effects

Parameter	Median	CI
(Intercept)	17.46	(2.06e-03, 34.54)
mom_iq	0.60	(0.49, 0.72)
mom_age	0.40	(-0.26, 1.05)

Wichtig

Die Regressionsgewichte unterscheiden sich (potenziell) von denen der jeweiligen univariaten Regressionen.

Bei einer multiplen Regression ist ein Regressionsgewicht jeweils “bereinigt” vom Zusammenhang mit dem (oder den) anderen Regressionsgewicht(en). Das bedeutet anschaulich, man betrachtet den Zusammenhang einer UV mit der AV, wobei man gleichzeitig den anderen Prädiktor konstant hält.

In Abbildung 12.16 ist das Modell m10.9 in 3D dargestellt via Plotly¹⁹.

Abbildung 12.16: 3D-Visualisierung von m10.9 (zwei Prädiktoren)

Abbildung 12.17 zeigt eine Visualisierung von m10.9, in der die 3. Dimension durch eine Farbschattierung ersetzt ist.

Abbildung 12.17: Modell m10.9; die Farbverläufe zeigen der Wert der abhängigen Variablen

Auf der Achse von mom_iq erkennt man deutlich (anhand der Farbänderung) die Veränderung für die AV (kid_score). Auf der Achse für mom_age sieht man, dass sich die AV kaum ändert, wenn sich mom_age ändert.

Abbildung 12.18 visualisiert den Zusammenhang von 10 Variablen untereinander.

Abbildung 12.18: So sieht der Zusammenhang im 10-dimensionalen Raum aus

Leider macht mein Hirn hier nicht mit. Unsere Schwäche, eine große Zahl an Dimensionen zu visualisieren, ist der Grund, warum wir mathematische Modelle brauchen. Daher kann man ein Modell verstehen als eine Zusammenfassung eines (ggf. hochdimensionalen) Variablenraums.

12.11 y ~ x1_z + x2_z

In diesem Abschnitt untersuchen wir ein Modell mit zwei z-standardisierten, metrischen Prädiktoren (und einer metrischen, nicht-standardisierten AV).

Zur Relevanz der Prädiktoren: Woher weiß man, welche UV am stärksten mit der AV zusammenhängt? Man könnte auch sagen: Welcher Prädiktor (welche UV) am “wichtigsten” ist oder den “stärksten Einfluss” auf die AV ausübt? Bei solchen kausal konnotierten Ausdrücken muss man vorsichtig sein: Die Regressionsanalyse als solche ist keine Kausalanalyse. Die Regressionsanalyse – wie jede statistische Methode – kann für sich nur Muster in den Daten, also Zusammenhänge bzw. Unterschiede, entdecken. Möchte man die Relevanz von Prädiktoren vergleichen, so ist ein Kausalmodell empfehlenswert.

Welcher Prädiktor ist nun “wichtiger” oder “stärker” in Bezug auf den Zusammenhang mit der AV, mom_iq oder mom_age (Modell m10.9)? Die Antwort hängt auch von der Streuung bzw. Skalierung der Variablen ab. mom_iq hat den größeren Koeffizienten und viel Streuung; mom_age hat weniger Streuung.

Um die Relevanz der Prädiktoren vergleichen zu können, müsste man vielleicht die Veränderung von kid_score betrachten, wenn man vom kleinsten zum größten Prädiktorwert geht. Allerdings sind Extremwerte meist instabil (da sie von einer einzigen Beobachtung bestimmt werden). Sinnvoller ist es daher, die Veränderung in der AV zu betrachten, wenn man den Prädiktor von “unterdurchschnittlich” auf “überdurchschnittlich” ändert. Das kann man mit z-Standardisierung erreichen.

kidiq2 <- 
  kidiq |> 
  mutate(mom_iq_z = ((mom_iq - mean(mom_iq)) / sd(mom_iq)))  |>  # z-Transformation
  select(mom_iq, mom_iq_z) 

kidiq2 |> 
  head()

Der Nutzen von Standardisieren (dieser Art) ist die bessere Vergleichbarkeit der Effekte von UV, die (zuvor) verschiedene Mittelwerte und Streuungen hatten²⁰. Die Standardisierung ist ähnlich zur Vergabe von Prozenträngen: “Dieser Messwert gehört zu den Top-3-Prozent”. Diese Aussage ist bedeutsam für Variablen mit verschiedenem Mittelwert und Streuung. So werden vergleichende Aussagen für verschiedene Verteilungen möglich.

Zu den Statistiken zu den z-transformierten Variablen: Tabelle 12.3 zeigt die Verteilung der (metrischen) Variablen im Datensatz kidiq.

Metrische Variablen in z-Werte zu transformieren, hat verschiedene Vorteile:

• der Achsenabschnitt ist einfacher zu interpretieren (da er sich dann auf ein Objekt mit mittlerer Ausprägung bezieht)
• Interaktionen sind einfacher zu interpretieren (aus dem gleichen Grund)
• Prioriwerte sind einfacher zu definieren (wieder aus dem gleichen Grund)
• die Effekte verschiedener Prädiktoren sind einfacher in ihrer Größe zu vergleichen, da sie dann die gleiche Skalierung/Streuung aufweisen
• kleine und ähnlich große Wertebereiche erleichtern dem Golem die Rechenarbeit

Man kann die z-Transformation (“Standardisierung”) mit standardize (aus easystats) durchführen, s. Tabelle 12.18.

kidiq_z <- 
  standardize(kidiq, append = TRUE)  # z-transformiert alle numerischen Werte

Tabelle 12.18: z-transformierte Variablen im Datensatz kidiq (erste paar Zeilen und Spalten)

kid_score	mom_hs	mom_iq	mom_age
65	1	121.12	27
98	1	89.36	25
85	1	115.44	27
83	1	99.45	25
115	1	92.75	27
98	0	107.90	18

Der Schalter append = TRUE sorgt dafür, dass die ursprünglichen Variablen beim z-Standardisieren nicht überschrieben werden, sondern angehängt werden (mit einem Suffix _z).

Man kann auch nur einzelne Variablen mit standardize standardisieren, indem man das Argument select nutzt.

kidiq |> 
  standardize(select = c("mom_iq", "mom_age", "kid_score"))

Man kann das Standardisieren auch von Hand machen, ohne ein Extra-Paket, s. Tabelle 12.19. Dazu verwendet man den Befehl scale().

kidiq |> 
  mutate(mom_iq_z2 = scale(mom_iq),
         mom_age_z2 = scale(mom_age),
         kid_score_z2 = scale(kid_score))

Tabelle 12.19: Z-Standardisierung ohne Extrapaket; ersten paar Spalten und Zeilen

Berechnen wir das Modell m10.10: y ~ x1_z + x2_z.

m10.10 <- stan_glm(kid_score ~ mom_iq_z + mom_age_z, 
                   data = kidiq_z, 
                   refresh = 0)
parameters(m10.10)

Tabelle 12.20: Parameter von m10.12

Fixed Effects

Parameter	Median	CI
(Intercept)	86.80	(85.11, 88.50)
mom_iq_z	9.06	(7.41, 10.75)
mom_age_z	1.05	(-0.68, 2.79)

12.12 y_z ~ x1_z + x2_z

In diesem Abschnitt berechnen wir ein Modell (Modell m10.12), in dem sowohl die Prädiktoren z-transformiert sind (standardisiert) als auch die AV. Das z-Standardisieren der AV, kid_score, ist zwar nicht nötig, um den Effekt der Prädiktoren (UV) auf die AV zu untersuchen. Standardisiert man aber die AV, so liefern die Regressionskoeffizienten (Betas) eine Aussage darüber, um wie viele SD-Einheiten sich die AV verändert, wenn sich ein Prädiktor um eine SD-Einheit verändert. Das kann auch eine interessante(re) Aussage sein.

Berechnen wir das Modell m10.12: y_z ~ x1_z + x2_z.

m10.12 <- stan_glm(kid_score_z ~ mom_iq_z + mom_age_z, 
                   data = kidiq_z, 
                   refresh = 0)
coef(m10.12)
## (Intercept)    mom_iq_z   mom_age_z 
##     0.00056     0.44320     0.05211

• Der Achsenabschnitt gibt den Mittelwert der AV (kid_score) an, da kid_score_z = 0 identisch ist zum Mittelwert von kid_score.
• Der Koeffizient für mom_iq_z gibt an, um wie viele SD-Einheiten sich kid_score (die AV) ändert, wenn sich mom_iq um eine SD-Einheit ändert.
• Der Koeffizient für mom_age_z gibt an, um wie viele SD-Einheiten sich kid_score (die AV) ändert, wenn sich mom_age um eine SD-Einheit ändert.

Jetzt sind die Prädiktoren in ihrer Relevanz (Zusammenhang mit der AV) vergleichbar. Man sieht, dass die Intelligenz der Mutter deutlich wichtiger ist als das Alter der Mutter (im Hinblick auf die Vorhersage bzw. den Zusammenhang mit der AV).

Mit parameters können wir uns ein PI für m10.12 ausgeben lassen, s. Abbildung 12.19; im Standard wird ein 95%-ETI berichtet²¹.

parameters(m10.12) 

Fixed Effects

Parameter	Median	CI
(Intercept)	5.64e-04	(-0.08, 0.08)
mom_iq_z	0.44	(0.36, 0.52)
mom_age_z	0.05	(-0.04, 0.14)

plot(eti(m10.12)) + scale_fill_okabeito()

Abbildung 12.19: Im Standard wird ein 95%-Intervall gezeigt bzw. berechnet; hier das ETI für m10.12

Wie ist es um die Modellgüte in m10.12 bestellt? Berechnen wir \(R^2\).

r2(m10.12)
## # Bayesian R2 with Compatibility Interval
## 
##   Conditional R2: 0.204 (95% CI [0.145, 0.266])

Ist dieser Wert von \(R^2\) “gut”? Diese Frage ist ähnlich zur Frage “Ist das viel Geld?”; man kann die Frage nur im Kontext beantworten.

Eine einfache Lösung ist immer, Modelle zu vergleichen. Dann kann man angeben, welches Modell die Daten am besten erklärt, z.B. auf Basis von \(R^2\).

Zu beachten ist, dass das Modell theoretisch fundiert sein sollte. Vergleicht man viele Modelle aufs Geratewohl, so muss man von zufällig hohen Werten der Modellgüte im Einzelfall ausgehen.

Wenn Sie aber unbedingt eine “objektive” Antwort auf die Frage “wie viel ist viel?” haben wollen, ziehen wir Herrn Cohen zu Rate, der eine Antwort auf die Frage “Wie viel ist viel?” gegeben hat (Cohen, 1992):

interpret_r2(0.2)  # aus `easystats`
## [1] "moderate"
## (Rules: cohen1988)

Danke, Herr Cohen!

Die Prioris für m10.12 kann man sich mit prior_summary(m10.12) ausgeben lassen. Danke, Stan!

prior_summary(m10.12)  # aus rstanarm
## Priors for model 'm10.12' 
## ------
## Intercept (after predictors centered)
##  ~ normal(location = -2.8e-16, scale = 2.5)
## 
## Coefficients
##  ~ normal(location = [0,0], scale = [2.5,2.5])
## 
## Auxiliary (sigma)
##  ~ exponential(rate = 1)
## ------
## See help('prior_summary.stanreg') for more details

🤖 Nix zu danken!

Wie gesagt, Stan nimmt dafür einfach die empirischen Mittelwerte und Streuungen her.²²

Stans Ausgabe kann man in Mathe-Sprech so darstellen, s. Gleichung 12.4.

\[ \begin{aligned} \text{kidscore}^z_i &\sim \mathcal{N}(\mu_i,\sigma)\\ \mu_i &= \beta_0 + \beta_1\text{momiq}_i^z + \beta_2\text{momage}_i^z \\ \beta_0 &\sim \mathcal{N}(0,2.5)\\ \beta_1 &\sim \mathcal{N}(0,2.5)\\ \beta_2 &\sim \mathcal{N}(0,2.5)\\ \sigma &\sim \mathcal{E}(1) \end{aligned} \tag{12.4}\]

Man beachte, dass der Achsenabschnitt zur Intelligenz der Kinder auf Null festgelegt wird: Bei mittlerer Intelligenz und mittlerem Alter der Mutter wird mittlere Intelligenz des Kindes erwartet in m10.12. Dadurch, dass nicht nur UV, sondern auch AV z-standardisiert (d.h. zentriert und in der Streuung auf 1 standardisiert) sind, ist der Mittelwert der AV Null. Schreibt man einen Bericht, so bietet es sich an, die Modelldefinition zumindest im Anhang aufzuführen.

Übungsaufgabe 12.3 (Anzahl der Modellparameter) Wie viele Modellparameter hat m10.12?²³

Unsere Antwort auf die Forschungsfrage:

Das Modell spricht sich klar für einen statistischen, linearen Effekt von Intelligenz der Mutter auf die Intelligenz des Kindes aus, wenn das Alter der Mutter statistisch kontrolliert wird (95%PI: [0.38, 0.51]). Pro Einheit Standardabweichung in der UV (Intelligenz der Mutter) ändert sich die AV um ca. 0.44 Standardabweichungseinheiten. Hingegen zeigt das Modell, dass das Alter der Mutter statistisch eher keine Rolle spielt (95%PI: [-0.02, 0.12]). Alle Variablen wurden z-transformiert. Insgesamt erklärt das Modell im Median einen Anteil von ca. 0.2 an der Varianz der Kinderintelligenz. Das Modell griff auf die Standard-Priori-Werte aus dem R-Paket rstanarm (Goodrich et al., 2020) zurück (s. Anhang für Details).

Wichtig

Hier wird von einem “statistischen Effekt” gesprochen, um klar zu machen, dass es sich lediglich um assoziative Zusammenhänge, und nicht um kausale Zusammenhänge, handelt. Kausale Zusammenhänge dürfen wir nur verkünden, wenn wir sie a) explizit untersuchen und b) sich in der Literatur Belege dafür finden oder c) wir ein Experiment fachgerecht durchgeführt haben.

12.13 Vertiefung

Sind X und Y z-standardisiert, so sind Korrelation und Regression identisch, s. Theorem 12.2.

Theorem 12.2 Man kann die Regression als Korrelation verstehen:

\[b = r \frac{sd_y}{sd_x}\quad \square\]

Berechnen wir dazu ein einfaches Modell mit z-standardisierten Variablen und betrachten die Punktschätzer für die Regressionskoeffizienten, s. m10.12.

m10.12 <- 
  stan_glm(kid_score_z ~ mom_iq_z , data = kidiq_z, refresh = 0)
coef(m10.12)
## (Intercept)    mom_iq_z 
##       0.001       0.446

Vergleichen Sie diese Werte mit der Korrelation, s. Tabelle 12.21.²⁴

kidiq_z |> 
  select(kid_score, mom_iq, kid_score_z, mom_iq_z) |> 
  correlation() |> 
  display()

Tabelle 12.21: Correlation Matrix (pearson-method)

Parameter1	Parameter2	r	95% CI	t(432)	p
kid_score	mom_iq	0.45	(0.37, 0.52)	10.42	< .001***
kid_score	kid_score_z	1.00	(1.00, 1.00)	Inf	< .001***
kid_score	mom_iq_z	0.45	(0.37, 0.52)	10.42	< .001***
mom_iq	kid_score_z	0.45	(0.37, 0.52)	10.42	< .001***
mom_iq	mom_iq_z	1.00	(1.00, 1.00)	Inf	< .001***
kid_score_z	mom_iq_z	0.45	(0.37, 0.52)	10.42	< .001***

p-value adjustment method: Holm (1979) Observations: 434

Korrelationen der z-transformierten Variablen im Datensatz kidiq

Zentrale Annahme eines linearen Modells: Die AV ist eine lineare Funktion der einzelnen Prädiktoren, \(y= \beta_0 + \beta_1x_1 + \beta_2 x_2 + \cdots\), vgl. Theorem 2.1.

Hingegen ist es weniger wichtig, dass die AV (y) normalverteilt ist. Zwar nimmt die Regression häufig normalverteilte Residuen an²⁵, aber diese Annahme ist nicht wichtig, wenn es nur darum geht, die Regressionskoeffizienten zu schätzen (Gelman et al., 2021).

Ist die Linearitätsannahme erfüllt, so sollte der Residualplot nur zufällige Streuung um \(y=0\) herum zeigen, s. Abbildung 12.20.

Ein Residuum \(e\) ist der Vorhersagefehler, also die Differenz zwischen vorhergesagtem und tatsächlichem Wert: \(e_i = y_i - \hat{y}_i\)

kidiq <-
  kidiq |> 
  mutate(m10.12_pred = predict(m10.12),  # vorhergesagten Werte
         m10.12_resid = resid(m10.12))  # Residuen

kidiq |> 
  ggplot(aes(x = m10.12_pred, y = m10.12_resid)) +
  geom_hline(color="white", yintercept = 0, size = 2) +
  geom_hline(color = "grey40", 
             yintercept = c(-1,1), 
             size = 1, 
             linetype = "dashed") +
  geom_point(alpha = .7) +
  geom_smooth()

Abbildung 12.20: Die Verteilung der Fehler scheint keinem starken Trend (in Abhängigkeit zum vorhergesagten Wert) zu folgen, was ein gutes Zeichen ist.

Hier erkennt man keine größeren Auffälligkeiten.

Mit der PPV kann man die Güte des Modells prüfen.

pp_check(m10.12)

Unser Modell – bzw. die Stichproben unserer Posteriori-Verteilung, \(y_{rep}\) – verfehlt den Mittelwert von \(y\) leider recht häufig.

Visualisieren wir noch die bereinigten Regressionskoeffizienten, s. Abbildung 12.21.

Abbildung 12.21: Bereinigte Regressionskoeffizienten

Abbildung 12.21 zeigt in der oberen Reihe die Regression eines Prädiktors auf den anderen Prädiktor. Untere Reihe: Regression der Residuen der oberen Reihe auf die AV, kid_score_z. Unten links (C): Die Residuen von mom_iq_c sind kaum mit der AV assoziiert. Das heißt, nutzt man den Teil von mom_age_z, der nicht mit mom_iq_z zusammenhängt, um kid_score vorherzusagen, findet man keinen (kaum) Zusammenhang. Unten rechts (D): Die Residuen von mom_age_c sind stark mit der AV assoziiert. Das heißt, nutzt man den Teil von mom_iq_z, der nicht mit mom_age_z zusammenhängt, um kid_score vorherzusagen, findet man einen starken Zusammenhang.

Eine multiple Regression liefert die gleichen Regressionskoeffizienten wie die Modelle aus Teildiagrammen (C) und (D).

Was ist eigentlich der numerische Unterschied zwischen Bayesianischen und Frequentistischen Modellen? Schauen wir es uns an.

Wie man sieht liefern stan_glm() und lm ähnliche Ergebnisse (bei schwach informativen Prioriwerten):

stan_glm(mpg ~ hp + cyl, data = mtcars, refresh = 0) |> coef()
## (Intercept)          hp         cyl 
##      36.901      -0.019      -2.255

lm(mpg ~ hp + cyl, data = mtcars) |> coef()
## (Intercept)          hp         cyl 
##      36.908      -0.019      -2.265

Wichtig

Wenn auch die Ergebnisse eines frequentistischen und eines Bayes-Modells numerisch ähnlich sein können, so ist doch die Interpretation grundverschieden. Bayesmodelle erlauben Wahrscheinlichkeitsaussagen zu den Parametern, Frequentistische Modelle nicht.

12.14 Fazit

Übungsaufgabe 12.4 (Ausstieg: Was haben Sie heute (nicht) verstanden?) Schreiben Sie zum Abschluss zwei Zettel:

• Zettel 1: Schreiben Sie eine Sache, die Sie gut verstanden haben.
• Zettel 2: Schreiben Sie eine Sache, die Sie NICHT gut verstanden haben.

Legen Sie beim Verlassen des Raumes Ihre Zettel vorne bei der Lehrkraft ab (anonym). \(\square\)

Eine Kurzdarstellung des Bayes-Inferenz findet sich in diesem Post und in diesem.

📺 Musterlösung und Aufgabe im Detail besprochen – Bayes-Modell: mtcars

📺 Musterlösung und Aufgabe im Detail besprochen – Bayes-Modell: CovidIstress

Ausblick: Binäre AV. Man kann eine Regression auch verwenden, wenn die AV binär ist.

Forschungsfrage: Kann man anhand des Spritverbrauchs vorhersagen, ob ein Auto eine Automatik- bzw. ein manuelle Schaltung hat? Anders gesagt: Hängen Spritverbrauch und Getriebeart? (Datensatz mtcars)

Dazu nutzen wir den Datensatz mtcars, wobei wir die Variablen z-standardisieren.

data(mtcars)
mtcars2 <-
  mtcars |> 
  standardize(append = TRUE)

Dann berechnen wir mit Hilfe von Stan ein Regressionsmodell: m14: am ~ mpg_z:

m14 <-
  stan_glm(am ~ mpg_z, 
           data = mtcars2, 
           refresh = 0)
coef(m14)
## (Intercept)       mpg_z 
##        0.41        0.30

Ab mpg_z = 0.41, 0.3 sagt das Modell am=1 (manuell) vorher. Ganz ok.

mtcars2 |> 
  ggplot(aes(x = mpg_z, y = am)) +
  geom_hline(yintercept = 0.5, color = "white", size = 2) +
  geom_point() +
  geom_abline(intercept = coef(m14)[1],
              slope = coef(m14)[2],
              color = "blue") 

neg_am <- predict(m14, newdata = tibble(mpg_z = -1.3))

Für kleine Werte von mpg_z (<1.3) sagt unser Modell negative Werte für am voraus. Das macht keinen Sinn: Es gibt keine negativen Werte von am, nur 0 und 1. Müssen wir mal bei Gelegenheit besser machen.

Genug für heute. Wir waren fleißig …

Quelle²⁶

Wichtig

Kontinuierliches Lernen ist der Schlüssel zum Erfolg.

Weitere Hinweise zu den Themen dieses Kapitels finden sich bei Gelman et al. (2021), Kap. 10, insbesondere 10.3.

Gelman et al. (2021) bieten einen Zugang mittleren Anspruchs zur Regressionsmodellierung. Das Buch ist von einem weltweit führenden Statistiker geschrieben und vermittelt tiefe Einblicke bei gleichzeitig überschaubarem mathematischen Aufwand.

Für das vorliegende Kapitel sind insbesondere daraus die Kapitel 6, 7, und 10 relevant.

12.15 Aufgaben

12.15.1 Papier-und-Bleistift-Aufgaben

1. anz-params
2. fofrage-regrformel2
3. modelldef-regrformel
4. finde-prior
5. Nullhyp-Beispiel
6. Griech-Buchstaben-Inferenz
7. Interaktionseffekt1
8. Regression2
9. Regression3
10. diamonds-nullhyp-mws
11. zwert-berechnen
12. stan_glm_parameterzahl
13. kausale-verben
14. zuordnung_regr_formel

12.15.2 Aufgaben, für die man einen Computer benötigt

1. Regr-Bayes-interpret
2. Regr-Bayes-interpret03
3. Regr-Bayes-interpret02
4. rope4

12.15.3 Vertiefende Aufgaben

1. Anova-skalenniveau
2. ttest-skalenniveau
3. stan_glm_prioriwerte

12.15.4 Quiz-Aufgaben

Hier finden Sie Single-Choice-Aufgaben zu diesem Kapitel. Wählen Sie eine Antwort aus und klicken Sie auf das Häkchen, um sie zu überprüfen; über das Fragezeichen erhalten Sie die ausführliche Lösung.

Für das Modell kid_score ~ mom_hs (m10.1) ergibt sich näherungsweise: \(\text{kid\_score} = 78 + 12\cdot\text{mom\_hs} + \text{error}\). Welche Aussage zur Interpretation von Achsenabschnitt und Regressionsgewicht ist korrekt?

Der Achsenabschnitt (78) ist der vorhergesagte mittlere IQ von Kindern, deren Mütter keinen Schulabschluss haben (\(\text{mom\_hs}=0\)); das Regressionsgewicht (12) ist der Unterschied im vorhergesagten IQ zwischen Kindern von Müttern mit und ohne Abschluss.Der Fehlerterm (error) ist bei jedem Kind exakt gleich groß.Bei \(\text{mom\_hs}=1\) sagt das Modell einen mittleren IQ von 78 voraus.Das Regressionsgewicht (12) ist der vorhergesagte mittlere IQ von Kindern, deren Mütter einen Schulabschluss haben.Der Achsenabschnitt gibt den mittleren IQ aller Kinder an, unabhängig vom Schulabschluss der Mutter.

Bei einem Modell der Form y ~ b mit binärem Prädiktor gibt der Achsenabschnitt den vorhergesagten Wert für die Referenzgruppe (\(\text{mom\_hs}=0\)) an: \(78 + 12\cdot 0 = 78\). Das Regressionsgewicht (12) ist die Differenz zur anderen Gruppe: Für \(\text{mom\_hs}=1\) ergibt sich \(78+12\cdot 1=90\), also 12 Punkte mehr als in der Referenzgruppe. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option ignoriert, dass der Achsenabschnitt sich spezifisch auf die Referenzgruppe (\(\text{mom\_hs}=0\)) bezieht, nicht auf den (unbedingten) Gesamtmittelwert aller Kinder.

• Die dritte Option verwechselt das Regressionsgewicht (12, eine Differenz) mit dem vorhergesagten absoluten Wert für die Gruppe \(\text{mom\_hs}=1\) (der tatsächlich \(78+12=90\) beträgt).

• Die vierte Option widerspricht der Definition des Fehlerterms: error bildet gerade die individuelle, von Kind zu Kind unterschiedliche Abweichung vom vorhergesagten Wert ab, ist also nicht konstant.

• Die fünfte Option verwechselt die beiden Gruppen: Ein vorhergesagter Wert von 78 (ohne Zuschlag) gilt für \(\text{mom\_hs}=0\), nicht für \(\text{mom\_hs}=1\).

• Wahr

• Falsch

• Falsch

• Falsch

• Falsch

Das Kapitel betont wiederholt, lieber von einem “statistischen Effekt” als von einem “kausalen Effekt” zu sprechen, wenn man Regressionskoeffizienten (z.B. den Effekt des mütterlichen Schulabschlusses auf den Kindes-IQ) interpretiert. Warum?

Ein “statistischer Effekt” ist per Definition immer größer als ein “kausaler Effekt”.Eine Regressionsanalyse deckt für sich genommen nur Muster (Zusammenhänge) in den Daten auf; eine kausale Interpretation erfordert eine zusätzliche Begründung, etwa durch Rückgriff auf Literatur, Theorie oder ein sauber durchgeführtes Experiment.“Statistischer Effekt” und “kausaler Effekt” bedeuten exakt dasselbe – die Unterscheidung ist reine Stilfrage ohne inhaltliche Konsequenz.Nur experimentelle Studien dürfen überhaupt Regressionsmodelle zur Auswertung verwenden.Kausale Effekte lassen sich in der Bayes-Statistik grundsätzlich nicht berechnen, im Gegensatz zur frequentistischen Statistik.

Das Kapitel weist explizit darauf hin: Ein Kausaleffekt ist eine starke Behauptung, die mehr Fundierung benötigt als ein bloßer statistischer Zusammenhang – man braucht ein zusätzliches Argument (Literatur, Theorie oder Experiment), um von Korrelation auf Kausalität zu schließen. Die Regressionsanalyse selbst liefert nur Muster in den Daten. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option widerspricht der zentralen Warnung des Kapitels: Die begriffliche Unterscheidung ist inhaltlich wichtig, weil “kausal” eine deutlich stärkere (und riskantere) Behauptung ist als “statistisch assoziiert”.

• Die dritte Option ist falsch: Kausale Effekte lassen sich grundsätzlich sowohl in der Bayes- als auch in der frequentistischen Statistik untersuchen (z.B. mit geeigneten Kausalmodellen/DAGs) – das Problem liegt nicht am statistischen Paradigma, sondern an der Notwendigkeit zusätzlicher (nicht rein statistischer) Annahmen.

• Die vierte Option stellt eine unsinnige, im Kapitel nicht getroffene quantitative Behauptung auf – die Größe eines Effekts hängt nicht davon ab, ob man ihn “statistisch” oder “kausal” nennt.

• Die fünfte Option ist falsch: Regressionsmodelle werden im Kapitel explizit auch für Beobachtungsdaten (wie den Kung- oder Kidiq-Datensatz) verwendet, ohne dass diese aus Experimenten stammen – nur die kausale Interpretation der Ergebnisse erfordert zusätzliche Vorsicht.

• Falsch

• Wahr

• Falsch

• Falsch

• Falsch

Im Modell kid_score ~ mom_iq + mom_hs (m10.3) wird der Koeffizient für mom_hs auf ca. 6 geschätzt (deutlich weniger als die ca. 12 im einfachen Modell kid_score ~ mom_hs). Wie ist der Wert 6 in diesem Mehrfach-Regressionsmodell korrekt zu interpretieren?

Der Wert 6 gibt die Korrelation zwischen mom_hs und kid_score an.Der Gesamtunterschied im IQ zwischen allen Kindern mit und ohne Mutter-Abschluss (unabhängig vom Mutter-IQ) beträgt laut diesem Modell 6 Punkte.Der Wert 6 ist der Anteil der durch mom_hs erklärten Varianz an kid_score.Vergleicht man Kinder von Müttern mit gleichem IQ, aber unterschiedlichem Schulabschluss, sagt das Modell im Schnitt einen Unterschied von 6 IQ-Punkten voraus.Der Wert 6 gibt an, um wie viele Punkte sich der IQ der Mütter im Schnitt zwischen den beiden Gruppen unterscheidet.

In einem Mehrfach-Regressionsmodell gibt der Koeffizient eines Prädiktors den Effekt an, wenn die übrigen Prädiktoren im Modell konstant gehalten werden (ceteris paribus) – hier also: den Unterschied im vorhergesagten IQ zwischen Kindern von Müttern mit vs. ohne Schulabschluss, bei gleichem mütterlichem IQ. Das ist auch der Grund, warum dieser Koeffizient (6) kleiner ausfällt als im einfachen Modell (12): Ein Teil des im einfachen Modell gemessenen Unterschieds wird durch den mit mom_hs zusammenhängenden mütterlichen IQ miterklärt. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option beschreibt eher das Ergebnis des einfachen Modells kid_score ~ mom_hs (unbedingter Gruppenunterschied, ca. 12 Punkte) – im kontrollierten Modell m10.3 ist der Effekt aber gerade bedingt auf gleichen mütterlichen IQ, daher der kleinere Wert.

• Die dritte Option verwechselt den zu erklärenden Zusammenhang (kid_score als AV): Der Koeffizient von mom_hs sagt etwas über Kindes-IQ aus, nicht über Unterschiede im Mütter-IQ zwischen den beiden Gruppen.

• Die vierte Option verwechselt den Regressionskoeffizienten mit einem Anteil erklärter Varianz (wie \(R^2\)) – das sind unterschiedliche Kennzahlen mit unterschiedlichen Einheiten und Wertebereichen.

• Die fünfte Option verwechselt den (unstandardisierten) Regressionskoeffizienten mit einem Korrelationskoeffizienten – Letzterer ist auf den Bereich \([-1,1]\) beschränkt, ein Regressionskoeffizient wie 6 (in IQ-Punkten) hingegen nicht.

• Falsch

• Falsch

• Falsch

• Wahr

• Falsch

Welche Aussage zum Konzept der Interaktion (Moderation) in einem Regressionsmodell wie kid_score ~ mom_iq + mom_hs + mom_iq:mom_hs ist korrekt?

Liegt eine Interaktion vor, müssen die Haupteffekte mom_iq und mom_hs zwingend aus dem Modell entfernt werden.Eine Interaktion kann nur zwischen zwei nominal skalierten Variablen bestehen, nicht zwischen einer metrischen und einer nominalen Variable.Eine Interaktion liegt vor, wenn sich die Steigung der Regressionsgeraden zwischen den Gruppen (z.B. \(\text{mom\_hs}=0\) vs. \(\text{mom\_hs}=1\)) unterscheidet; sind die Geraden parallel, liegt keine Interaktion vor.Der Interaktionsterm gibt den quadrierten Effekt von mom_iq an.Eine Interaktion bedeutet, dass die beiden Prädiktoren perfekt miteinander korreliert sind.

Wie im Kapitel festgehalten: “Liegt eine Interaktion vor, so unterscheidet sich die Steigung der Geraden in den Gruppen. Liegt keine Interaktion vor, so sind die Geraden parallel.” Die Interaktion misst also, ob (und wie stark) sich der Zusammenhang zwischen einem Prädiktor und der AV je nach Ausprägung eines anderen Prädiktors unterscheidet. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option verwechselt Interaktion mit Multikollinearität: Eine Interaktion sagt nichts über die Korrelation zwischen den beiden Prädiktoren aus, sondern darüber, wie der Effekt des einen Prädiktors vom Wert des anderen abhängt.

• Die dritte Option beschreibt einen quadratischen Term (\(\text{mom\_iq}^2\)), nicht einen Interaktionsterm (Produkt zweier verschiedener Prädiktoren, hier \(\text{mom\_iq}\cdot\text{mom\_hs}\)).

• Die vierte Option ist eine gängige Fehlvorstellung: Auch bei vorhandener Interaktion bleiben die Haupteffekte im Modell (und behalten dann eine spezifische, meist bedingte Interpretation, z.B. “Effekt bei Referenzwert des anderen Prädiktors”), sie werden nicht entfernt.

• Die fünfte Option ist falsch: Das Beispiel im Kapitel (mom_iq:mom_hs) zeigt gerade eine Interaktion zwischen einer metrischen Variable (mom_iq) und einer nominalen/binären Variable (mom_hs) – Interaktionen sind für beliebige Kombinationen von Skalenniveaus definierbar.

• Falsch

• Falsch

• Wahr

• Falsch

• Falsch

Im Interaktionsmodell wird empfohlen, den metrischen Prädiktor mom_iq zu zentrieren (mom_iq_c). Was ändert sich dadurch, verglichen mit dem unzentrierten Modell?

Zentrieren verändert die Regressionskoeffizienten (Steigungen) des Modells.Durch das Zentrieren verschwindet der Interaktionseffekt aus dem Modell.Zentrieren ist ausschließlich bei nominal skalierten Prädiktoren wie mom_hs sinnvoll, nicht bei metrischen wie mom_iq.Zentrieren verändert die vorhergesagten Werte des Modells für die beobachteten Datenpunkte.Mit zentriertem Prädiktor bezieht sich der Achsenabschnitt auf ein Kind einer Mutter mit durchschnittlichem IQ (und ohne Schulabschluss); ohne Zentrierung bezöge er sich auf die kaum sinnvoll interpretierbare Situation einer Mutter mit IQ\(=0\).

Das Kapitel zeigt explizit: “Mit zentrierten Prädiktoren gibt der Achsenabschnitt den Y-Wert für eine Beobachtung mit mittlerem X-Wert an; daher ist der Achsenabschnitt besser zu interpretieren als ohne Zentrierung.” Ohne Zentrierung müsste man den Achsenabschnitt auf eine (unrealistische) Mutter mit IQ von 0 beziehen. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option widerspricht dem Kapitel direkt: “Zentrieren ändert nichts an den Vorhersagen” – die Vergleichsrechnung zwischen zentriertem und unzentriertem Modell zeigt (bis auf Rundungsfehler) identische vorhergesagte Werte.

• Die dritte Option widerspricht ebenfalls dem Kapitel: “Das Zentrieren ändert auch nicht die Regressionskoeffizienten, da die Streuungen dieser Variablen nicht verändert wurden durch das Zentrieren” – Zentrieren ist eine reine Verschiebung, keine Reskalierung.

• Die vierte Option ist falsch: Der Interaktionseffekt bleibt im zentrierten Modell (m10.5) genauso enthalten wie im unzentrierten Modell (m10.4), nur die Interpretation des Achsenabschnitts (und der Haupteffekte bei Referenzwert Null) ändert sich.

• Die fünfte Option kehrt die im Kapitel gegebene Empfehlung um: Dort wird explizit empfohlen, metrische Prädiktoren (wie mom_iq) zu zentrieren; bei nominalen Prädiktoren wie mom_hs ist eine Zentrierung hingegen für eine einfache Interpretation meist nicht sinnvoll.

• Falsch

• Falsch

• Falsch

• Falsch

• Wahr

Warum werden im Kapitel die Prädiktoren mom_iq und mom_age z-standardisiert, bevor ihre “Relevanz” für die Vorhersage von kid_score verglichen wird?

Z-standardisierte Variablen erlauben automatisch eine kausale Interpretation der Effekte, unstandardisierte hingegen nicht.Z-Standardisierung verändert, ob ein Nulleffekt anhand des 95%-Schätzbereichs ausgeschlossen werden kann oder nicht.Unstandardisierte Prädiktoren können von stan_glm() technisch gar nicht verarbeitet werden.Der Vergleich der unstandardisierten Regressionskoeffizienten kann irreführend sein, wenn die Prädiktoren unterschiedliche Streuungen haben; die z-Standardisierung macht die Effekte in Einheiten der jeweiligen Streuung vergleichbar.Durch die z-Standardisierung verschwindet ein etwaiger Interaktionseffekt zwangsläufig aus dem Modell.

Das Kapitel stellt fest: “mom_iq hat den größeren Koeffizienten und viel Streuung; mom_age hat weniger Streuung.” Ein Vergleich der rohen (unstandardisierten) Koeffizienten wäre daher irreführend, da unterschiedliche Skalen/Streuungen die Größe der Koeffizienten mitbestimmen. Die z-Standardisierung bringt beide Prädiktoren auf eine gemeinsame Skala (Mittelwert 0, SD 1), sodass ihre Koeffizienten direkt in “SD-Einheiten” vergleichbar werden. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option ist falsch: stan_glm() verarbeitet unstandardisierte Prädiktoren problemlos, wie an den zahlreichen Beispielen mit unstandardisierten Variablen im Kapitel gezeigt wird – die Standardisierung ist eine Frage der Interpretierbarkeit, keine technische Notwendigkeit.

• Die dritte Option widerspricht der wiederholten Warnung des Kapitels: Auch standardisierte Regressionskoeffizienten sind zunächst nur statistische (assoziative), keine automatisch kausalen Effekte – die Standardisierung selbst rechtfertigt keine Kausalaussage.

• Die vierte Option ist falsch: Die Standardisierung der Prädiktoren ändert nichts an der grundsätzlichen Existenz eines Interaktionseffekts im Modell, sofern dieser weiterhin spezifiziert wird.

• Die fünfte Option ist falsch: Die z-Standardisierung ist eine reine (lineare) Reskalierung der Variablen und ändert daher nicht, ob ein Effekt statistisch von Null unterscheidbar ist – nur die Einheit, in der der Koeffizient ausgedrückt wird, ändert sich.

• Falsch

• Falsch

• Falsch

• Wahr

• Falsch

Im vollständig standardisierten Modell kid_score_z ~ mom_iq_z + mom_age_z (m10.12) wird der Koeffizient für mom_iq_z auf ca. 0,44 geschätzt. Wie ist dieser Wert korrekt zu interpretieren?

Der Koeffizient 0,44 kann nicht sinnvoll mit dem Koeffizienten für mom_age_z verglichen werden, da beide Prädiktoren unterschiedliche ursprüngliche Einheiten (IQ-Punkte vs. Jahre) haben.Der Wert 0,44 bedeutet, dass 44% der Kinder einen höheren IQ als ihre Mutter haben.0,44 ist der Anteil erklärter Varianz (\(R^2\)) des Gesamtmodells.Ändert sich der IQ der Mutter um eine Standardabweichung, ändert sich der IQ des Kindes im Schnitt um etwa 0,44 Standardabweichungen.Der IQ des Kindes ändert sich um 0,44 IQ-Punkte, wenn sich der IQ der Mutter um einen einzigen IQ-Punkt ändert.

Da sowohl AV als auch UV z-standardisiert sind (Mittelwert 0, SD 1), liefert der Regressionskoeffizient eine Aussage in “SD-Einheiten”: Eine Veränderung von mom_iq um eine Standardabweichung geht im Schnitt mit einer Veränderung von kid_score um etwa 0,44 Standardabweichungen einher. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option beschreibt die Interpretation eines unstandardisierten Koeffizienten (in IQ-Punkten pro IQ-Punkt) – im vollständig standardisierten Modell beziehen sich beide Variablen aber auf SD-Einheiten, nicht auf die ursprüngliche Punkteskala.

• Die dritte Option verwechselt den Regressionskoeffizienten eines einzelnen Prädiktors mit der Gesamt-Modellgüte \(R^2\), die eine andere (im Kapitel separat berichtete) Kennzahl ist.

• Die vierte Option interpretiert den Koeffizienten fälschlich als einen Anteilswert von Personen – ein standardisierter Regressionskoeffizient ist aber ein Zusammenhangsmaß, keine Prozentangabe über Personen.

• Die fünfte Option widerspricht gerade dem zentralen Vorteil der z-Standardisierung: Genau weil beide Prädiktoren auf eine gemeinsame Skala (SD-Einheiten) gebracht wurden, sind ihre Koeffizienten nun direkt miteinander vergleichbar – das Kapitel nutzt dies explizit, um zu zeigen, dass mütterlicher IQ “deutlich wichtiger” ist als mütterliches Alter.

• Falsch

• Falsch

• Falsch

• Wahr

• Falsch

Für das Modell kid_score ~ mom_hs gilt näherungsweise: \(\text{kid\_score} = 78 + 12\cdot\text{mom\_hs}\). Welchen mittleren IQ-Wert sagt das Modell für ein Kind vorher, dessen Mutter einen Schulabschluss hat (\(\text{mom\_hs}=1\))?

\(78\)\(12\)\(78 - 12 = 66\)\(78 + 12\cdot 1 = 90\)\(78 \cdot 12 = 936\)

Man setzt \(\text{mom\_hs}=1\) in die Modellgleichung ein: \(\text{kid\_score} = 78 + 12\cdot 1 = 90\).

Warum die Distraktoren falsch sind:

• Die zweite Option gibt den vorhergesagten Wert für die andere Gruppe zurück (\(\text{mom\_hs}=0\): \(78+12\cdot 0=78\)), nicht für die gefragte Gruppe mit Schulabschluss.

• Die dritte Option gibt nur das Regressionsgewicht (die Differenz zwischen den Gruppen) zurück, ohne den Achsenabschnitt (den Basiswert) zu addieren.

• Die vierte Option verwendet das falsche Vorzeichen (Subtraktion statt Addition) für den Term \(12\cdot\text{mom\_hs}\).

• Die fünfte Option multipliziert Achsenabschnitt und Regressionsgewicht miteinander, statt das Regressionsgewicht mit dem Prädiktorwert zu multiplizieren und zum Achsenabschnitt zu addieren – das widerspricht der linearen Modellgleichung.

• Falsch

• Falsch

• Falsch

• Wahr

• Falsch

Ein Regressionsmodell mit Interaktion ist gegeben durch \(\text{kid\_score} = \beta_0 + \beta_1\cdot\text{mom\_hs} + \beta_2\cdot\text{mom\_iq} + \beta_3\cdot\text{mom\_hs}\cdot\text{mom\_iq}\), mit den Punktschätzern \(\beta_0=20\), \(\beta_1=8\), \(\beta_2=0{,}5\), \(\beta_3=-0{,}1\). Welchen kid_score-Wert sagt das Modell für ein Kind vorher, dessen Mutter einen Schulabschluss hat (\(\text{mom\_hs}=1\)) und einen IQ von 100 hat?

\(20 + 8\cdot 1 + 0{,}5\cdot 100 + (-0{,}1)\cdot 1\cdot 100 = 20+8+50-10 = 68\)\(20 + 8\cdot 1 + 0{,}5\cdot 100 = 20+8+50 = 78\)\(20 + 0{,}5\cdot 100 + (-0{,}1)\cdot 1\cdot 100 = 20+50-10 = 60\)\(20 + 8\cdot 1 + 0{,}5\cdot 100 - (-0{,}1)\cdot 1\cdot 100 = 20+8+50+10 = 88\)\(20 + 8 + 0{,}5 + (-0{,}1) = 28{,}4\)

Man setzt alle gegebenen Werte in die vollständige Modellgleichung ein: \(20 + 8\cdot 1 + 0{,}5\cdot 100 + (-0{,}1)\cdot 1\cdot 100 = 20+8+50-10 = 68\).

Warum die Distraktoren falsch sind:

• Die zweite Option vergisst den Interaktionsterm \(\beta_3\cdot\text{mom\_hs}\cdot\text{mom\_iq}\) komplett – das entspricht der (falschen) Annahme, es gäbe keine Interaktion im Modell.

• Die dritte Option vergisst den Haupteffekt-Term \(\beta_1\cdot\text{mom\_hs}\), obwohl dieser Term im Modell explizit vorhanden ist und mit eingerechnet werden muss.

• Die vierte Option verwendet ein falsches Vorzeichen beim Interaktionsterm (Subtraktion des bereits negativen Werts \(\beta_3=-0{,}1\) statt dessen Addition), was den Term versehentlich positiv statt negativ werden lässt.

• Die fünfte Option setzt fälschlich alle Prädiktorwerte (mom_hs, mom_iq) auf 1 und addiert nur die Koeffizienten selbst, statt die tatsächlich gegebenen Werte (\(\text{mom\_hs}=1\), \(\text{mom\_iq}=100\)) korrekt in die Modellgleichung einzusetzen.

• Wahr

• Falsch

• Falsch

• Falsch

• Falsch

Was prüft das ROPE-Verfahren (Region of Practical Equivalence), wie es im Kapitel zur Prüfung der “Praktisch-Null-Hypothese” (z.B. beim Effekt des mütterlichen Schulabschlusses) verwendet wird?

Der ROPE-Bereich ist ein feststehender, gesetzlich vorgeschriebener Wertebereich und darf nicht vom Forschenden selbst festgelegt werden.ROPE prüft, wie stark sich die Post-Verteilung eines Effekts mit einem vorab definierten Bereich “praktisch vernachlässigbar kleiner” Werte um Null überlappt; ist die Überlappung gering oder gibt es keine, wird die Praktisch-Null-Hypothese verworfen.ROPE berechnet einen klassischen frequentistischen p-Wert für den Regressionskoeffizienten.ROPE ersetzt die gesamte Post-Verteilung durch eine einzige Punktschätzung des Effekts.ROPE prüft ausschließlich, ob der Effekt exakt Null ist, nicht ob er praktisch bedeutsam (relevant) ist.

Beim ROPE-Verfahren definiert man vorab einen Bereich “praktisch vernachlässigbarer” Effektgrößen (im Beispiel des Kapitels: \(\pm 5\) IQ-Punkte) und prüft, wie viel der Post-Verteilung des Effekts innerhalb dieses Bereichs liegt. Überlappt die Post-Verteilung kaum oder gar nicht mit diesem Bereich, spricht das dafür, dass der Effekt nicht nur statistisch von Null verschieden, sondern auch praktisch bedeutsam ist – die Praktisch-Null-Hypothese wird dann verworfen. Die erste Option ist korrekt.

Warum die Distraktoren falsch sind:

• Die zweite Option verwechselt das Bayes’sche ROPE-Verfahren mit einem frequentistischen p-Wert-Test – ROPE arbeitet direkt mit der Post-Verteilung (bzw. deren HDI) und deren Überlappung mit dem ROPE-Bereich, nicht mit einem p-Wert.

• Die dritte Option ist falsch: ROPE nutzt gerade die vollständige Post-Verteilung (bzw. deren HDI), um den Anteil der Überlappung mit dem ROPE-Bereich zu bestimmen – eine einzelne Punktschätzung würde diese Information verlieren.

• Die vierte Option widerspricht dem Kapitel direkt: Der Grenzwert für “praktisch gleich” bzw. “kaum unterschiedlich” wird explizit als subjektiv bezeichnet und muss von den Forschenden inhaltlich begründet festgelegt werden (im Beispiel: 5 IQ-Punkte).

• Die fünfte Option beschreibt eher die exakte Nullhypothese (Effekt \(=0\)), die die klassische Signifikanztestung prüft – der Witz am ROPE-Verfahren ist gerade, dass es über die exakte Null hinausgeht und einen ganzen Bereich “praktisch bedeutungsloser” Werte berücksichtigt.

• Falsch

• Wahr

• Falsch

• Falsch

• Falsch

12.16 —

Cohen, J. (1992). A Power Primer. Psychological Bulletin, 112(1), 155–159.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and Other Stories. Cambridge University Press.

Goodrich, B., Gabry, J., Ali, I., & Brilleman, S. (2020). Rstanarm: Bayesian Applied Regression Modeling via Stan. https://mc-stan.org/rstanarm

McElreath, R. (2020). Statistical Rethinking: A Bayesian Course with Examples in R and Stan (2. Aufl.). Taylor and Francis, CRC Press.

---

1. https://statistik1.netlify.app/080-regression1

2. https://statistik1.netlify.app/090-regression2

3. https://raw.githubusercontent.com/sebastiansauer/Lehre/main/data/kidiq.csv

4. Hier ist eine kurze Erklärung dazu: https://statistik1.netlify.app/010-rahmen#sec-arten-variablen

5. Hier finden Sie eine kurze Erklärung zur Interaktion: https://statistik1.netlify.app/090-regression2#interaktion

6. auch DAG genannt, s. Kapitel 13

7. Häufig erlaubt uns unser Vorwissen eine gerichtete Hypothese – “größer als/kleiner als” – zu formulieren, anstelle der “empirisch ärmeren” einfachen, ungerichteten Ungleichheit

8. unknown, sozusagen der unbekannte Gott, also für alle sonstigen Einflüsse; man kann das “u” ohne Schaden weglassen, da wir es sowieso nicht modellieren. Hier ist es nur aufgeführt, um zu verdeutlichen, dass wir nicht so verwegen sind, zu behaupten, es gäbe keine anderen Einflüsse als mom_hs auf die IQ des Kindes.

9. Genauer gesagt wird geprüft, wie wahrscheinlich es auf Basis des Modells ist, noch extremere Ergebnisse zu beobachten unter der Annahme, dass die (exakte) Nullhypothese wahr ist. Es ist etwas kompliziert.

10. Genauer gesagt, erlaubt der t-Test in der Form des Welch-Tests auch Abweichungen von der Varianzhomogenität der getesteten Gruppen. Die Regression geht hingegen von Varianzhomogenität (Homoskedastizität) aus. Allerdings ist diese Annahme nicht von besonderer Bedeutung, wenn es um die Regressionskoeffizienten geht.

11. DAG, s. Kapitel 13

12. Dabei nehmen wir an, dass x und u nicht voneinander abhängen, was man daran erkennt, dass es keine Pfeile zwischen den beiden Variablen gibt.

13. Vgl. Abschnitt “UV zentrieren” (https://statistik1.netlify.app/090-regression2.html#uv-zentrieren) im Kursbuch Statistik1 (https://statistik1.netlify.app/).

14. https://vincentarelbundock.github.io/Rdatasets/csv/palmerpenguins/penguins.csv

15. https://vincentarelbundock.github.io/Rdatasets/doc/palmerpenguins/penguins.html

16. https://jfly.uni-koeln.de/color/

17. s. Details unter https://mc-stan.org/rstanarm/articles/priors.html#how-to-specify-flat-priors-and-why-you-typically-shouldnt.

18. https://xkcd.com/552/

19. https://plotly.com/r/

20. am nützlichsten ist diese Standardisierung bei normal verteilten Variablen.

21. Zumindest zur Zeit als ich diese Zeilen schreibe. Achtung: Voreinstellungen können sich ändern. Am besten in der Dokumentation nachlesen: ?parameters.

22. Nicht unbedingt die feine bayesianische Art, denn die Prioris sollten ja eigentlich apriori, also vor Kenntnis der Daten, bestimmt werden. Auf der anderen Seite behauptet Stan, von uns zur Rede gestellt, dass die empirischen Mittelwerte ja doch gute Schätzer der echten Parameter sein müssten, wenn die Stichprobe, die wir ihm angeschleppt hätten, tatsächlich gut ist…

23. 4: \(\beta_0, \beta_1, \beta_2, \sigma\)

24. Ignorieren Sie die Zeile mit dem Befehl display(). Dieser Befehl dient nur dazu, die Ausgabe zu verschönern in Markdown-Dokumenten, wie im Quelltext dieses Kapitels.

25. was auf normal verteilte AV hinauslaufen kann aber nicht muss

26. https://giphy.com/gifs/XIqCQx02E1U9W