---
title: "Estimating trial-level effects"
output: rmarkdown::html_vignette
bibliography: citations.bib
csl: apa.csl
link-citations: true
vignette: >
  %\VignetteIndexEntry{Estimating trial-level effects}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<!-- categorical.Rmd is generated from _categorical.Rmd.src Please edit that file -->




## Introduction
By default, the `hmetad` package uses aggregated data (i.e., counts of the  number of trials with the same stimulus, type 1 response, and type 2 response.) This is because data aggregation makes model fitting and simulation much more efficient. But sometimes researchers will be interested in trial-level effects.

One common example would be what is often called "crossed random effects". For example, in a design where all participants make responses to the same set of items, a researcher might want to estimate both participant-level and item-level effects on their model parameters.

We can simulate data from such a design like so:

``` r
library(tidyverse)
library(tidybayes)
library(hmetad)

## average model parameters
K <- 3 ## number of confidence levels
mu_log_M <- -0.5
mu_dprime <- 1.5
mu_c <- 0
mu_c2_0 <- rep(-1, K - 1)
mu_c2_1 <- rep(-1, K - 1)

## participant-level standard deviations
sd_log_M_participant <- 0.25
sd_dprime_participant <- 0.5
sd_c_participant <- 0.33
sd_c2_0_participant <- cov_matrix(rep(0.25, K - 1), diag(K - 1))
sd_c2_1_participant <- cov_matrix(rep(0.25, K - 1), diag(K - 1))

## item-level standard deviations
sd_log_M_item <- 0.1
sd_dprime_item <- 0.5
sd_c_item <- 0.75
sd_c2_0_item <- cov_matrix(rep(0.1, K - 1), diag(K - 1))
sd_c2_1_item <- cov_matrix(rep(0.1, K - 1), diag(K - 1))


## simulate data
d <- expand_grid(
  participant = 1:50,
  item = 1:25
) |>
  ## simulate participant-level differences
  group_by(participant) |>
  mutate(
    z_log_M_participant = rnorm(1, sd = sd_log_M_participant),
    z_dprime_participant = rnorm(1, sd = sd_dprime_participant),
    z_c_participant = rnorm(1, sd = sd_c_participant),
    z_c2_0_participant = list(rmulti_normal(1, mu = rep(0, K - 1), Sigma = sd_c2_0_participant)),
    z_c2_1_participant = list(rmulti_normal(1, mu = rep(0, K - 1), Sigma = sd_c2_1_participant))
  ) |>
  ## simulate item-level differences
  group_by(item) |>
  mutate(
    z_log_M_item = rnorm(1, sd = sd_log_M_item),
    z_dprime_item = rnorm(1, sd = sd_dprime_item),
    z_c_item = rnorm(1, sd = sd_c_item),
    z_c2_0_item = list(rmulti_normal(1, mu = rep(0, K - 1), Sigma = sd_c2_0_item)),
    z_c2_1_item = list(rmulti_normal(1, mu = rep(0, K - 1), Sigma = sd_c2_1_item))
  ) |>
  ungroup() |>
  ## compute model parameters
  mutate(
    log_M = mu_log_M + z_log_M_participant + z_log_M_item,
    dprime = mu_dprime + z_dprime_participant + z_dprime_item,
    c = mu_c + z_c_participant + z_c_item,
    c2_0_diff = map2(
      z_c2_0_participant, z_c2_0_item,
      ~ exp(mu_c2_0 + .x + .y)
    ),
    c2_1_diff = map2(
      z_c2_1_participant, z_c2_1_item,
      ~ exp(mu_c2_1 + .x + .y)
    )
  ) |>
  ## simulate two trials per participant/item (stimulus = 0 and stimulus = 1)
  mutate(trial = pmap(list(dprime, c, log_M, c2_0_diff, c2_1_diff), sim_metad, N_trials = 2)) |>
  select(participant, item, trial) |>
  unnest(trial)
```


```
#> # A tibble: 2,500 × 16
#>    participant  item trial stimulus response correct confidence dprime      c meta_dprime      M
#>          <int> <int> <int>    <int>    <int>   <int>      <int>  <dbl>  <dbl>       <dbl>  <dbl>
#>  1           1     1     1        0        1       0          3  1.25  -0.672      -0.326 -0.261
#>  2           1     1     1        1        1       1          3  1.25  -0.672      -0.326 -0.261
#>  3           1     2     1        0        0       1          2  1.03  -0.988      -0.246 -0.239
#>  4           1     2     1        1        1       1          3  1.03  -0.988      -0.246 -0.239
#>  5           1     3     1        0        1       0          3  0.604 -0.305      -0.203 -0.336
#>  6           1     3     1        1        1       1          3  0.604 -0.305      -0.203 -0.336
#>  7           1     4     1        0        0       1          1  1.97   0.786      -0.479 -0.244
#>  8           1     4     1        1        0       0          3  1.97   0.786      -0.479 -0.244
#>  9           1     5     1        0        0       1          2  1.16   0.479      -0.349 -0.300
#> 10           1     5     1        1        1       1          1  1.16   0.479      -0.349 -0.300
#> # ℹ 2,490 more rows
#> # ℹ 5 more variables: meta_c2_0 <list>, meta_c2_1 <list>, theta <dbl>, theta_1 <dbl>, theta_2 <dbl>
```

Don't worry about the details of the simulation code- what matters is that we have a data set with repeated measures for participants:


``` r
count(d, participant)
#> # A tibble: 50 × 2
#>    participant     n
#>          <int> <int>
#>  1           1    50
#>  2           2    50
#>  3           3    50
#>  4           4    50
#>  5           5    50
#>  6           6    50
#>  7           7    50
#>  8           8    50
#>  9           9    50
#> 10          10    50
#> # ℹ 40 more rows
```
And repeated measures for items:

``` r
count(d, item)
#> # A tibble: 25 × 2
#>     item     n
#>    <int> <int>
#>  1     1   100
#>  2     2   100
#>  3     3   100
#>  4     4   100
#>  5     5   100
#>  6     6   100
#>  7     7   100
#>  8     8   100
#>  9     9   100
#> 10    10   100
#> # ℹ 15 more rows
```

## Standard model with data aggregation
If we would like, we can use the `fit_metad` function on this data with participant-level and item-level effects. However, if we aggregate the data ourselves, we can see that the aggregation doesn't really help us here:


``` r
aggregate_metad(d, participant, item)
#> `hmetad` has inferred that there are K=3 confidence levels in the data. If this is incorrect, please set this manually using the argument `K=<K>`
#> # A tibble: 1,250 × 5
#>    participant  item   N_0   N_1 N[,"N_0_1"] [,"N_0_2"] [,"N_0_3"] [,"N_0_4"] [,"N_0_5"] [,"N_0_6"]
#>          <int> <int> <int> <int>       <int>      <int>      <int>      <int>      <int>      <int>
#>  1           1     1     1     1           0          0          0          0          0          1
#>  2           1     2     1     1           0          1          0          0          0          0
#>  3           1     3     1     1           0          0          0          0          0          1
#>  4           1     4     1     1           0          0          1          0          0          0
#>  5           1     5     1     1           0          1          0          0          0          0
#>  6           1     6     1     1           0          0          1          0          0          0
#>  7           1     7     1     1           1          0          0          0          0          0
#>  8           1     8     1     1           0          0          1          0          0          0
#>  9           1     9     1     1           0          0          0          0          0          1
#> 10           1    10     1     1           0          0          1          0          0          0
#> # ℹ 1,240 more rows
#> # ℹ 1 more variable: N[7:12] <int>
```

As you can see, the aggregated data set has 1250 rows (with two observations per row), which is not much smaller than the trial-level data that we started with! So, in this case, it will probably be easier *not* to aggregate our data. Nevertheless, there is nothing stopping us from fitting the model like normal:^[Note that in practice, fitting hierarchical models will usually require setting informed priors and adjusting the Stan sampler settings.]

``` r
# Priors are chosen arbitrarily for this example.
# Please choose your own wisely!
priors <- prior(normal(0, .25), class = Intercept) +
  set_prior(
    "normal(0, .25)",
    class = "Intercept",
    dpar = c("dprime", "c")
  ) +
  set_prior("normal(-0.5, .1)", class = "Intercept", dpar = metac2_parameters(K = 3)) +
  prior(normal(0, 1), class = sd) +
  set_prior("normal(0, 1)", class = "sd", dpar = c("dprime", "c", metac2_parameters(K = 3)))

m.multinomial <- fit_metad(
  bf(
    N ~ 1 + (1 | participant) + (1 | item),
    dprime + c +
      metac2zero1diff + metac2zero2diff +
      metac2one1diff + metac2one2diff ~
      1 + (1 | participant) + (1 | item)
  ),
  data = d, init = 0, cores = 4, prior = priors
)
```


```
#>  Family: metad__3__normal__absolute__multinomial 
#>   Links: mu = log; dprime = identity; c = identity; metac2zero1diff = log; metac2zero2diff = log; metac2one1diff = log; metac2one2diff = log 
#> Formula: N ~ 1 + (1 | participant) + (1 | item) 
#>          dprime ~ 1 + (1 | participant) + (1 | item)
#>          c ~ 1 + (1 | participant) + (1 | item)
#>          metac2zero1diff ~ 1 + (1 | participant) + (1 | item)
#>          metac2zero2diff ~ 1 + (1 | participant) + (1 | item)
#>          metac2one1diff ~ 1 + (1 | participant) + (1 | item)
#>          metac2one2diff ~ 1 + (1 | participant) + (1 | item)
#>    Data: data.aggregated (Number of observations: 1250) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Multilevel Hyperparameters:
#> ~item (Number of levels: 25) 
#>                               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)                     3.77      0.61     2.63     5.03 1.00     3090     2051
#> sd(dprime_Intercept)              0.62      0.16     0.37     0.99 1.01      818     1551
#> sd(c_Intercept)                   0.60      0.09     0.44     0.82 1.01      705     1089
#> sd(metac2zero1diff_Intercept)     0.24      0.10     0.05     0.46 1.01      623      414
#> sd(metac2zero2diff_Intercept)     0.07      0.06     0.00     0.21 1.00     1406     1055
#> sd(metac2one1diff_Intercept)      0.09      0.07     0.00     0.25 1.00     1348     2050
#> sd(metac2one2diff_Intercept)      0.11      0.07     0.00     0.28 1.00     1201     1405
#> 
#> ~participant (Number of levels: 50) 
#>                               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)                     0.33      0.26     0.01     0.96 1.00     2344     1967
#> sd(dprime_Intercept)              0.34      0.11     0.12     0.55 1.01      747      448
#> sd(c_Intercept)                   0.28      0.04     0.20     0.37 1.00     1394     2236
#> sd(metac2zero1diff_Intercept)     0.28      0.09     0.08     0.45 1.01      445      224
#> sd(metac2zero2diff_Intercept)     0.13      0.08     0.01     0.31 1.00      808     1241
#> sd(metac2one1diff_Intercept)      0.31      0.09     0.13     0.48 1.00      890      994
#> sd(metac2one2diff_Intercept)      0.21      0.09     0.03     0.38 1.00      759      738
#> 
#> Regression Coefficients:
#>                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept                    -0.47      0.26    -0.96     0.04 1.00     3426     3142
#> dprime_Intercept              1.16      0.16     0.80     1.45 1.00      743     1227
#> c_Intercept                  -0.16      0.11    -0.37     0.07 1.00      493      877
#> metac2zero1diff_Intercept    -0.78      0.07    -0.91    -0.64 1.01     1568     2256
#> metac2zero2diff_Intercept    -0.76      0.05    -0.86    -0.67 1.00     3849     3025
#> metac2one1diff_Intercept     -0.84      0.06    -0.96    -0.72 1.00     2262     2593
#> metac2one2diff_Intercept     -0.76      0.06    -0.86    -0.64 1.00     2839     2443
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
```


## Data preparation
Fitting the trial-level model does not require data aggregation, however it still requires a small amount of data preparation. To fit the model, we will need two things:

  * a column with the stimulus per trial (`0` or `1`), and
  * a column containing the joint type 1/type 2 responses per trial (between `1` and `2*K`).

Our data already has a `stimulus` column but separate columns for the two responses. So, we can add in a joint response column now:

``` r
d <- d |>
  mutate(joint_response = joint_response(response, confidence, K)) |>
  relocate(joint_response, .after = "stimulus")
```


```
#> # A tibble: 2,500 × 17
#>    participant  item trial stimulus joint_response response correct confidence dprime      c
#>          <int> <int> <int>    <int>          <int>    <int>   <int>      <int>  <dbl>  <dbl>
#>  1           1     1     1        0              6        1       0          3  1.25  -0.672
#>  2           1     1     1        1              6        1       1          3  1.25  -0.672
#>  3           1     2     1        0              2        0       1          2  1.03  -0.988
#>  4           1     2     1        1              6        1       1          3  1.03  -0.988
#>  5           1     3     1        0              6        1       0          3  0.604 -0.305
#>  6           1     3     1        1              6        1       1          3  0.604 -0.305
#>  7           1     4     1        0              3        0       1          1  1.97   0.786
#>  8           1     4     1        1              1        0       0          3  1.97   0.786
#>  9           1     5     1        0              2        0       1          2  1.16   0.479
#> 10           1     5     1        1              4        1       1          1  1.16   0.479
#> # ℹ 2,490 more rows
#> # ℹ 7 more variables: meta_dprime <dbl>, M <dbl>, meta_c2_0 <list>, meta_c2_1 <list>, theta <dbl>,
#> #   theta_1 <dbl>, theta_2 <dbl>
```


## Model fitting
Now that we have our data, we can fit the trial-level model using `joint_response` as our response variable, `stimulus` as an extra variable passed to `brms`, and the argument `categorical=TRUE` to tell `fit_metad` not to aggregate the data:

``` r
m.categorical <- fit_metad(
  bf(
    joint_response | vint(stimulus) ~ 1 + (1 | participant) + (1 | item),
    dprime + c +
      metac2zero1diff + metac2zero2diff +
      metac2one1diff + metac2one2diff ~
      1 + (1 | participant) + (1 | item)
  ),
  data = d, categorical = TRUE, init = 0, cores = 4, prior = priors
)
```


```
#>  Family: metad__3__normal__absolute__categorical 
#>   Links: mu = log; dprime = identity; c = identity; metac2zero1diff = log; metac2zero2diff = log; metac2one1diff = log; metac2one2diff = log 
#> Formula: joint_response | vint(stimulus) ~ 1 + (1 | participant) + (1 | item) 
#>          dprime ~ 1 + (1 | participant) + (1 | item)
#>          c ~ 1 + (1 | participant) + (1 | item)
#>          metac2zero1diff ~ 1 + (1 | participant) + (1 | item)
#>          metac2zero2diff ~ 1 + (1 | participant) + (1 | item)
#>          metac2one1diff ~ 1 + (1 | participant) + (1 | item)
#>          metac2one2diff ~ 1 + (1 | participant) + (1 | item)
#>    Data: data.aggregated (Number of observations: 2500) 
#>   Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
#>          total post-warmup draws = 4000
#> 
#> Multilevel Hyperparameters:
#> ~item (Number of levels: 25) 
#>                               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)                     3.77      0.63     2.59     5.08 1.00     3879     2953
#> sd(dprime_Intercept)              0.63      0.16     0.36     1.00 1.01      725     1359
#> sd(c_Intercept)                   0.60      0.10     0.44     0.82 1.01      802     1648
#> sd(metac2zero1diff_Intercept)     0.24      0.10     0.04     0.44 1.00      703      722
#> sd(metac2zero2diff_Intercept)     0.07      0.06     0.00     0.21 1.00     1862     2128
#> sd(metac2one1diff_Intercept)      0.09      0.07     0.00     0.25 1.00     1105     1535
#> sd(metac2one2diff_Intercept)      0.10      0.07     0.00     0.27 1.00     1320     1739
#> 
#> ~participant (Number of levels: 50) 
#>                               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sd(Intercept)                     0.33      0.26     0.01     0.96 1.00     2513     1621
#> sd(dprime_Intercept)              0.34      0.10     0.12     0.54 1.00      967      659
#> sd(c_Intercept)                   0.28      0.04     0.20     0.36 1.00     1643     2665
#> sd(metac2zero1diff_Intercept)     0.29      0.08     0.13     0.45 1.00     1046     1038
#> sd(metac2zero2diff_Intercept)     0.13      0.08     0.01     0.31 1.01     1023     1512
#> sd(metac2one1diff_Intercept)      0.31      0.09     0.15     0.49 1.00      922      651
#> sd(metac2one2diff_Intercept)      0.21      0.09     0.02     0.38 1.01      762      882
#> 
#> Regression Coefficients:
#>                           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> Intercept                    -0.47      0.25    -0.98     0.03 1.00     4213     2674
#> dprime_Intercept              1.15      0.17     0.77     1.43 1.01      874     1592
#> c_Intercept                  -0.15      0.12    -0.38     0.08 1.00      608     1119
#> metac2zero1diff_Intercept    -0.79      0.07    -0.92    -0.65 1.00     2215     2501
#> metac2zero2diff_Intercept    -0.76      0.05    -0.86    -0.66 1.00     4528     2893
#> metac2one1diff_Intercept     -0.84      0.06    -0.96    -0.72 1.00     2906     2269
#> metac2one2diff_Intercept     -0.76      0.05    -0.86    -0.65 1.00     3840     3076
#> 
#> Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
#> and Tail_ESS are effective sample size measures, and Rhat is the potential
#> scale reduction factor on split chains (at convergence, Rhat = 1).
```

As you can see, aside from the way that the data is formatted, this model is exactly the same as the multinomial model above.





## Extracting model estimates
Obtaining posterior estimates over model parameters, predictions, and other estimates is very similar to the multinomial model (for details, see `vignette("hmetad")`). So, here we will focus on the type 1 ROC curves, this time using `roc1_rvars` instead of `roc1_draws` for increased efficiency.

To get the posterior estimates, we need to specify a dataset to make predictions for, as well as a random effects formula to use in model predictions. For example, to estimate a ROC averaging over participants and items, we can use an empty data set with `re_formula=NA`:

``` r
roc1_rvars(m.multinomial, tibble(.row = 1), re_formula = NA)
#> # A tibble: 5 × 6
#> # Groups:   .row, joint_response, response, confidence [5]
#>    .row joint_response response confidence          p_fa         p_hit
#>   <int>          <int>    <int>      <int>    <rvar[1d]>    <rvar[1d]>
#> 1     1              1        0          3  0.73 ± 0.046  0.94 ± 0.016
#> 2     1              2        0          2  0.54 ± 0.053  0.87 ± 0.029
#> 3     1              3        0          1  0.34 ± 0.050  0.77 ± 0.042
#> 4     1              4        1          1  0.21 ± 0.040  0.59 ± 0.051
#> 5     1              5        1          2  0.11 ± 0.027  0.39 ± 0.053
```

The process is exactly the same for the categorical model:

``` r
roc1_rvars(m.categorical, tibble(.row = 1), re_formula = NA)
#> # A tibble: 5 × 6
#> # Groups:   .row, joint_response, response, confidence [5]
#>    .row joint_response response confidence          p_fa         p_hit
#>   <int>          <int>    <int>      <int>    <rvar[1d]>    <rvar[1d]>
#> 1     1              1        0          3  0.73 ± 0.049  0.94 ± 0.017
#> 2     1              2        0          2  0.54 ± 0.057  0.87 ± 0.031
#> 3     1              3        0          1  0.34 ± 0.052  0.76 ± 0.044
#> 4     1              4        1          1  0.21 ± 0.042  0.59 ± 0.053
#> 5     1              5        1          2  0.11 ± 0.028  0.39 ± 0.054
```

Next, to get participant-level ROCs (averaging over items), we can use a dataset with one row per participant and only the participant-level random effects:

``` r
roc1_rvars(m.categorical, distinct(d, participant), re_formula = ~ (1 | participant))
#> # A tibble: 250 × 7
#> # Groups:   .row, participant, joint_response, response, confidence [250]
#>     .row participant joint_response response confidence          p_fa         p_hit
#>    <int>       <int>          <int>    <int>      <int>    <rvar[1d]>    <rvar[1d]>
#>  1     1           1              1        0          3  0.71 ± 0.086  0.90 ± 0.042
#>  2     2           2              1        0          3  0.60 ± 0.097  0.84 ± 0.059
#>  3     3           3              1        0          3  0.75 ± 0.079  0.95 ± 0.027
#>  4     4           4              1        0          3  0.81 ± 0.071  0.97 ± 0.018
#>  5     5           5              1        0          3  0.87 ± 0.056  0.98 ± 0.013
#>  6     6           6              1        0          3  0.78 ± 0.076  0.95 ± 0.024
#>  7     7           7              1        0          3  0.71 ± 0.084  0.91 ± 0.041
#>  8     8           8              1        0          3  0.77 ± 0.087  0.97 ± 0.017
#>  9     9           9              1        0          3  0.67 ± 0.089  0.89 ± 0.048
#> 10    10          10              1        0          3  0.64 ± 0.095  0.92 ± 0.038
#> # ℹ 240 more rows
```


We can use a similar process to get item-level ROCs (averaging over participants):

``` r
roc1_rvars(m.categorical, distinct(d, item), re_formula = ~ (1 | item))
#> # A tibble: 125 × 7
#> # Groups:   .row, item, joint_response, response, confidence [125]
#>     .row  item joint_response response confidence          p_fa          p_hit
#>    <int> <int>          <int>    <int>      <int>    <rvar[1d]>     <rvar[1d]>
#>  1     1     1              1        0          3  0.89 ± 0.028  0.99 ± 0.0066
#>  2     2     2              1        0          3  0.91 ± 0.024  1.00 ± 0.0034
#>  3     3     3              1        0          3  0.82 ± 0.037  0.95 ± 0.0168
#>  4     4     4              1        0          3  0.50 ± 0.055  0.80 ± 0.0434
#>  5     5     5              1        0          3  0.63 ± 0.051  0.89 ± 0.0309
#>  6     6     6              1        0          3  0.89 ± 0.028  0.97 ± 0.0116
#>  7     7     7              1        0          3  0.60 ± 0.053  0.83 ± 0.0393
#>  8     8     8              1        0          3  0.69 ± 0.049  0.95 ± 0.0183
#>  9     9     9              1        0          3  0.87 ± 0.031  0.98 ± 0.0095
#> 10    10    10              1        0          3  0.79 ± 0.040  0.93 ± 0.0218
#> # ℹ 115 more rows
```


## Other benefits
Aside from representing the data in a more convenient format, the trial-level model should be more useful for things like model comparison using the `loo` package, multivariate models, and mediation models. These features should mostly work out of the box but they are still under active development, so stay tuned!