Übung 6

DDM Parameter schätzen.

Published

2022-05-10

Note

Die Daten für diese Übung finden Sie hier:

👉 Download Data

👉 Download True Parameters

Laden Sie bitte Ihre Lösung als Word Dokument oder R Skript bis Dienstag, 24.5.2022, um 00:30 Uhr, in den Order für Übung 6 auf ILIAS.

Nennen Sie Ihr File Matrikelnummer_Nachname_uebung-6.R oder Matrikelnummer_Nachname_uebung-6.docx.

Aufgabenstellung

In dieser Aufgabe geht es darum, dass Sie einmal selber die Parameter eines DDM an Daten fitten können. Es geht nicht darum, jede Zeile Code zu verstehen; es geht vor allem darum, zu erleben, wie schwierig es sein kann, solche Modelle an Daten zu fitten.

Bitte führen Sie den Code, in dem die Parameter geschätzt werden, mehrmals aus, und kopieren Sie die geschätzten Parameter in ein File. Dies muss abgegeben werden.

Daten einlesen

Downloaden Sie die Daten, und die wahren Parameter (zum Vergleich).

Wir laden die Daten in ein DataFrame, d. Ich habe ein RT Experiment simuliert, in dem 4 Versuchspersonen in 2 Bedingungen (“A” und “B”) getestet wurden. Die experimentelle Manipulation sollte sich vor allem auf einen Parameter unterschieden. In dieser Übung geht es darum, herauszufinden, auf welchen Parameter die Manipulation einen Einfluss hat.

Aufgaben

Important
  1. Führen Sie die Parameterschätzung mehrmals aus, und berichten Sie die resultierenden Parameter. Es kann sein, dass nicht immer dasselbe rauskommt (deshlab machen wir es mehrmals). Welcher Parameter wurde beeinflusst?

  2. Überlegen Sie sich, wie Sie weiterführen würden. Wie können Sie zeigen, dass es in einem Parameter einen Unterschied zwischen den Bedingungen gibt? Wie können Sie zeigen, dass es in den anderen Parameter keine Unterschied gibt? Beschreiben Sie in einem Paragraphen, was Sie sich überlegt haben.

library(tidyverse)
library(rtdists)

d <- read_csv("ddm-data.csv")
participant_params <- read_csv("participant-params.csv")

Daten vorbereiten

Wir schauen uns die Daten an:

d |> glimpse()
Rows: 9,600
Columns: 4
$ ID        <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ condition <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", …
$ rt        <dbl> 0.6888584, 0.3411717, 1.6588621, 0.8845731, 0.7925147, 0.440…
$ response  <chr> "lower", "upper", "upper", "upper", "upper", "upper", "lower…

Die Variablen ID, condition und response sollten Faktoren sein.

d <- d |> 
  mutate(across(c(ID, condition, response), ~as_factor(.)))

d |> glimpse()
Rows: 9,600
Columns: 4
$ ID        <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ condition <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, …
$ rt        <dbl> 0.6888584, 0.3411717, 1.6588621, 0.8845731, 0.7925147, 0.440…
$ response  <fct> lower, upper, upper, upper, upper, upper, lower, upper, uppe…

Weil es nur 4 VP sind, können wir die RTs von allen in einem Plot anschauen. Die upper Responses können wir hier als korrekte Antworten auffassen, die lower Responses als inkorrekte.

d |>
  ggplot(aes(rt, response, fill = response)) +
  geom_violin() +
  geom_jitter(height = 0.1, alpha = 0.2, size = 0.25) +
  scale_fill_viridis_d(option = "B", direction = 1,
                       begin = 1/2, end = 2/2) +
  xlim(c(0, 1.5)) +
  facet_grid(condition ~ ID)
Warning: Removed 488 rows containing non-finite values (stat_ydensity).
Warning: Removed 488 rows containing missing values (geom_point).

Daten zusammenfassen

UM uns einen Überblick zu verschaffen, fassen wir die Daten durch mittlere RT und “Accuracy” zusammen.

summary <- d |> group_by(ID, condition) |> 
  summarise(mean_rt = mean(rt),
            median_rt = median(rt),
            accuracy = mean(response == "upper"))
`summarise()` has grouped output by 'ID'. You can override using the `.groups`
argument.
summary
# A tibble: 8 × 5
# Groups:   ID [4]
  ID    condition mean_rt median_rt accuracy
  <fct> <fct>       <dbl>     <dbl>    <dbl>
1 1     A           0.716     0.589    0.615
2 1     B           0.643     0.536    0.863
3 2     A           0.725     0.591    0.61 
4 2     B           0.617     0.519    0.842
5 3     A           0.843     0.696    0.598
6 3     B           0.685     0.578    0.864
7 4     A           0.715     0.591    0.602
8 4     B           0.617     0.518    0.85 
summary |> 
  ggplot(aes(mean_rt, accuracy, color = condition)) +
  geom_line(aes(group = ID), color = "black", linetype = "dotted") +
  geom_point(size = 4) +
  scale_color_viridis_d(end = 0.8)

In der Grafik sehen wir, dass die Accuracy in Bedingung “B” höher ist. Gleichzeitig ist die mittlere RT niedriger. Das ist ein Hinweis, dass hier kein Speed-Accuracy Tradeoff vorliegt.

Negative loglikelihood Funktion definieren

Nun wollen wir für jede Person in den beiden Bedingungen die Parameter des DDM fitten. Wir definieren dafür eine Funktion, welche die negative log likelihood als Output hat.

diffusionloglik <- function(pars, condition, rt, response) {

  conditions <- levels(condition)
  likelihoods <- vector("numeric", length(rt))


  likelihoods <- ddiffusion(rt = rt,
                             response = response,
                             a = pars["a"],
                             v =  pars["v"],
                             t0 = pars["t0"],
                             z = pars["z"] * pars["a"],
                             s = 1.0)
  
  if (any(likelihoods == 0)) return(1e6)
  return(-sum(log(likelihoods)))
}

Startwerte für Maxmimum Likelihood Schätzung

Für die Minimierung brauchen wir Anfangswerte. Diese werden für jeden zu schätzenden Parameter zufällig gewählt.

init_params <- function() {
  params <- c(a = runif(1, 0.2, 1.2),
              v = rnorm(1, 0.5, 0.5),
              z = runif(1, 0.45, 0.55),
              t0 = runif(1, 0.01, 0.3))
  params
}

Maxmimum Likelihood Schätzung

Nun definieren wir zuerst ein paar Variablen.

participants <- levels(d$ID)
n_participants <- length(participants)
conditions <- levels(d$condition)
n_conditions <- length(conditions)

# no. parameters (a, v, z, t0)
n_pars <- length(init_params())

Und nun schreiben wir eine zweifache for-Loop, über die Versuchspersonen, und über die Bedingungen innerhalb der Personen.

Führen Sie diesen Teil mehrmals aus.

p <- vector("list", n_participants)

for (i in seq_along(participants)) {
  
  estimates <- array(NA, c(n_conditions, n_pars))
  colnames(estimates) <- c("a", "v", "z", "t0")
  rownames(estimates) <- c("A", "B")
  
  for (j in seq_along(conditions)) {
    data <- filter(d, ID == i, condition == conditions[j])
  
    fit <- nlminb(init_params(),
                diffusionloglik,
                lower = 0,
                condition = data$condition,
                rt = data$rt,
                response = data$response)
    
    estimates[j, ] <- fit$par |> round(3)
  }
  p[[i]] <- estimates
}

Geschätzte Parameterwerte

Die geschätzten Parameter sind in einer List gespeichert.

p
[[1]]
      a     v     z    t0
A 1.439 0.261 0.523 0.206
B 1.513 1.231 0.495 0.190

[[2]]
      a     v     z    t0
A 1.462 0.293 0.505 0.200
B 1.427 1.180 0.499 0.202

[[3]]
      a     v     z    t0
A 1.631 0.284 0.484 0.188
B 1.583 1.268 0.469 0.192

[[4]]
      a     v     z    t0
A 0.664 0.259 0.519 0.292
B 0.514 0.751 0.527 0.285

Vergleich mit wahren Werten

Die geschätzten Parameter können nun mit den “wahren” Werten verglichen werden.

participant_params
# A tibble: 4 × 6
     ID     a    v1    v2     z    t0
  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1     1  1.48 0.301  1.19 0.515 0.199
2     2  1.48 0.299  1.23 0.495 0.199
3     3  1.57 0.277  1.21 0.484 0.197
4     4  1.47 0.296  1.23 0.498 0.195

Reuse

Citation

BibTeX citation:
@online{ellis2022,
  author = {Andrew Ellis},
  title = {Übung 6},
  date = {2022-05-10},
  url = {https://kogpsy.github.io/neuroscicomplabFS22//pages/exercises/exercise_06.html},
  langid = {en}
}
For attribution, please cite this work as:
Andrew Ellis. 2022. “Übung 6.” May 10, 2022. https://kogpsy.github.io/neuroscicomplabFS22//pages/exercises/exercise_06.html.