8  Einfache lineare Modelle

8.1 Lernsteuerung

8.1.1 Position im Modulverlauf

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

8.1.2 Lernziele

Nach Absolvieren des jeweiligen Kapitels sollen folgende Lernziele erreicht sein.

Sie können …

• die Post-Verteilung für einfache lineare Modelle in R berechnen
• zentrale Informationen zu Modellparametern - wie Lage- oder Streuungsmaße und auch Schätzintervalle - aus der Post-Verteilung herauslesen
• künftige, laut Modell zu erwartende Beobachtungen mit der PPV simulieren

8.1.3 Begleitliteratur

Der Stoff dieses Kapitels orientiert sich an McElreath (2020), Kap. 4.4.

8.1.4 Vorbereitung im Eigenstudium

• Statistik1, Kap. “Geradenmodelle 1”

8.1.5 Benötigte R-Pakete

In diesem Kapitel benötigen Sie folgende R-Pakete.

library(tidyverse)
library(easystats)
library(rstanarm)  # Bayes-Golem
library(ggpubr)  # Datenvisualisierung

Da wir in diesem Kapitel immer mal wieder eine Funktion aus dem R-Paket {easystats} verwenden: Hier finden Sie eine Übersicht aller Funktionen des Pakets.¹

8.1.6 Benötigte Daten

In diesem Kapitel benötigen wir den Datensatz zu den !Kung-Leuten, Howell1a, McElreath (2020). Sie können ihn hier herunterladen.

Kung_path <- "https://raw.githubusercontent.com/sebastiansauer/Lehre/main/data/Howell1a.csv"

d <- read.csv(Kung_path)

d2 <- d %>% filter(age > 18)

1

Pfad zum Datensatz; Sie müssen online sein, um die Daten herunterzuladen.

2

Daten einlesen

3

Auf Erwachsene Personen begrenzen (d.h. Alter > 18)

8.1.7 Einstieg

Beispiel 8.1 (Grundkonzepte der linearen Regression) Fassen Sie die Grundkonzepte der linearen Regression kurz zusammen! \(\square\)

Beispiel 8.2 (Was ist eine Post-Verteilung und wozu ist sie gut?) Erklären Sie kurz, was eine Post-Verteilung ist - insbesondere im Zusammenhang mit den Koeffizienten einer einfachen Regression - und wozu sie gut ist. \(\square\)

8.1.8 Überblick

Dieses Kapitel stellt ein einfaches Regressionsmodell vor, wo die Körpergröße auf das Gewicht zurückgeführt wird; also ein sehr eingängiges Modell.
Neu ist dabei lediglich, dass die Parameter des Modells - \(\beta_0\), \(\beta_1\), \(\sigma\) - jetzt über eine Post-Verteilung verfügen. Die Post-Verteilung ist der Zusatznutzen der Bayes-Statistik. Die “normale” Regression hat uns nur einzelne Werte für die Modellparameter geliefert (“Punktschätzer”). Mit Bayes haben wir eine ganz Verteilung pro Parameter.

8.2 Post-Verteilung der Regression

8.2.1 Einfache Regression

Die (einfache) Regression prüft, inwieweit zwei Variablen, \(Y\) und \(X\) linear zusammenhängen. Je mehr sie zusammenhängen, desto besser kann man \(X\) nutzen, um \(Y\) vorherzusagen (und umgekehrt). Hängen \(X\) und \(Y\) zusammen, heißt das nicht (unbedingt), dass es einen kausalen Zusammenhang zwischen \(X\) und \(Y\) gibt. Linear ist ein Zusammenhang, wenn der Zuwachs in \(Y\) relativ zu \(X\) konstant ist: wenn \(X\) um eine Einheit steigt, steigt \(Y\) immer um \(b\) Einheiten (nicht kausal, sondern deskriptiv gemeint).²

Laden wir die !Kung-Daten und visualisieren wir uns den Zusammenhang zwischen Gewicht (X) und Größe (Y), Abbildung 8.1.

d2 %>% 
  ggplot(
       aes(x = weight, y = height)) +
  geom_point(alpha = .7) +
  geom_smooth(method = "lm")

ggscatter(d2,
          x = "weight", y = "height",
          add = "reg.line") 

Mit ggplot2

Mit ggpubr

Abbildung 8.1: Der Zusammenhang zwischen Gewicht (X) und Größe (Y)

8.2.2 Statistiken zum !Kung-Datensatz

Die Daten können Sie hier herunterladen.

Tabelle 8.1 zeigt die zentralen deskriptiven Statistiken zum !Kung-Datensatz.

Kung_path <- "data/Howell1a.csv"  
d <- read_csv(Kung_path)  

d2 <- d %>% filter(age > 18)

describe_distribution(d2)

Tabelle 8.1: Verteiung der (metrischen) Variablen im !Kung-Datensatz

Variable	Mean	SD	IQR	Range	Skewness	Kurtosis	n	n_Missing
height	154.64	7.77	12.06	(136.53, 179.07)	0.14	-0.50	346	0
weight	45.05	6.46	9.14	(31.52, 62.99)	0.14	-0.53	346	0
age	41.54	15.81	22.00	(19.00, 88.00)	0.68	-0.20	346	0
male	0.47	0.50	1.00	(0.00, 1.00)	0.10	-2.00	346	0

Wie aus Tabelle 8.1 abzulesen ist, beträgt das mittlere Körpergewicht (weight) liegt ca. 45kg (sd 7 kg).

8.2.3 Etwas mehr EDA

Wir brauchen die EDA hier nicht wirklich, aber es ist praktisch. Das Paket DataExplorer hat ein paar nette Hilfen zur explorativen Datenanalyse.

library(DataExplorer)

8.2.3.1 Gibt es fehlende Werte?

Nein, s. Abb. Abbildung 8.2.

d2 %>% plot_missing()

Abbildung 8.2: Fehlende Werte - fehlen.

8.2.3.2 Verteilung der numerischen Variablen

Betrachten wir die Verteilung der numerischen Variablen des Datensatzes, s. Abbildung 8.3.

d2 %>% plot_histogram()

Abbildung 8.3: Verteilung (als Histogramme dargestellt) der numerischen Variablen des Datensatzes

8.2.3.3 Verteilung der kategorialen Variablen

Betrachten wir die Verteilung der kategorialen Variablen des Datensatzes, s. Abbildung 8.4.

d2 %>% plot_bar()

Abbildung 8.4: Verteilung (als Balkendiagramme dargestellt) der kategorialen Variablen des Datensatzes

8.2.3.4 Korrelationen

Die Korrelationen der (numerischen) Variablen sind in Abbildung 8.5 dargestellt.

d2 %>% plot_correlation()

Abbildung 8.5: Korrelationsmatrix

Übungsaufgabe 8.1 (EDA-Bericht) Probieren Sie mal die folgende Funktion aus, die Ihnen einen Bericht zur EDA erstellt: create_report(d2). \(\square\)

8.2.4 Prädiktor zentrieren

Zieht man von jedem Gewichtswert den Mittelwert ab, so bekommt man die Abweichung des Gewichts vom Mittelwert (Prädiktor “zentrieren”, engl. to center). Wenn man den Prädiktor (weight) zentriert hat, ist der Achsenabschnitt, \(\beta_0\), einfacher zu verstehen. In einem Modell mit zentriertem Prädiktor (weight) gibt der Achsenabschnitt die Größe einer Person mit durchschnittlichem Gewicht an. Würde man weight nicht zentrieren, gibt der Achsenabschnitt die Größe einer Person mit weight=0 an, was nicht wirklich sinnvoll zu interpretieren ist. Vgl. Gelman et al. (2021), Kap. 10.4, 12.2.

Man zentriert eine Variable \(X\), indem man von \(x_i\) den Mittelwert \(\bar{x}\) abzieht: \(x_i - \b ar{x}\).

d3 <-
  d2 %>% 
  mutate(weight_c = weight - mean(weight))

Mit Hilfe von center() aus {easystats} kann man sich das Zentrieren auch erleichtern.

d3 <- 
  d2 %>% 
  mutate(weight_c = as.numeric(center(weight)))

height	weight	age	male	weight_c
152	48	63	1	3
140	36	63	0	−9
137	32	65	0	−13

Wie man sieht, wird die Verteilung von weight durch die Zentrierung “zur Seite geschoben”: Der Mittelwert von weight_c (das zentrierte Gewicht) liegt jetzt bei 0, s. Abbildung 8.6.

Abbildung 8.6: Das Zentrieren ändert die Verteilungsform nicht, sondern “schiebt” die Verteilung nur zur Seite

Das schwierigste ist dabei, nicht zu vergessen, dass d3 die Tabelle mit zentriertem Prädiktor ist, nicht d2.

8.3 Modell m43: zentrierter Prädiktor

📺 Prädiktoren zentrieren

Einige Regressionskoeffizienten, wie der Achsenabschnitt (Intercept) sind schwer zu interpretieren: Bei einem (erwachsenen) Menschen mit Gewicht 0, was wäre wohl die Körpergröße? Hm, Philosophie steht heute nicht auf der Tagesordnung.

Da wäre es schön, wenn wir die Daten so umformen könnten, dass der Achsenabschnitt eine sinnvolle Aussage macht. Zum Glück geht das leicht: Wir zentrieren den Prädiktor (Gewicht)!

Wichtig

Durch Zentrieren kann man die Ergebnisse einer Regression einfacher interpretieren.

8.3.1 Modelldefinition von m43

Für jede Ausprägung des Prädiktors (weight_centered), \(wc_i\), wird eine Post-Verteilung für die abhängige Variable (height, \(h_i\)) berechnet. Der Mittelwert \(\mu\) für jede Post-Verteilung ergibt sich aus dem linearen Modell (unserer Regressionsformel). Die Post-Verteilung berechnet sich auf Basis der Priori-Werte und des Likelihood (Bayes-Formel). Wir brauchen Priori-Werte für die Steigung \(\beta_1\) und den Achsenabschnitt \(\beta_0\) der Regressionsgeraden. Außerdem brauchen wir einen Priori-Wert, der die Streuung \(\sigma\) der Größe (height) angibt; dieser Wert wird als exonentialverteilt angenommen. Der Likelihood gibt an, wie wahrscheinlich ein Wert height ist, gegeben \(\mu\) und \(\sigma\). Gleichung 8.1 stellt die Modelldefinition dar.

\[\begin{align*} \color{red}{\text{height}_i} & \color{red}\sim \color{red}{\operatorname{Normal}(\mu_i, \sigma)} && \color{red}{\text{Likelihood}} \\ \color{green}{\mu_i} & \color{green}= \color{green}{\beta_0 + \beta_1\cdot \text{weightcentered}_i} && \color{green}{\text{Lineares Modell} } \\ \color{blue}{\beta_0} & \color{blue}\sim \color{blue}{\operatorname{Normal}(178, 20)} && \color{blue}{\text{Priori}} \\ \color{blue}{\beta_1} & \color{blue}\sim \color{blue}{\operatorname{Normal}(0, 10)} && \color{blue}{\text{Priori}}\\ \color{blue}\sigma & \color{blue}\sim \color{blue}{\operatorname{Exp}(0.1)} && \color{blue}{\text{Priori}} \end{align*} \tag{8.1}\]

Hinweis

Der Achsenabschnitt (engl. intercept) eines Regressionsmodell wird in der Literatur oft mit \(\beta_0\) bezeichnet, aber manchmal auch mit \(\alpha\). Und manchmal mit noch anderen Buchstaben, das Alphabet ist weit. 🤷

8.3.2 Likelihood, m43

\[ \begin{aligned} \color{red}{\text{height}_i} & \color{red}\sim \color{red}{\operatorname{Normal}(\mu_i, \sigma)} && \color{red}{\text{Likelihood}} \end{aligned} \]

Der Likelihood von m43 ist ähnlich zu den vorherigen Modellen (m41, m42). Nur gibt es jetzt ein kleines “Index-i” am \(\mu\) und am \(h\) (h wie heights). Es gibt jetzt nicht mehr nur einen Mittelwert \(\mu\), sondern für jede Beobachtung (Zeile) einen Mittelwert \(\mu_i\). Lies etwa so:

“Die Wahrscheinlichkeit, eine bestimmte Größe bei Person \(i\) zu beobachten, gegeben \(\mu\) und \(\sigma\) ist normalverteilt (mit Mittelwert \(\mu\) und Streuung \(\sigma\))”.

8.3.3 Regressionsformel, m43

\[ \begin{aligned} \color{green}{\mu_i} & \color{green}= \color{green}{\beta_0 + \beta_1\cdot \text{weightcentered}_i} && \color{green}{\text{Lineares Modell} } \\ \end{aligned} \]

\(\mu\) ist jetzt nicht mehr ein Parameter, der (stochastisch) geschätzt werden muss. \(\mu\) wird jetzt (deterministisch) berechnet. Gegeben \(\beta_0\) und \(\beta_1\) ist \(\mu\) ohne Ungewissheit bekannt. \(\text{weight}_i\) ist der Prädiktorwert (weight) der \(i\)ten Beobachtung, also einer !Kung-Person (Zeile \(i\) im Datensatz). Lies etwa so:

“Der Mittelwert \(\mu_i\) der \(i\)ten Person berechnet sich als Summe von \(\beta_0\) und \(\beta_1\) mal \(\text{weight}_i\)”.

\(\mu_i\) ist eine lineare Funktion von weight. \(\beta_1\) gibt den Unterschied in height zweier Beobachtung an, die sich um eine Einheit in weight unterscheiden (Steigung der Regressionsgeraden). \(\beta_0\) gibt an, wie groß \(\mu\) ist, wenn weight Null ist (Achsenabschnitt, engl. intercept).

8.3.4 Priori-Werte des Modells m43

\[\begin{align*} \color{blue}\beta_1 & \color{blue}\sim \color{blue}{\operatorname{Normal}(178, 20)} && \color{blue}{\text{Priori Achsenabschnitt}} \\ \color{blue}\beta_1 & \color{blue}\sim \color{blue}{\operatorname{Normal}(0, 10)} && \color{blue}{\text{Priori Regressionsgewicht}}\\ \color{blue}\sigma & \color{blue}\sim \color{blue}{\operatorname{Exp}(0.1)} && \color{blue}{\text{Priori Sigma}} \end{align*}\]

Parameter sind hypothetische Kreaturen: Man kann sie nicht beobachten, sie existieren nicht wirklich. Ihre Verteilungen nennt man Priori-Verteilungen. \(\beta_0\) wurde in m41 als \(\mu\) bezeichnet, da wir dort eine “Regression ohne Prädiktoren” berechnet haben. \(\sigma\) ist uns schon als Parameter bekannt und behält seine Bedeutung aus dem letzten Kapitel. Da height nicht zentriert ist, der Mittelwert von \(\beta_0\) bei 178 und nicht 0. \(\beta_1\) fasst unser Vorwissen, ob und wie sehr der Zusammenhang zwischen Gewicht und Größe positiv (gleichsinnig) ist. Die Anzahl der Prioris entspricht der Anzahl der Parameter des Modells.

8.4 Die Post-Verteilung befragen

📺 Post-Verteilung auslesen 1

📺 Post-Verteilung auslesen 2

8.4.1 m43a

Sagen wir, auf Basis gut geprüfter Evidenz haben wir folgendes Modell festgelegt: height ~ weight_c, s. Gleichung 8.2.

Prioris:

\[\beta_1 \sim N(5,3); \\ \beta_0 \sim N(178, 20); \\ \sigma \sim E(0.1) \tag{8.2}\]

Wir nennen das Modell m43a³, s. Listing 8.1.

Listing 8.1: Modelldefinition von m43a in R

m43a <-
  stan_glm(
    height ~ weight_c,  # Regressionsformel
    prior = normal(5, 3),  # Regressionsgewicht (beta 1)
    prior_intercept = normal(178, 20),  # beta 0
    prior_aux = exponential(0.1),  # sigma
    refresh = 0,  # zeig mir keine Details
    seed = 42,  # lege die Zufallszahlen fest für Reproduzierbarkeit
    data = d3)

Hinweis

Mit seed kann man die Zufallszahlen fixieren, so dass jedes Mal die gleichen Werte resultieren. So ist die Nachprüfbarkeit der Ergebnisse (“Reproduzierbarkeit”) sichergestellt⁴. Welche Wert für seed man verwendet, ist egal, solange alle den gleichen verwenden. Der Autor verwendet z.B. oft den Wert 42. Zur Erinnerung: Der Golem zieht Zufallszahlen, damit erstellt er Stichproben, die die Postverteilung schätzen.

8.4.2 Mittelwerte von \(\beta_0\) und \(\beta_1\) aus der Post-Verteilung

Die ersten paar Zeilen:

id	(Intercept)	weight_c	sigma
1	155.1	0.9	5.0
2	155.5	0.8	5.1
3	155.5	0.9	5.1

Hier sind die Zusammenfassungen der Stichproben aus der Post-Verteilung, komfortabel zu erhalten mit dem Befehle parameters, s. Tabelle 8.2.

Tabelle 8.2: Parameter von m43a

parameters(m43a)

Parameter	Median	95% CI	pd	Rhat	ESS	Prior
(Intercept)	154.65	(154.14, 155.19)	100%	0.999	3214.00	Normal (178 +- 20)
weight_c	0.91	(0.82, 0.99)	100%	1.001	4134.00	Normal (5 +- 3)

Definition 8.1 (Effektwahrscheinlichkeit) Die Kennzahl pd (propability of direction) gibt die Effektwahrscheinlichkeit an: Die Wahrscheinlichkeit, dass der Effekt positiv (also größer als Null) oder negativ ist (je nachdem ob der Median des Effekts positiv oder negativ ist). pd gibt aber nicht an, wie stark der Effekt ist, nur ob er klar auf einer Seite der Null liegt. Damit ist er so etwas (grob!) Ähnliches wie der p-Wert in der Frequentistischen Statistik (Makowski et al., 2019).

Am besten das Diagramm dazu anschauen, s Abbildung 8.7.

plot(p_direction(m43a))

Abbildung 8.7: Diagramm zur Probability of Direction, Modell m43a

Rhat und ESS sind Kennzahlen, die untersuchen, ob mit der Stichprobenziehung im Bayes-Modell alles gut funktioniert hat. Bei einfachen Modellen (die wir hier berechnen) sollte da in der Regel alles in Ordnung sein. Rhat sollte nicht (viel) größer als 1 oder 1,01 sein. ESS (effective sample size) gibt die Anzahl der effektiv nutzbaren Stichproben an (im Standard werden 4000 berechnet). Die Zahl sollte nicht deutlich geringer sein.

Wir werden uns aber mit diesen beiden Kennwerten nicht weiter beschäftigen in diesem Kurs.

8.4.3 Visualisieren der “mittleren” Regressiongeraden

Zur Erinnerung: Die Bayes-Analyse liefert uns viele Stichproben zu den gesuchten Parametern, hier \(\beta_0\), \(\beta_1\) und \(\sigma\). Überzeugen wir uns mit einem Blick in die Post-Verteilung von m43a:

m43a %>% 
  as_tibble() %>% 
  head()

Wir können z.B. ein Lagemaß wie den Median hernehmen, um die “mittlere” Regressionsgerade zu betrachten.

Mit ggplot

d3 %>% 
  ggplot() +
  aes(x = weight_c, y = height) +
  geom_point() +
  geom_abline(
    slope = 0.9,  # Median beta 1
    intercept = 154,  # Median beta 0
    color = "blue")

Mit easystats

Einfacher ist die Syntax vielleicht, wenn man die Funktion estimate_expectation benutzt, s. Abbildung 8.8. Mit “expectation” sind hier die erwarteten Werte, also die Regressionsgerade, gemeint.

m43_expect <- estimate_expectation(m43a)   # aus {easystats}
plot(m43_expect)

Abbildung 8.8: Erwartete Werte des Modell m43a, sprich, die Regressionsgerade

8.4.4 Zentrale Statistiken zu den Parametern

In diesem Modell gibt es drei Parameter: \(\beta_0, \beta_1, \sigma\).⁵ Hier folgen einige Beispiele an Fragen, die wir an unser Modell bzw. die Post-Verteilung stellen können.

8.4.4.1 Lagemaße zu den Parametern

• Was ist die mittlere Größe einer !Kung-Person? (\(\beta_0\))
• Was ist der Schätzwert für den Zusammenhang von Gewicht und Größe? (\(\beta_1\))
• Was ist der Schätzwert für Ungewissheit in der Schätzung der Größe? (\(\sigma\))
• Was ist der wahrscheinlichste Wert für z.B: \(\beta_1\)?

Eine nützliche Zusammenfassung der Post-Verteilung bekommt man mit parameters(modell), s. Tabelle 8.2.

Wandelt man das Ausgabe-Objekt der Bayes-Regression, d.h. m43a, mit as_tibble() in eine Tabelle um, so bekommt man eine Tabelle mit den Stichproben der Post-Verteilung:

m43a_post <- 
  m43a %>% 
  as_tibble()

m43a_post %>% 
  head()

Wie wir gesehen haben, nutzen wir diese Tabelle der Post-Verteilung immer wieder. Speichern wir uns sie also als ein Objekt ab, m43_post. Jetzt haben wir wieder eine schöne Tabelle mit Stichproben aus der Post-Verteilung, die wir wie gewohnt befragen können. Eine Visualisierung zeigt gut sowohl Lage- als auch Streuungsmaße der Parameter, zumindest grob.,

Oder man erstellt selber ein Diagramm mit ggplot oder ggpubr, s. Abbildung 8.9.

m43a_post %>% 
  ggplot(aes(x = weight_c)) +
  geom_density(fill = "orange")

Abbildung 8.9: Postverteilung für den Parameter Gewicht (zentriert)

Abbildung 8.9 zeigt, dass Mittelwert, Median und Modus eng zusammenliegen. Zur Erinnerung: Der Modus gibt den häufigsten, d.h. hier also den wahrscheinlichsten, Wert an. Der Modus wird hier auch Maximum a Posteriori (MAP) genannt, daher:

m43a_post %>% 
  summarise(map_b1 = map_estimate(weight_c))

Hier ist die Verteilung von \(\sigma\) visualisiert, s. Abbildung 8.10.

m43a_post %>% 
  ggplot(aes(x = sigma)) +
  geom_density(fill = "orange")

Abbildung 8.10: Die Post-Verteilung für den Parameter sigma, m43a

Alternativ kann man sich die Verteilung eines Parameters auch so ausgeben lassen, gleich mit Intervallgrenzen, z.B. 95%, s. Abbildung 8.11.

m43a_hdi <- hdi(m43a_post)  # analog mit eti(m43a)

plot(m43a_hdi)

Abbildung 8.11: Die Parameter Gewicht (zentriert) und sigma des Modells m43a

Ergänzt man bei plot() noch show_intercept = TRUE wird auch der Achsenabschnitt angezeigt.

8.4.5 Streuungsmaße zu den Parametern

• Wie unsicher sind wir uns in den Schätzungen der Parameter?

Diese Frage wird durch die Ungewissheitsintervalle in der Ausgabe beantwortet.

Hinweis

An einigen Stellen wird empfohlen, anstelle eines (gebräuchlichen) 95%-Intervalls auf ein 90%- oder 89%-Intervall auszuweichen, aufgrund der besseren numerischen Stabilität.

8.4.6 Ungewissheit von \(\beta_0\) und \(\beta_1\) aus der Post-Verteilung visualisiert

10

Die ersten 10 Stichproben:

d3 %>% 
  ggplot(aes(x = weight_c, 
             y = height)) +
  geom_point() +
  geom_abline(
    data = m43a_post %>% 
      slice_head(n = 10),
    aes(slope = weight_c,
        intercept = `(Intercept)`),
    alpha = .3)

8.5 100

Die ersten 100 Stichproben:

d3 %>% 
  ggplot(aes(x = weight_c, 
             y = height)) +
  geom_point() +
  geom_abline(
    data = m43a_post %>% 
      slice_head(n = 100),
     aes(slope = weight_c,
        intercept = `(Intercept)`),
    alpha = .1)

1000

Die ersten 1e3 Stichproben:

d3 %>% 
  ggplot(aes(x = weight_c, 
             y = height)) +
  geom_point() +
  geom_abline(
    data = m43a_post %>% 
      slice_head(n = 1e3),
     aes(slope = weight_c,
        intercept = `(Intercept)`),
    alpha = .01)

1000000

Die ersten 1000000⁶ … okay, lassen wir es gut sein⁷.

Einfacher ist die Visualisierung mit estimate_expectation, s. Abbildung 8.12.

estimate_expectation(m43a, seed = 42) %>% plot()

Abbildung 8.12: Schätzbereich für die bedingten mittleren Körpergröße, also die Regressionsgerade mit Unsicherheitsintervall

8.5.1 Fragen zu Quantilen des Achsenabschnitts

Hinweis

Zur Erinnerung: Bei einem zentrierten Prädiktor misst der Achsenabschnitt die mittlere Größe⁸.

• Welche mittlere Größe wird mit einer Wahrscheinlichkeit von 50%, 90% bzw. 95% Wahrscheinlichkeit nicht überschritten?
• Welche mittlere Größe mit Wahrscheinlichkeit von 95% nicht unterschritten?
• Von wo bis wo reicht der innere 50%-Schätzbereich der mittleren Größe?

Quantile:

m43a_post %>% 
  summarise(
    q_50 = quantile(`(Intercept)`, prob = .5),
    q_90 = quantile(`(Intercept)`, prob = .9),
    q_05 = quantile(`(Intercept)`, prob = .95))

50%-PI:

m43a %>% 
  eti(ci = .5)

8.5.2 Fragen zu Wahrscheinlichkeitsmassen des Achsenabschnitts

Wie wahrscheinlich ist es, dass die mittlere Größe bei mind. 155 cm liegt?

m43a_post %>% 
  count(gross = `(Intercept)` >= 155) %>% 
  mutate(prop = n / sum(n))

Die Wahrscheinlichkeit beträgt 0.1.

Wie wahrscheinlich ist es, dass die mittlere Größe höchstens 154.5 cm beträgt?

m43a_post %>% 
  count(klein = (`(Intercept)` <= 154.5)) %>% 
  mutate(prop = n / sum(n))

Die Wahrscheinlichkeit beträgt 0.29.

8.5.3 Typischer Bayes-Nutzer in Aktion

Typischer Bayes-Nutzer, der ein Ungewissheitsintervall berechnet. Bildquelle: Easystats, bayestestR

Quelle

8.6 Post-Verteilung bedingt auf einen Prädiktorwert

8.6.1 Bei jedem Prädiktorwert eine Post-Verteilung für \(\mu\)

Komfort pur: Unser Modell erlaubt uns für jeden beliebigen Wert des Prädiktors eine Post-Verteilung (von \(\mu\)) zu berechnen.

Hier am Beispiel von m42, s. Abbildung 8.13.

Abbildung 8.13: Für jeden beliebigen Prädiktorwert kann man eine Post-Verteilung bekommen. A: Regressionsmodell mit einigen ausgewählten Gewichtswerten. Es ist jeweils die Wahrscheinlichkeitsverteilung für den vorhergesagten Y-Wert dargestellt (hier sind die Verteilungen zu groß dargestellt zur besseren Sichtbarkeit). B: Für jeden beliebigen Gewichtswert (Y) bekommt man eine (auf den jeweiligen X-Wert bedingten) Post-Verteilung.

8.6.2 Visualisierung

Was ist wohl die Wahrscheinlichkeit der Körpergröße bei einem bestimmten Gewicht?

Angenommen wir wissen, dass das Gewicht bei, sagen wir 45 kg liegt. Welche Körpergröße ist (im Schnitt) zu erwarten? Wie unsicher sind wir uns über diesen Mittelwert?

Etwas formaler ausgedrückt:

\(\mu|\text{weight}=45\)

45 kg entspricht genau dem Mittelwert von weight. Geht man von zentrierten Prädiktorwerten aus, gilt in dem Fall weight_c = 0. Erstellen wir uns dazu eine Tabelle:

mu_at_45 <-
  m43a_post %>% 
  mutate(mu_at_45 = `(Intercept)`)

Und plotten diese, s. Abbildung 8.14.

mu_at_45 %>% 
  ggplot(aes(x = mu_at_45)) +
  geom_density()

Abbildung 8.14: Post-Verteilung der Größe (laut unserem Modell) bei einem Gewicht von 45kg

Analog können wir fragen, wie groß wohl eine Person mit 50 kg im Mittelwert sein wird und wie (un)gewiss wir uns über diesen Mittelwert sind.

50 kg, das sind 5 über dem Mittelwert, in zentrierten Einheiten ausgedrückt also weight_c = 5. Auch dazu erstellen wir uns eine Tabelle, s. Tabelle 8.3.

mu_at_50 <-
  mu_at_45 %>% 
  mutate(mu_at_50 = `(Intercept)` + 5 * weight_c)

head(mu_at_50)

Tabelle 8.3: Die Verteilung von mu bedingt auf ein Gewicht von 50kg.

Die Verteilung der mittleren Größe bei einem Gewicht von 50kg ist weiter “rechts” (Richtung höhere Größe) zentriert, s. Abbildung 8.15.

mu_at_50 %>% 
  ggplot(aes(x = mu_at_50)) +
  geom_density()

Abbildung 8.15: Post-Verteilung der mittleren Größe (laut unserem Modell) bedingt auf ein Gewicht von 50 kg

8.6.3 Lagemaße und Streuungen

Befragen wir die bedingte Post-Verteilung. Eine erste Frage zielt nach den typischen deskriptiven Statistiken, also nach Lage und Streuung der Verteilung der Körpergröße.

Was ist das 90% PI für \(\mu|w=50\) ?

mu_at_50 %>% 
  eti(mu_at_50, ci = .9)

Die mittlere Größe - gegeben \(w=50\) - liegt mit 90% Wahrscheinlichkeit zwischen den beiden Werten (ca.) 159cm und 160cm.

Welche mittlere Größe wird mit 95% Wahrscheinlichkeit nicht überschritten, wenn die Person 45kg wiegt?

mu_at_45 %>% 
  summarise(q_95 = quantile(mu_at_45, prob = .95))

8.7 Prior-Prädiktiv-Verteilung

🏎🏎️ VERTIEFUNG (nicht prüfungsrelevant ) 🏎🏎

8.7.1 Moment

🤔 Moment. Dieser Prior, \(\beta_1\) in m43 erachtet positive und negative Zusammenhang als gleich wahrscheinlich?!

Sind wir wirklich indifferent, ob der Zusammenhang von Gewicht und Größe positiv oder negativ ist? Nein, sind wir nicht.

8.7.2 Priori-Prädiktiv-Verteilung für m43

Was denkt wir bzw. unser Golem apriori über den Zusammenhang von Größe und Gewicht? Um diese Frage zu beantworten ziehen wir Stichproben aus den Priori-Verteilungen des Modells, also für \(\beta_0\), \(\beta_1\) und \(\sigma\).

m43_prior_pred <-
    stan_glm(height ~ weight_c, 
             prior = normal(0, 10),
             prior_intercept = normal(178, 20),  # mu
             prior_aux = exponential(0.1),  # sigma
             refresh = FALSE, 
             prior_PD = TRUE,  # Schalter für Prior-Pred-Verteilung
             data = d3)


m43_prior_pred_draws <- 
  m43_prior_pred %>% 
  as_tibble() %>% 
  rename(a = `(Intercept)`,
         b = weight_c) %>% 
  slice_sample(n = 50)

a	b	sigma
170.3	15.0	7.0
140.2	6.0	6.5
155.5	−4.2	6.5
159.6	13.3	20.9
178.1	2.5	4.0

Jede Zeile definiert eine Regressionsgerade.

8.7.3 Prior-Prädiktiv-Simulation für m43 mit stan_glm()

m43_prior_pred <-
    stan_glm(height ~ weight_c, 
             prior = normal(0, 10),  # beta
             prior_intercept = normal(178, 20),  # alpha
             prior_aux = exponential(0.1),  # sigma
             refresh = FALSE, 
             prior_PD = TRUE,  # DIESER Schalter macht's
             data = d3)

m43_prior_pred_draws <- 
  m43_prior_pred %>% 
  as_tibble() %>% 
  rename(a = `(Intercept)`,
         b = weight_c) %>% 
  slice_sample(n = 50)

8.7.4 Visualisieren der Prior-Prädiktiv-Verteilung

d3 %>% ggplot() +
  geom_point(aes(x = weight_c, y = height)) + 
  geom_abline(data = m43_prior_pred_draws,
aes(intercept = a, slope = b), color = "skyblue", size = 0.2) +
  scale_y_continuous(limits = c(0, 500)) +
  geom_hline(yintercept = 272, size = .5) +
  geom_hline(yintercept = 0, linetype = "dashed")

🤯 Einige dieser Regressionsgeraden sind unsinnig!

d3 %>% ggplot() +
  geom_point(aes(x = weight_c, y = height)) + 
  geom_abline(data = m43_prior_pred_draws,
aes(intercept = a, slope = b), color = "skyblue", size = 0.2) +
  scale_y_continuous(limits = c(0, 500)) +
  geom_hline(yintercept = 272, size = .5) +
  geom_hline(yintercept = 0, linetype = "dashed")

Die durchgezogene horizontale Linie gibt die Größe des größten Menschens, Robert Pershing Wadlow, an.

8.7.5 Ein positiver Wert für \(\beta_1\) ist plausibler

8.7.5.1 Oh no

Eine Normalverteilung mit viel Streuung:

👎 \(\beta=-20\) wäre mit diesem Prior gut möglich: Pro kg Gewicht sind Menschen im Schnitt 20cm kleiner, laut dem Modell. Quatsch.

8.7.5.2 Oh yes

Wir bräuchten eher so eine Verteilung, mit mehr Masse auf der positiven Seite (x>0):

👍 Vermutlich besser: Ein Großteil der Wahrscheinlichkeitsmasse ist \(X>0\). Allerdings gibt’s keine Gewähr, dass unser Prior “richtig” ist.

8.7.6 Priori-Prädiktiv-Simulation, 2. Versuch

m43a_prior_pred <-
    stan_glm(
      height ~ weight_c, 
      prior = normal(2, 2),  # Regressionsgewicht
      prior_intercept = normal(178, 20),  # mu
      prior_aux = exponential(0.1),  # sigma
      refresh = FALSE, 
      # Schalter für Prior-Pred-Verteilung:
      prior_PD = TRUE, 
      data = d3)


m43a_prior_pred_draws <- 
  m43a_prior_pred %>% 
  as_tibble() %>% 
  # Spaltennamen kürzen: 
  rename(a = `(Intercept)`) %>%  
  rename(b = weight_c,
         s = sigma)

a	b	s
204.2	4.1	14.7
175.1	−0.6	2.3
149.4	0.3	1.7
172.3	0.3	31.2
173.7	1.0	12.0

Das Argument prior_PD = TRUE sorgt dafür, dass keine Posteriori-Verteilung, sondern eine Prior-Prädiktiv-Verteilung berechnet wird.

8.7.7 Visualisieren der Prior-Prädiktiv-Verteilung, m43a

Unsere Priori-Werte scheinen einigermaßen vernünftige Vorhersagen zu tätigen. Allerdings erwartet unser Golem einige Riesen.

d3 %>% 
  ggplot(aes(x = weight_c, y = height)) +
  geom_point() +
  geom_abline(data = {m43a_prior_pred_draws %>% slice_head(n=50)},
              aes(slope = b,
                  intercept = a),
              color = "skyblue",
              size = .2,
              alpha = .7) +
  geom_hline(yintercept = 272, size = .5) +
  geom_hline(yintercept = 0, linetype = "dashed")+
  scale_y_continuous(limits = c(0, 500)) 

Die durchgezogene horizontale Linie gibt die Größe des größten Menschens, Robert Pershing Wadlow, an.

8.7.8 Moment, kann hier jeder machen, was er will?

Es doch den einen, richtigen, objektiven Priori-Wert geben?!

Kann denn jeder hier machen, was er will?! Wo kommen wir da hin?!

This is a mistake. There is no more a uniquely correct prior than there is a uniquely correct likelihood. Statistical models are machines for inference. Many machines will work, but some work better than others. Priors can be wrong, but only in the same sense that a kind of hammer can be wrong for building a table.

McElreath (2020), p. 96.

8.7.9 Hier ist unser Modell, m43a

\[\begin{align} \text{height}_i &\sim \operatorname{Normal}(\mu_i, \sigma) \\ \mu_i &= \alpha + \beta \cdot \text{weight}_i\\ \alpha &\sim \operatorname{Normal}(178, 20)\\ \beta &\sim \operatorname{Normal}(5,3)\\ \sigma &\sim \operatorname{Exp}(0.1) \end{align}\]

# Posteriori-Vert. berechnen:
m43a <-
  stan_glm(
    height ~ weight_c,  # Regressionsformel
    prior = normal(5, 3),  # Regressionsgewicht (beta 1)
    prior_intercept = normal(178, 20),  # mu
    prior_aux = exponential(0.1),  # sigma
    refresh = 0,  # zeig mir keine Details
    seed = 42,  # Zufallszahlen festlegen
    data = d3)

8.7.10 Eine Zusammenfassung der Posteriori-Verteilung für m43a

m43a %>% 
  parameters()

Parameter	Median	95% CI	pd	Rhat	ESS	Prior
(Intercept)	154.65	(154.14, 155.19)	100%	0.999	3214.00	Normal (178 +- 20)
weight_c	0.91	(0.82, 0.99)	100%	1.001	4134.00	Normal (5 +- 3)

Unser Modell m43a schätzt die typische Körpergröße einer !Kung-Person mittleren Gewichts (weight_c = 0) auf knapp 155 cm, und ist sich dieses Werts ziemlich sicher. Pro Kilogramm kommt (laut unserem Modell) ein knapper Zentimeter hinzu, typischerweise; auch hier ist sich das Modell ziemlich sicher, da dass zugehörige 95%-CI keine 20 Zentimenter umfasst.

8.8 Die PPV befragen

🏎🏎️ VERTIEFUNG (nicht prüfungsrelevant ) 🏎🏎

Die Posterior-Prädiktiv-Verteilung (PPV) gibt uns die Möglichkeit, nach der Wahrscheinlichkeit tatsächlicher Körpergrößen zu fragen - und nicht nur nach mittleren Körpergrößen anhand der Post-Verteilung.

Wichtig

Die Post-Verteilung macht nur Aussagen zur mittleren Körpergröße, denn das ist was wir modellieren wollten. Möchten wir Aussagen zur Wahrscheinlichkeit tatsächlicher Größen treffen, brauchen wir die PPV. Allgemeiner gesagt: Die PPV macht Vorhersagen auf Basis eines Modells. Für jede Vorhersage gibt es eine Verteilung, die wir zu einem Punktschätzer (z.B. Median) und einem Schätzbereich (z.B. 89%-HDI) zusammenfassen können.

An dieser Stelle sollten wir uns vor Augen führen, dass die PPV mehr Ungewissheit beinhaltet, denn sie berücksichtigt derer zweier Arten:

1. Ungewissheit bzgl. der Modellparameter (Steigung und Achsenabschnitt der Regressionsgeraden)
2. Ungewissheit der Vorhersagen (das Modell macht keine perfekten Vorhersagen)

Die Post-Verteilung berücksichtigt nur die Ungewissheit in den Modellparametern, macht also nur Aussagen zur Regressionsgeraden.

Die PPV macht Aussagen für konkrete Beobachtungen. Der Unterschied ist in Abbildung 8.16 dargestellt; die Funktionen stammen übrigens aus {easystats}.

estimate_prediction(m43a) %>% plot()
estimate_relation(m43a) %>% plot()

(a) PPV mit viel Ungewissheit

(b) Post-Verteilung mit wenig(er) Ungewissheit

Abbildung 8.16: PPV vs. Post-Verteilung

8.8.1 Perzentil-Intervalle für bestimmte Prädiktor-Werte

Wie groß ist ein !Kung-Mann mit mittlerem Gewicht?

set.seed(42)
estimate_prediction(m43a, data = tibble(weight_c = 0), seed = 42)

Unser Modell, ma43a schätzt ca. 145cm bis 165cm.

Wir können uns auch eine Sequenz an Prädiktorwerten, die uns interessieren, erstellen, s. weight_df:

weight_df <- tibble(weight_c = seq(-20,20, by = 5))

Für diese Werte lassen wir uns dann die Perzentil-Intervalle (PI) ausgeben:

mus <- 
  estimate_prediction(m43a, data = weight_df) 

head(mus)

Um die Perzentilintervalle zu erstellen, wird von estimate_prediction() für jeden Prädiktorwert eine PPV erstellt und (in der Voreinstellung) das 5%- sowie 95%-Quantil dafür berechnet. Sie können die Voreinstellung ändern mittels des Arguments ci; um ein 89%-PI zu berechnen, würde man z.B. schreiben ci = .89.

Um Reproduzierbarkeit sicherzustellen, haben wir mit set.seed(42) die Zufallszahlen fixiert.

Hoppla! Das ist ja viel ungenauer, als die Angaben der Post-Verteilung oben. Ja, denn die Post-Verteilung hat die Ungewissheit zum Mittelwert ausgedrückt; die PPV gibt die Ungewissheit tatsächlicher beobachtbarer Körpergrößen aus, nicht nur die Ungewissheit zum Mittelwert.

Berechnen wir die PPV für die bestehenden Beobachtungen aus m43a:

ppv_m43a <- estimate_prediction(
  m43a,
  data = weight_df)

mus 

8.8.2 Perzentilintervalle für verschiedenen Prädiktorwerte visualisiert

Abbildung 8.17 visualisiert die Ungewissheit von Vorhersagen laut der PPV. Die Ungewissheit in Abbildung 8.17 ist die Antwort auf die Frage: “Wie sicher sind wir uns, zur Größe einer !Kung-Person, gegeben dass die z.B. 10 kg mehr als der Durchschnitt wiegt?” Eine Vorhersage bezeichnet man auch als “bedingte Verteilung”, da man den Wert einer Verteilung voraussagt, gegeben einer Bedingung, z.B. weight_c = 10.

Abbildung 8.17: Visualisierung der Ungewissheit der Vorhersagen laut PPV: Die Vorhersage von Beobachtungen beinhaltet mehr Ungewissheit als die Vorhersage von zu erwartenden Werten, daher sind die Ungewissheitsintervalle der PPV größer als die der Post-Verteilung.

Die vertikalen Balken geben die 95%-KI wieder, die wir jeweils zu erwarten haben.

Noch eine andere Visualisierung, s. Abbildung 8.18; je dicker die “Katzenaugen”, desto mehr Stichproben (samples) liegen vor an der Stelle, und umso genauer ist die Schätzung.

Abbildung 8.18: Die PPV für bestimmte Gewichtswerte w, visualisiert mit Katzenaugen-Diagrammen

Also: Je dicker die Violine, desto wahrscheinlicher \(\mu | w_i\).

8.8.3 Die PPV über alle Beobachtungen visualisiert

Gerade eben haben wir bedingte PPVen angeschaut: Also eine PPV für einen bestimmten Prädiktorwert, z.B. bei einer Person mittleren Gewichts. Wir können auch den Mittelwert über alle bedingten PPV anschauen, sozusagen die “Master-PPV” oder “unbedingte PPV” oder schlicht PPV. Vergleichen wir die echten Werte für height, \(y\), mit den von der PPV simulierten Werten für height, \(y_{rep}\), s. Abbildung 8.19.

check_predictions(m43a)  # aus easystatss

Abbildung 8.19: Vergleich der Vorhersagen für y (leichte, blaue Linien) mit der beobachteten Verteilung von y

?check_predictions zeigt Hilfe für diese Funktion. Die Funktion zeigt die Vorhersagen für die AV laut der Posteriori-Verteilung.

Die zwei Gipfel hat unser Modell nicht mitgekriegt, ansonsten decken sich die Vorhersagen der PPV gut mit den echten Daten.

8.8.4 Fragen an die Master-PPV

• Wie groß sind die !Kung im Schnitt?
• Welche Größe wird von 90% der Personen nicht überschritten?
• Wie groß sind die 10% kleinsten?

ppv_m43a <- posterior_predict(
  m43a,
  newdata = weight_df,
  draws = 100) %>% 
  as_tibble() %>% 
  pivot_longer(
    cols = everything(),
    names_to = "weight_condition",
    values_to = "height")
head(ppv_m43a)

ppv_m43a <-
  ppv_m43a <- posterior_predict(
  m43a,
  newdata = weight_df,
  draws = 100) %>% 
  as_tibble() %>% 
  pivot_longer(
    cols = everything(),
    names_to = "weight_condition",
    values_to = "height")

head(ppv_m43a)

ppv_m43a %>% 
  summarise(
    q_10 = quantile(height, prob = .1),
    height_mean = mean(height),
    q_50 = quantile(height, prob = .5),
    q_90 = quantile(height, prob = .9)
  )

Was ist der 50% Bereich der Körpergröße?

ppv_m43a %>% 
  eti(ci = .5)

8.9 Fazit

8.9.1 Ausstieg

Beispiel 8.3 (Fassen Sie das Wesentliche zusammen!) Schreiben Sie 5-10 Sätze zum Wesentlichen Stoff dieses Kapitels und reichen Sie bei der von Lehrkraft vorgegebenen Stelle ein! \(\square\)

8.9.2 Vertiefung

McElreath (2020) bietet eine tiefere Darstellung von linearen Modellen auf Basis der Bayes-Statistik, insbesondere Kapitel 4 daraus vertieft die Themen dieses Kapitels. Kurz (2021) greift die R-Inhalte von McElreath (2020) auf und setzt sie mit anderen R-Methoden um; ein interessanter Blickwinkel, wenn man tiefer in die R-Umsetzung einsteigen möchte. Gelman et al. (2021) bieten ebenfalls viele erhellende Einblicke in das Thema Regressionsanalyse, sowohl aus einem frequentistischen als auch aus einer Bayes-Perspektive.

8.10 Aufgaben

8.10.1 Papier-und-Bleistift-Aufgaben

1. Bayes-Ziel1
2. Bayesmod-bestimmen01
3. mtcars-post2a
4. Bed-Post-Wskt1
5. mtcars-post3a
6. regression1a
7. regression1b
8. Regression2
9. Priorwahl1
10. Bayesmod-bestimmen02
11. Aussagen-einfache-Regr
12. Likelihood-identifizieren
13. Priorwahl2
14. Likelihood2
15. interpret-koeff
16. bed-post-wskt1

8.10.2 Computer-Aufgaben

1. Post-befragen1
2. penguins-stan-01

8.11 —

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

Kurz, S. (2021). Statistical Rethinking with Brms, Ggplot2, and the Tidyverse: Second Edition. https://bookdown.org/content/4857/

Makowski, D., Ben-Shachar, M. S., Chen, S. H. A., & Lüdecke, D. (2019). Indices of Effect Existence and Significance in the Bayesian Framework. Frontiers in Psychology, 10. https://www.frontiersin.org/article/10.3389/fpsyg.2019.02767

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

---

1. Da es viele Funktionen sind, bietet es sich an mit Strg-F auf der Webseite nach Ihrem Lieblingsbefehl zu suchen.

2. Datenquelle, McElreath (2020).

3. Wer ist hier für die Namensgebung zuständig? Besoffen oder was?

4. oder zumindest besser sichergestellt

5. In manchen Lehrbüchern wird \(\beta_0\) auch als \(\alpha\) bezeichnet.

6. 1e6

7. Im Standard beschert uns stan_glm() 4000 Stichproben.

8. \(\mu\)