9  Einfache lineare Modelle

9.1 Lernsteuerung

9.1.1 Position im Modulverlauf

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

9.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

9.1.3 Begleitliteratur

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

9.1.4 Vorbereitung im Eigenstudium

• Statistik1, Kap. “Geradenmodelle 1”

9.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.¹

9.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"

kung <- read.csv(Kung_path)

kung_erwachsen <- kung %>% 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)

9.1.7 Einstieg

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

Beispiel 9.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\)

9.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.

9.2 Post-Verteilung der Regression

9.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, UV) und Größe (Y, AV), Abbildung 9.1.

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

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

Mit ggplot2

Mit ggpubr

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

9.2.2 Statistiken zum !Kung-Datensatz

Die Daten können Sie hier herunterladen.

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

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

kung_erwachsen <- kung %>% filter(age > 18)

describe_distribution(kung_erwachsen)

Tabelle 9.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 9.1 abzulesen ist, beträgt das mittlere Körpergewicht (weight) liegt ca. 45kg (sd 7 kg).

9.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)

9.2.3.1 Gibt es fehlende Werte?

Nein, s. Abb. Abbildung 9.2.

kung_erwachsen %>% plot_missing()

Abbildung 9.2: Fehlende Werte - fehlen.

9.2.3.2 Verteilung der numerischen Variablen

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

kung_erwachsen %>% plot_histogram()

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

9.2.3.3 Verteilung der kategorialen Variablen

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

kung_erwachsen %>% plot_bar()

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

9.2.3.4 Korrelationen

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

kung_erwachsen %>% plot_correlation()

Abbildung 9.5: Korrelationsmatrix

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

9.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 - \bar{x}\).

kung_zentriert <-
  kung_erwachsen %>% 
  mutate(weight_c = weight - mean(weight))

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

kung_zentriert <- 
  kung_erwachsen %>% 
  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 9.6.

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

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

9.3 Modell m_kung_gewicht_c: 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.

9.3.1 Modelldefinition von m_kung_gewicht_c

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 exponentialverteilt angenommen. Der Likelihood gibt an, wie wahrscheinlich ein Wert height ist, gegeben \(\mu\) und \(\sigma\). Theorem 9.1 stellt die Modelldefinition dar. ?fig-kruschke_regr_one_predctor zeigt, wie die drei Parameter zusammen die Likelihood definieren.

Theorem 9.1 (Modelldefinition m_kung_gewicht_c) \[\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*}\quad \square\]

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. 🤷

9.3.2 Likelihood, m_kung_gewicht_c

\[ \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 m_kung_gewicht_c 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\))”.

9.3.3 Regressionsformel, m_kung_gewicht_c

\[ \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).

9.3.4 Priori-Werte des Modells m_kung_gewicht_c

\[\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.

9.3.5 Prioris plus Daten gleich Post

![Parameter ]

9.4 Die Post-Verteilung befragen

📺 Post-Verteilung auslesen 1

📺 Post-Verteilung auslesen 2

9.4.1 m_kung_starkes_beta1

Sagen wir, auf Basis gut geprüfter Evidenz haben wir folgendes Modell festgelegt: height ~ weight_c, s. Gleichung 9.1. Dabei haben wir folgende Prioris gewählt:

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

Wir nennen das Modell m_kung_starkes_beta1³, s. Listing 9.1.

Listing 9.1: Modelldefinition von m_kung_starkes_beta1 in R

m_kung_starkes_beta1 <-
  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
    data = kung_zentriert)

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.

9.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 9.2.

Tabelle 9.2: Parameter von m_kung_starkes_beta1

parameters(m_kung_starkes_beta1)

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

Definition 9.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 9.7.

plot(p_direction(m_kung_starkes_beta1))

Abbildung 9.7: Diagramm zur Probability of Direction, Modell m_kung_starkes_beta1

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.

9.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 m_kung_starkes_beta1:

m_kung_starkes_beta1 %>% 
  as_tibble() %>% 
  head()

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

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

estimate_expectation(m_kung_starkes_beta1, by = "weight_c") |> plot()  # aus {easystats}

Mit ggplot

Mit easystats

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

(a) Erwartete Werte des Modell m_kung_starkes_beta1, sprich, die Regressionsgerade

Abbildung 9.8

9.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.

9.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 9.2.

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

m_kung_starkes_beta1_post <- 
  m_kung_starkes_beta1 %>% 
  as_tibble()

m_kung_starkes_beta1_post %>% 
  head()

Wie wir gesehen haben, nutzen wir diese Tabelle der Post-Verteilung immer wieder. Speichern wir uns sie also als ein Objekt ab, m_kung_starkes_beta1_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 9.9.

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

Abbildung 9.9: Die Postverteilung für den Parameter Gewicht (zentriert); die plausiblen Werte liegen zwischen 0.8 und 1.0 cm pro Kilogramm Gewicht

Abbildung 9.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:

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

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

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

Abbildung 9.10: Die Post-Verteilung für den Parameter sigma, m_kung_starkes_beta1; die plausiblen Werte liegen zwischen 4.7 cm und 5.5 cm.

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

m_kung_starkes_beta1_hdi <- hdi(m_kung_starkes_beta1_post)  # analog mit eti(m_kung_starkes_beta1)

plot(m_kung_starkes_beta1_hdi)

Abbildung 9.11: Die Parameter Gewicht (zentriert) und sigma des Modells m_kung_starkes_beta1

Ergänzt man bei plot() noch show_intercept = TRUE wird auch der Achsenabschnitt angezeigt. Den Parameter der Vorhersage-Genauigkeit, \(\sigma\) bekommt man mit get_sigma:

get_sigma(m_kung_starkes_beta1)
## [1] 5.1
## attr(,"class")
## [1] "insight_aux" "numeric"

9.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.

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

Abbildung 9.12 stellt die Ungewissheit der Post-Verteilung dar, in dem einige Stichproben aus der Post-Verteilung visualisiert werden.

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

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

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

10

Die ersten 10 Stichproben:

100

Die ersten 100 Stichproben:

1000

Die ersten 1e3 Stichproben:

1000000

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

Abbildung 9.12

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

estimate_expectation(m_kung_starkes_beta1, by = "weight_c") %>% plot()

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

9.4.7 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:

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

50%-PI:

m_kung_starkes_beta1 %>% 
  eti(ci = .5)

9.4.8 Fragen zu Wahrscheinlichkeitsmassen des Achsenabschnitts

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

m_kung_starkes_beta1_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?

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

Die Wahrscheinlichkeit beträgt 0.29.

9.4.9 Typischer Bayes-Nutzer in Aktion

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

Quelle

9.5 Post-Verteilung bedingt auf einen Prädiktorwert

9.5.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 m_kung_post, s. Abbildung 9.14.

Abbildung 9.14: 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.

9.5.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 <-
  m_kung_gewicht_c_post %>% 
  mutate(mu_at_45 = `(Intercept)`)

Und plotten diese, s. Abbildung 9.15.

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

Abbildung 9.15: 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 9.3.

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

head(mu_at_50)

Tabelle 9.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 9.16.

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

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

9.5.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(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))

9.6 Fazit

9.6.1 Ausstieg

Beispiel 9.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\)

9.6.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 Tidyverse-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 einer frequentistischen als auch aus einer Bayes-Perspektive.

9.7 Aufgaben

9.7.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

9.7.2 Computer-Aufgaben

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

9.8 —

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\)