pupil-size2

Aufgabe

Pupillendaten sind ein verbreiteter Analysegegenstand in Bereichen wie Psychologie, Marktforschung und Marketing.

Betrachten wir dazu ein R-Paket (zum Vorverbarbeitung, preprocessing) und einen Datensatz der Uni Münster.

library(tidyverse)
library(PupilPre)  # installieren, einmalig, nicht vergessen
data("Pupildat")
d <-
  Pupildat %>% 
  select(size = RIGHT_PUPIL_SIZE,
         time = TIMESTAMP) %>% 
  mutate(size = size / 100) # in millimeter

Ein Forschungsteam untersucht den Datensatz und möchte ein Modell zur Schätzung der Pupillengröße aufstellen.

Mit dem R-Paket eaystats kann man sich bequem typische Statistiken ausgeben lassen, was das Forschungsteam auch macht.

library(easystats)
d %>% 
  describe_distribution()

Variable	Mean	SD	IQR	Min	Max	Skewness	Kurtosis	n	n_Missing
size	1.000562e+01	5.107151e+00	3.88	1.04	25.01	1.2474408	0.3199426	45343	1607
time	2.990837e+06	9.874538e+05	1947161.00	1443974.00	4062110.00	-0.4111739	-1.6982498	46950	0

Ds Forschungsteam verzichtet hier auf eine Aufbereitung der Daten (was eigentlich nötig wäre, aber nicht Gegenstand dieser Übung ist). Stattdessen konzentrieren wir uns auf die Posteriori-Verteilung zur Pupillengröße.

Das Forschungsteam ist interessiert an einem Modell zur Schätzung der (Verteilung der) Pupillengröße; die Posteriori-Verteilung bildet das ab.

Zuerst definiert das Forschungsteam ein Modell:

\[\begin{aligned} s_i &\sim \mathcal{N}(\mu, \sigma)\qquad \text{| s wie size }\\ \mu &\sim \mathcal{N}(10, 5)\\ \sigma &\sim \mathcal{E}(.2) \end{aligned}\]

Für das Modell wird folgende Begründung vom Forschungsteam gegeben:

\(s_i\): Pupillengrößen sind normalverteilt, da viele Gene additiv auf die Größe hin zusammenwirken

\(\mu\): Da wir nicht viel wissen über die mittlere Pupillengröße, entscheiden wir uns für Normalverteilung für diesen Parameter, da dies keine weiteren Annahmen (außer dass Mittelwert und Streuung endlich sind) hinzufügt. Ein Modell mit wenig Annahmen nennt man “sparsam” oder konservativ. Es ist wünschenswert, dass Modelle mit so wenig wie möglich Annahmen auskommt (aber so vielen wie nötig).

\(\sigma\): Die Streuung muss positiv sein, daher kommt keine Normalverteilung in Frage. Eine Exponentialverteilung ist eine von mehreren denkbaren Verteilungen.

Aber welche Werte von lambda kommen in Frage? Das Forschungsteam probiert etwas herum:

qexp(p = .5, rate = 1)

[1] 0.6931472

Mit \(\lambda = 1\) liegt der Median der Streuung der Pupillengrößen (p = .5) bei ca. 0.7 mm. Dieser Wert erscheint etwas klein.

qexp(p = .5, rate = 0.2)

[1] 3.465736

Hm. Eine Streuung der Pupillengrößen von ca. 3.5mm (im Median) um den Mittelwert herum; das könnte passen.

Die große Stichprobe wird den Priori-Wert vermutlich überstimmen, überlegt das Forschungsteam (und hat im Großen und Ganzen Recht).

Die Modelle wie stan_glm() tun sich leichter, wenn man nur die relevanten Daten, ohne fehlende Werte und schon schön fertig vorverarbeitet, zur Analyse in die Modellberechnung gibt:

d3 <-
  d %>% 
  select(size) %>% 
  drop_na()

Die Posteriori-Verteilung kann man mit dem Paket {rstanarm} d.h. mit der Funktion stan_glm() berechnen:

library(rstanarm)
m_pupil <- stan_glm(size ~ 1,
                    data = d3,
                    seed = 42,
                    refresh = 0)

Die Daten sind groß, es kann ein paar Sekunden brauchen…

Hier ist eine nützliche Zusammenfassung der Post-Verteilung.

parameters(m_pupil)

Parameter	Median	CI	CI_low	CI_high	pd	Rhat	ESS	Prior_Distribution	Prior_Location	Prior_Scale
(Intercept)	10.00556	0.95	9.958776	10.05205	1	1.000263	1991.503	normal	10.00562	12.76788

Hier eine Visualisierung der Parameter:

plot(parameters(m_pupil), show_intercept = TRUE)

Natürlich kann man auch die Post-Verteilung plotten (z.B: HDI):

m_hdi <- hdi(m_pupil, ci = c(0.5, 0.95))

plot(m_hdi, show_intercept = TRUE)  # Im Default wird der Intercept nicht gezeigt

Hier zur Info die ersten paar Zeilen des Post-Verteilung:

(Intercept)	sigma
10.04	5.13
10.00	5.07
9.99	5.08
10.00	5.08
9.99	5.11

Das Forschungsteam lässt sich folgende Statistiken zum Modell ausgeben:

eti(m_pupil)

Parameter	CI	CI_low	CI_high	Effects	Component
(Intercept)	0.95	9.958776	10.05205	fixed	conditional

hdi(m_pupil)

Parameter	CI	CI_low	CI_high	Effects	Component
(Intercept)	0.95	9.95854	10.05182	fixed	conditional

eti(m_pupil, ci = .89)

Parameter	CI	CI_low	CI_high	Effects	Component
(Intercept)	0.89	9.966516	10.04339	fixed	conditional

prior_summary(m_pupil)

Priors for model 'm_pupil' 
------
Intercept (after predictors centered)
  Specified prior:
    ~ normal(location = 10, scale = 2.5)
  Adjusted prior:
    ~ normal(location = 10, scale = 13)

Auxiliary (sigma)
  Specified prior:
    ~ exponential(rate = 1)
  Adjusted prior:
    ~ exponential(rate = 0.2)
------
See help('prior_summary.stanreg') for more details

Aufgabe Berichten Sie die Breite des 95% Equal Tails Intervalls.

Hinweise:

• Verwenden Sie die Defaults von rstanarm für Ihr Modell.
• Beachten Sie die üblichen Hinweise des Datenwerks
• Falls Sie Teile der Aufgabe nicht lösen können, weil Ihnen der Stoff dazu fehlt: Einfach ignorieren 😄.

Lösung

Man kann die Breite des 95% ETI aus der Ausgabe von eti() ablesen:

[1] 0.09327839

eti() verwendet im Default eine 95%-Intervall.

---

Categories:

• probability
• bayes
• regression
• string