penguins-stan-01

bayes
regression
string
qm2
qm2-pruefung
Published

July 12, 2023

Aufgabe

Wir untersuchen Einflussfaktoren bzw. Prädiktoren auf das Körpergewicht von Pinguinen. In dieser Aufgabe untersuchen wir in dem Zusammenhang den Zusammenhang von Schnabellänge (als UV) und Körpergewicht (als AV).

Wie groß ist der statistische Einfluss der UV auf die AV?

  1. Berechnen Sie den Punktschätzer des Effekts!
  2. Wie viele Parameter hat das Modell?
  3. Geben Sie die Breite eines 90%-HDI an (zum Effekt)!
  4. Wie groß ist die Wahrscheinlichkeit, dass der Effekt vorhanden ist (also größer als Null ist), die “Effektwahrscheinlichkeit”?
  5. Wie groß ist das 95%-HDI, wenn Sie nur die Spezies Adelie untersuchen?
  6. Geben Sie die Prioris an für m1 für die Regressionskoeffizienten!

Hinweise:

  • Nutzen Sie den Datensatz zu den Palmer Penguins.
  • Verwenden Sie Methoden der Bayes-Statistik und die Software Stan.
  • Fixieren Sie die Zufallszahlen auf den Startwert 42!
  • Sie können den Datensatz z.B. hier beziehen oder über das R-Paket palmerpenguins.
  • Geben Sie keine Prozentzahlen, sondern stets Anteile an.
  • Beachten Sie die übrigen Hinweise.











Lösung

Zentrieren ist eigentlich immer nützlich, aber hier streng genommen nicht unbedingt nötig. Der Hauptgrund ist, dass Stan für uns den Prior für den Intercept festlegt, und zwar auf Basis der Daten, wir uns also nicht um die komische Frage zu kümmern brauchen, welchen Prior wir für den unzentrierten Achsenabschnitt vergeben wollten: Wie schwer sind Pinguins der Schnabellänge Null? Mit zentrierten Prädiktoren ist die Frage nach dem Prior viel einfacher zu beantworten: Wie schwer ist ein Pinguin mit mittelgroßem Schnabel?

Setup:

library(tidyverse)
library(easystats)
library(rstanarm)

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

Es wird in dieser Aufgabe vorausgesetzt, dass Sie den Datensatz selbständig importieren können. Tipp: Kurzes Googeln hilft ggf., den Datensatz zu finden.

Alternativ könnten Sie den Datensatz als CSV-Datei importieren:

d_path <- "https://vincentarelbundock.github.io/Rdatasets/csv/palmerpenguins/penguins.csv"
penguins <- data_read(d_path)  # oder z.B. mit read_csv

Ein Blick in die Daten zur Kontrolle, ob das Importieren richtig funktioniert hat:

glimpse(penguins)
Rows: 344
Columns: 8
$ species           <fct> Adelie, Adelie, Adelie, Adelie, Adelie, Adelie, Adel…
$ island            <fct> Torgersen, Torgersen, Torgersen, Torgersen, Torgerse…
$ bill_length_mm    <dbl> 39.1, 39.5, 40.3, NA, 36.7, 39.3, 38.9, 39.2, 34.1, …
$ bill_depth_mm     <dbl> 18.7, 17.4, 18.0, NA, 19.3, 20.6, 17.8, 19.6, 18.1, …
$ flipper_length_mm <int> 181, 186, 195, NA, 193, 190, 181, 195, 193, 190, 186…
$ body_mass_g       <int> 3750, 3800, 3250, NA, 3450, 3650, 3625, 4675, 3475, …
$ sex               <fct> male, female, female, NA, female, male, female, male…
$ year              <int> 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007, 2007…

Vertrauen ist gut, aber - was Golems betrifft - ist Kontrolle eindeutig besser ;-)

  1. Punktschätzer
m1 <- stan_glm(body_mass_g ~  bill_length_mm,  # Regressionsgleichung
               data = penguins, #  Daten
               seed = 42,  # Reproduzierbarkeit
               refresh = 0)  # nicht so viel Output
parameters(m1, ci_method = "hdi", ci = .9)
Parameter      | Median |            90% CI |     pd |  Rhat |     ESS |                       Prior
----------------------------------------------------------------------------------------------------
(Intercept)    | 361.85 | [-108.79, 834.20] | 89.55% | 0.999 | 3920.00 | Normal (4201.75 +- 2004.89)
bill_length_mm |  87.45 | [  77.17,  98.51] |   100% | 0.999 | 3931.00 |     Normal (0.00 +- 367.22)

Uncertainty intervals (highest-density) and p-values (two-tailed)
  computed using a MCMC distribution approximation.
  1. Anzahl Parameter

Das Modell hat 3 Paramter:

  • \(\beta_0\) (oder \(\alpha\))
  • \(\beta_01\)
  • \(\sigma\)
  1. Breite des Intervalls

Dazu liest man die Intervallgrenzen (90% CI) in der richtigen Zeile ab (Tabelle parameters):

 97.70  - 76.24
[1] 21.46

Einheit: mm

  1. Effektwahrscheinlichkeit
m1_post <-
  m1 %>% 
  as_tibble()

m1_post %>% 
  count(bill_length_mm > 0)
# A tibble: 1 × 2
  `bill_length_mm > 0`     n
  <lgl>                <int>
1 TRUE                  4000

Also: 100% oder 1 (4000 von 4000 Stichproben finden dieses Ergebnis in unserem Modell).

Man kann diesen Wert aus der Tabelle oben (Ausgabe von parameters()) einfach in der Spalte pd ablesen. pd steht für probability of direction, s. Details hier.

Oder so, ist auch einfach:

pd_m1 <- p_direction(m1) # aus Paket easystats
pd_m1
Probability of Direction

Parameter      |     pd
-----------------------
(Intercept)    | 89.55%
bill_length_mm |   100%

Und plotten ist meist hilfreich: plot(pd_m1).

Man kann sich auch ein “Dashboard” mit allen Ergebnissen des Modells ausgeben lassen:

model_dashboard(m1)
  1. Nur Adelie:

Welche Spezies gibt es im Datensatz?

penguins %>% 
  count(species)
# A tibble: 3 × 2
  species       n
  <fct>     <int>
1 Adelie      152
2 Chinstrap    68
3 Gentoo      124

Filtern:

penguins_adelie <-
  penguins %>% 
  filter(species == "Adelie")

Modell berechnen:

m2 <- stan_glm(body_mass_g ~  bill_length_mm,  # Regressionsgleichung
               data = penguins_adelie, #  Daten
               seed = 42,  # Repro.
               refresh = 0)  # nicht so viel Output

Das Modell ist - bis auf die Daten - identisch zu m1.

parameters(m2)
Parameter      | Median |            95% CI |     pd |  Rhat |     ESS |                       Prior
----------------------------------------------------------------------------------------------------
(Intercept)    |  18.20 | [-881.48, 930.13] | 51.75% | 1.000 | 3988.00 | Normal (3700.66 +- 1146.42)
bill_length_mm |  94.71 | [  71.21, 118.31] |   100% | 1.000 | 3962.00 |     Normal (0.00 +- 430.43)

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a MCMC distribution approximation.
hdi(m2, parameters = "bill_length_mm")
Highest Density Interval

Parameter      |         95% HDI
--------------------------------
bill_length_mm | [70.88, 117.70]

S. auch Tabelle oben.

118.09 - 71.86
[1] 46.23
  1. Prioris
describe_prior(m1, component = "auxiliary")
       Parameter Prior_Distribution Prior_Location Prior_Scale
1    (Intercept)             normal       4201.754   2004.8863
2 bill_length_mm             normal          0.000    367.2233

Steht auch in der Tabelle, die von parameters ausgegeben wird.

Categories:

  • bayes
  • regression
  • string