penguins-interact

Aufgabe

Eine Forscherin untersucht, ob das Geschlecht eines Pinguins den Einfluss der Flossenlänge (Flipper, mm) auf das Körpergewicht (g) des Tieres moderiert.

Auf Basis der Daten: Liegt ein (substanzieller) Interaktionseffekt vor?

Hinweise:

• Nutzen Sie die folgende Analyse als Grundlage Ihrer Antworten.
• Beachten Sie die Hinweise des Datenwerks.
• Unter “substanziell” sei ein Effekt von mind. 100 g verstanden.

Setup:

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

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

Dafür ist folgende Analyse gegeben.

Setup

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

Modell und Hypothese

Die Forschungsfrage kann man wie folgt als Hypothese formalisieren:

\[\beta_{m} \ne 0\]

“Der Regressionskoeffizient der Moderation (\(m\), d.h. Interaktion) ist ungleich Null.”

Testet man nicht eine exakte, sondern einen Mindestwert (ROPE), so kann man die Hypothese so formulieren:

\[\beta_{m} > 100\]

Die Prioris übernehmen wir vom Stan-Golem.🤖

🤖 Beep, beep!

👩‍🏫 An die Arbeit, Stan-Golem!

Wir könnten den Datensatz auch 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…

Wir entfernen noch alle fehlenden Werte:

penguins_nona <- 
  penguins |> 
  filter(sex == "female" | sex == "male")

penguins_nona$sex |> unique()

[1] male   female
Levels: female male

Zur besseren Interpretierbarkeit standardisieren wir die (metrische) UV:

penguins_nona_z <- 
  penguins_nona |> 
  standardise(select = "flipper_length_mm",
              append = TRUE)

m_interaction <- stan_glm(body_mass_g ~  sex + flipper_length_mm_z + sex:flipper_length_mm_z,  # Regressionsgleichung
               data = penguins_nona_z, #  Daten
               seed = 42,  # Reproduzierbarkeit
               refresh = 0)  # nicht so viel Output

m_interaction_params <- parameters(m_interaction, ci_method = "hdi", ci = .9)
m_interaction_params

Parameter	Median	CI	CI_low	CI_high	pd	Rhat	ESS	Prior_Distribution	Prior_Location	Prior_Scale
(Intercept)	4032.207724	0.9	3986.27445	4081.67307	1.00000	1.001110	3591.836	normal	4207.057	2013.040
sexmale	347.079159	0.9	282.11779	415.10070	1.00000	1.001413	3910.908	normal	0.000	4020.192
flipper_length_mm_z	659.450595	0.9	608.70219	713.30490	1.00000	1.000452	2446.063	normal	0.000	2013.040
sexmale:flipper_length_mm_z	-3.583776	0.9	-70.50816	65.73824	0.53025	1.000672	2640.031	normal	0.000	2695.055

Lösung

Punktschätzer

Der Punktschätzer ist in der Spalte Median in der Tabelle parameters zu finden. Sein Wert ist:

[1] -3.583776

Hier ist die Post-Verteilung des Effekts:

m_interaction_params |> plot()

Alternative Visualisierung:

hdi(m_interaction, ci = .9) |> plot()

Breite des Intervalls

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

Obere Grenze: 65.738242.

Untere Grenze: -70.5081637.

Differenz = Obere_Grenze - Untere_Grenze:

[1] 136.2464

Einheit: mm

Effektwahrscheinlichkeit

Man erkennt schon im Diagramm zum Konfidenzintervall, dass 100% des Intervalls positiv ist. Daher ist die Effektwahrscheinlichkeit auch positiv.

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_m_interaction_params <- p_direction(m_interaction) # aus Paket easystats
pd_m_interaction_params

	Parameter	pd	Effects	Component
1	(Intercept)	1.00000	fixed	conditional
3	sexmale	1.00000	fixed	conditional
2	flipper_length_mm_z	1.00000	fixed	conditional
4	sexmale:flipper_length_mm_z	0.53025	fixed	conditional

Und plotten ist meist hilfreich: plot(pd_m_interaction_params).

plot(pd_m_interaction_params)

Substanzielle Effektwahrscheinlichkeit

Die Frage ist nichts anderes als nach ROPE zu fragen.

rope_m_interact <- rope(m_interaction, range = c(-100,+100))
rope_m_interact

	Parameter	CI	ROPE_low	ROPE_high	ROPE_Percentage	Effects	Component
1	(Intercept)	0.95	-100	100	0	fixed	conditional
3	sexmale	0.95	-100	100	0	fixed	conditional
2	flipper_length_mm_z	0.95	-100	100	0	fixed	conditional
4	sexmale:flipper_length_mm_z	0.95	-100	100	1	fixed	conditional

plot(rope_m_interact)

Das 95%-Post-Intervall ist komplett innerhalb des ROPE.

Wir können die ROPE-Hypothese daher bestätigen.