class: center, middle, inverse, title-slide # Datenanalyse ## Teil 6 <br/> Binäre Daten ### Andrew Ellis ### Methodenkurs Neurowissenschaft im Computerlab ### April 27, 2021 --- layout: true <!-- Home icon --> <div class="my-footer"> <span> <a href="https://kogpsy.github.io/neuroscicomplab" target="_blank"><svg aria-hidden="true" role="img" viewBox="0 0 576 512" style="height:1em;width:1.12em;vertical-align:-0.125em;margin-left:auto;margin-right:auto;font-size:inherit;fill:#0F4C81;overflow:visible;position:relative;"><path d="M280.37 148.26L96 300.11V464a16 16 0 0 0 16 16l112.06-.29a16 16 0 0 0 15.92-16V368a16 16 0 0 1 16-16h64a16 16 0 0 1 16 16v95.64a16 16 0 0 0 16 16.05L464 480a16 16 0 0 0 16-16V300L295.67 148.26a12.19 12.19 0 0 0-15.3 0zM571.6 251.47L488 182.56V44.05a12 12 0 0 0-12-12h-56a12 12 0 0 0-12 12v72.61L318.47 43a48 48 0 0 0-61 0L4.34 251.47a12 12 0 0 0-1.6 16.9l25.5 31A12 12 0 0 0 45.15 301l235.22-193.74a12.19 12.19 0 0 1 15.3 0L530.9 301a12 12 0 0 0 16.9-1.6l25.5-31a12 12 0 0 0-1.7-16.93z"/></svg></a> Methodenkurs Neurowissenschaft im Computerlab </span> </div> <!-- Name (left) --> <!-- <div class="my-footer"> --> <!-- <span> --> <!-- Andrew Ellis - <a href="https://kogpsy.github.io/neuroscicomplab" target="_blank">kogpsy.github.io/neuroscicomplab</a> --> <!-- </span> --> <!-- </div> --> <!-- slide separator (for xaringan) --> --- ## Binäre Daten Sehr häufig anzutreffen in der experimentellen Psychologie/Neuroscience - Entscheidungsexperimente (von Neuroeconomics bis Recognition Memory) - Psychophysik (heller/dunkler, links/rechts) - Accuracy (richtig/falsch) Wir nehmen an, Antworten `\(y_{ijk}\)` erfolgen mit einer bestimmten Wahrscheinlichkeit `\(\theta_{ijk}\)`. Diese hängt wiederum ab von experimentellen Manipulationen und personenspezifischen Faktoren. `\(y_{ijk}\)`: `\(i\)`-te Antwort von Person `\(k\)` in Bedingung `\(j\)`. Oft werden auch Reaktionszeiten miterhoben. Diese schauen wir uns später an. --- ## Beipiel [Übung 3, Aufgabe 3](https://kogpsy.github.io/neuroscicomplab/solution-3.html#aufgabe-3). .panelset[ .panel[.panel-name[Studie] In einem Experiment von Ganis and Kievit (2015) mussten 54 Vpn entscheiden ob zwei 3-D Stimuli, welche zuneinander um 0, 50, 100, oder 150 Grad rotiert waren, identisch oder nicht waren. ] .panel[.panel-name[Variablen] - id: Versuchspersonennummer - angle: Rotationswinkel - rt: Reaktionszeiten - accuracy: Indikator für korrekte Antwort ] .panel[.panel-name[Fragestellung] Sie interessieren sich dafür, ob die Personen bei steigendem Rotationswinkel länger brauchen, um eine korrekte Antwort zu geben, und ob die Personen bei steigendem Rotationswinkel mehr Fehler machen. ] .panel[.panel-name[Daten] `\(y_{ijk}\)`: Accuracy (Korrektheit) der Antwort `\(i\)` von Person `\(k\)` bei Winkel `\(j\)`. ]] --- .panelset[ .panel[.panel-name[Daten laden] ```r library(tidyverse) d3 <- read_csv("https://raw.githubusercontent.com/kogpsy/neuroscicomplab/main/data/mental-rotation.csv") d3 <- d3 %>% mutate(id = as_factor(id)) ``` ] .panel[.panel-name[Daten anschauen] ```r d3 ``` ``` ## # A tibble: 5,184 x 4 ## id angle rt accuracy ## <fct> <dbl> <dbl> <dbl> ## 1 1 0 1355 1 ## 2 1 150 2079 1 ## 3 1 150 1834 1 ## 4 1 100 4780 1 ## 5 1 50 1685 1 ## 6 1 50 1237 1 ## 7 1 100 2961 1 ## 8 1 0 1275 1 ## 9 1 150 4716 1 ## 10 1 150 3794 0 ## # … with 5,174 more rows ``` ]] --- ## Traditionelle Analyse - Two-level model + Level 1: Individuelle Personen + Level 2: Gruppenebene - Daten werden zuerst pro Person pro Bedingung zusammengefasst. - Danach: Group-level Statistik mit den über Trials aggregierten Daten. --- ## Level 1: Individual level .panelset[ .panel[.panel-name[Level 1] ```r rts <- d3 %>% group_by(id, angle, accuracy) %>% summarise(rt = mean(rt)) %>% filter(accuracy == 1) %>% select(-accuracy) acc <- d3 %>% group_by(id, angle) %>% summarise(accuracy = round(mean(accuracy), 2)) d3sum <- acc %>% left_join(rts) %>% mutate(angle = as_factor(angle)) ``` ] .panel[.panel-name[Level 1 Daten] ```r d3sum ``` ``` ## # A tibble: 216 x 4 ## # Groups: id [54] ## id angle accuracy rt ## <fct> <fct> <dbl> <dbl> ## 1 1 0 1 1512. ## 2 1 50 0.83 1977. ## 3 1 100 0.83 2771. ## 4 1 150 0.88 2931. ## 5 2 0 1 1939. ## 6 2 50 0.88 3044. ## 7 2 100 0.71 4740. ## 8 2 150 0.54 3927. ## 9 3 0 1 1808. ## 10 3 50 1 2689. ## # … with 206 more rows ``` ]] --- ## Level 2: Group level .panelset[ .panel[.panel-name[Deskriptiv] ```r se <- function(x) sd(x)/sqrt(length(x)) d3sum_agg <- d3sum %>% group_by(angle) %>% summarise(mean = mean(accuracy), sd = sd(accuracy), se = se(accuracy)) d3sum_agg ``` ``` ## # A tibble: 4 x 4 ## angle mean sd se ## <fct> <dbl> <dbl> <dbl> ## 1 0 0.948 0.0786 0.0107 ## 2 50 0.928 0.0648 0.00882 ## 3 100 0.867 0.102 0.0139 ## 4 150 0.813 0.122 0.0166 ``` ] .panel[.panel-name[Grafik Code] ```r d3sum_agg %>% ggplot(aes(angle, mean)) + geom_line(aes(group = 1), linetype = 3) + geom_errorbar(aes(ymin = mean-se, ymax = mean+se), width = 0.2, size=1, color="blue") + geom_point(size = 4) + theme_bw(base_size = 14) ``` ] .panel[.panel-name[Grafik] <img src="06-choices-slides_files/figure-html/unnamed-chunk-8-1.png" width="100%" /> ] .panel[.panel-name[Level 2 Modell] ```r library(brms) m_l2 <- brm(accuracy ~ angle + (angle | id), data = d3sum, file = "models/06-m_l2", file_refit = "on_change") ``` ] .panel[.panel-name[Level 2 Ouput] ```r fixef(m_l2) ``` ``` ## Estimate Est.Error Q2.5 Q97.5 ## Intercept 0.94739656 0.01045058 0.9270760 0.968323078 ## angle50 -0.01917463 0.01311616 -0.0454511 0.006417585 ## angle100 -0.08048942 0.01506370 -0.1104279 -0.050822549 ## angle150 -0.13445853 0.01858680 -0.1712514 -0.097525181 ``` ] .panel[.panel-name[Level 2 Parameter] ```r m_l2 %>% mcmc_plot("b") ``` <img src="06-choices-slides_files/figure-html/unnamed-chunk-11-1.png" width="100%" /> ] ] --- .your-turn[ Was könnten Vor- und Nachteile dieser 2-Level Analyse sein (dies könnte hier genausogut eine frequentistische Analyse sein)? ] <div class="countdown" id="timer_6087de4c" style="right:0;bottom:0;" data-warnwhen="0"> <code class="countdown-time"><span class="countdown-digits minutes">02</span><span class="countdown-digits colon">:</span><span class="countdown-digits seconds">00</span></code> </div> --- ## Probleme - `accuracy` liegt im Bereich `\([0, 1]\)`. Nahe bei `\(0\)` und `\(1\)` sind die Daten nicht approximativ normalverteilt. - Bei der Level 1 Analyse mitteln wir über Trials: `mean(accuracy)`. Hier wird die Unsicherheit beim Schätzen der Accuracy vernachlässigt. - Dies führt dazu, dass die Parameterschätzung auf Level 2 fälschlicherweise zu wenig unsicher ist, d.h. die Genauigkeit der Schätzung wird überschätzt. .alert[ Moscatelli, A., M. Mezzetti, and F. Lacquaniti. “Modeling Psychophysical Data at the Population-Level: The Generalized Linear Mixed Model.” Journal of Vision 12, no. 11 (October 25, 2012): 26–26. https://doi.org/10.1167/12.11.26. ] --- ## Generalized Linear Model (GLM) - Anstelle der aggregierten Accuracy können wir auch die Antworten selber modellieren + Wir wissen, ob eine Antwort richtig oder falsch ist. + Wir sagen die Wahrscheinlichkeit einer korrekten Antwort mit einem linearen Modell vorher. - Ein lineares Modell: `\(b_0 + b_1 \cdot X_1 + ...\)` hat einen Wertebereich `\([-\infty, \infty]\)`. - Eine Wahrscheinlichkeit liegt in `\([0, 1]\)`. Also muss der lineare Prädiktor so transformiert werden, dass das Resultat auch in `\([0, 1]\)` liegt. - Funktionen, welche so etwas machen: kumulative Verteilungsfunktionen --- ## Generalized Linear Model (GLM) .pull-left[ Die kumulative Verteilungsfunktionen der logistischen Verteilung (eine logistische Funktion) heisst in R `plogis()`. Die Formel lautet $$ F(x; \mu, s) = \frac{1}{1 + e^{-(x-\mu / s)}} $$ ] .pull-right[ ```r theme_set(theme_grey(base_size = 16) + theme(panel.grid = element_blank())) d <- tibble(x = seq(-5, 5, by = 0.01), y = plogis(x)) d %>% ggplot(aes(x, y)) + geom_hline(yintercept = 0.5, linetype = 3) + geom_vline(xintercept = 0, linetype = 3) + geom_line(size = 2) ``` <img src="06-choices-slides_files/figure-html/unnamed-chunk-13-1.png" width="100%" /> ] --- ## Generalized Linear Model (GLM) Vielleicht etwas verwirrenderweise nennt man die Umkehrfunktion die `Link` function. Die inverse logistische Funktion heisst `logit` (oder `log-odds`). `\(p\)` ist eine Wahrscheinlichkeit, welche mit einem linearen Prädiktor vorhergesagt werden soll. $$ logit(p) = b_0 + b_1 \cdot X_1$$ $$ p = logistic(b_0 + b_1 \cdot X_1)$$ Sie kennen ein solches Modell als logistische Regression. --- .panelset[ .panel[.panel-name[Verteilungsfunktion] ```r d1 <- tibble(x = seq(-5, 5, by = 0.01), y = plogis(x)) d1 %>% ggplot(aes(x, y)) + geom_hline(yintercept = 0.5, linetype = 3) + geom_vline(xintercept = 0, linetype = 3) + geom_line(size = 2, color = "steelblue") ``` <img src="06-choices-slides_files/figure-html/unnamed-chunk-14-1.png" width="100%" /> ] .panel[.panel-name[Quantilfunktion] ```r d2 <- tibble(y = seq(0, 1, by = 0.01), x = qlogis(y)) d2 %>% ggplot(aes(y, x)) + geom_hline(yintercept = 0, linetype = 3) + geom_vline(xintercept = 0.5, linetype = 3) + geom_line(size = 2, color = "steelblue") + scale_y_continuous(limits = c(-5, 5)) ``` <img src="06-choices-slides_files/figure-html/unnamed-chunk-15-1.png" width="100%" /> ]] --- ## Mulitlevel Logistische Regression Wir wollen nun die Rohdaten modellieren, d.h. die Anworten. Diese folgen einer Bernoulli-Verteilung: `$$y_{ijk} \sim Bernoulli(\theta_{ijk})$$` `$$logit(\theta_{ijk}) = b_0^{k} + b_{angle}^{k} \cdot angle_{j}$$` oder `$$\theta_{ijk} = logistic(b_0^{k} + b_{angle}^{k} \cdot angle_{j})$$` --- ## Logistische Regression in brms ```r d3 <- d3 %>% mutate(angle = as_factor(angle)) p <- get_prior(accuracy ~ 1 + angle + (1 + angle | id), family = bernoulli(link = logit), data = d3) as_tibble(p) %>% select(1:4) ``` ``` ## # A tibble: 13 x 4 ## prior class coef group ## <chr> <chr> <chr> <chr> ## 1 "" b "" "" ## 2 "" b "angle100" "" ## 3 "" b "angle150" "" ## 4 "" b "angle50" "" ## 5 "lkj(1)" cor "" "" ## 6 "" cor "" "id" ## 7 "student_t(3, 0, 2.5)" Intercept "" "" ## 8 "student_t(3, 0, 2.5)" sd "" "" ## 9 "" sd "" "id" ## 10 "" sd "angle100" "id" ## 11 "" sd "angle150" "id" ## 12 "" sd "angle50" "id" ## 13 "" sd "Intercept" "id" ``` --- ## Logistische Regression in brms Nur die ersten 20 Personen (damit es nicht so lange dauert) ```r d3_20 <- d3 %>% filter(id %in% 1:20) priors <- prior(normal(0, 1), class = b) m1 <- brm(accuracy ~ 1 + angle + (1 + angle | id), family = bernoulli(link = logit), prior = priors, data = d3_20, control = list(adapt_delta = 0.9), file = "models/06-m1", file_refit = "on_change") ``` --- ```r fixef(m1) ``` ``` ## Estimate Est.Error Q2.5 Q97.5 ## Intercept 3.2467205 0.2730034 2.722703 3.79682284 ## angle50 -0.6505253 0.3298081 -1.262615 0.02722595 ## angle100 -1.4767602 0.2665032 -1.994347 -0.95807806 ## angle150 -1.9144632 0.2552711 -2.420739 -1.43763655 ``` --- ```r mcmc_plot(m1, "b") ``` <img src="06-choices-slides_files/figure-html/unnamed-chunk-19-1.png" width="100%" /> --- ```r conditional_effects(m1) ``` <img src="06-choices-slides_files/figure-html/unnamed-chunk-20-1.png" width="100%" /> --- ```r d3_20 %>% filter(id %in% c(11, 13, 15, 17)) %>% ggplot(aes(angle, accuracy)) + geom_violin() + geom_jitter() + facet_wrap(~id) ``` <img src="06-choices-slides_files/figure-html/unnamed-chunk-21-1.png" width="100%" />