Hello,
I have run a model with brms using the Wiener family following the tutorial by Henrik Singmann. I would like to extract subject-specific drift rate estimates per condition from this model and plot them. My idea was to combine the fixed effects estimates with the random effects estimates.
Here is the brms output of the model provided in the tutorial as a reference:
library("brms")
library("tidyverse")
tmp <- tempdir()
download.file("https://singmann.github.io/files/brms_wiener_example_fit.rda",
file.path(tmp, "brms_wiener_example_fit.rda"))
load(file.path(tmp, "brms_wiener_example_fit.rda"))
summary(fit_wiener)
Family: wiener
Links: mu = identity; bs = identity; ndt = identity; bias = identity
Formula: rt | dec(response2) ~ 0 + condition:frequency + (0 + condition:frequency | p | id)
bs ~ 0 + condition + (0 + condition | p | id)
ndt ~ 0 + condition + (0 + condition | p | id)
bias ~ 0 + condition + (0 + condition | p | id)
Data: speed_acc (Number of observations: 10462)
Draws: 4 chains, each with iter = 1000; warmup = 500; thin = 1;
total post-warmup draws = 2000
Group-Level Effects:
~id (Number of levels: 17)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(conditionaccuracy:frequencyhigh) 0.74 0.15 0.52 1.10 1.00 938 1145
sd(conditionspeed:frequencyhigh) 0.83 0.15 0.59 1.18 1.00 697 1104
sd(conditionaccuracy:frequencynw_high) 0.56 0.10 0.40 0.77 1.00 874 1374
sd(conditionspeed:frequencynw_high) 0.68 0.12 0.49 0.95 1.00 793 1185
sd(bs_conditionaccuracy) 0.60 0.12 0.42 0.87 1.00 1146 1348
sd(bs_conditionspeed) 0.32 0.07 0.21 0.49 1.00 884 1324
sd(ndt_conditionaccuracy) 0.07 0.01 0.05 0.10 1.00 1211 1225
sd(ndt_conditionspeed) 0.06 0.01 0.04 0.09 1.00 696 843
sd(bias_conditionaccuracy) 0.04 0.01 0.02 0.06 1.00 1098 1380
sd(bias_conditionspeed) 0.04 0.01 0.03 0.07 1.00 959 1383
cor(conditionaccuracy:frequencyhigh,conditionspeed:frequencyhigh) 0.33 0.19 -0.07 0.66 1.00 1116 1421
cor(conditionaccuracy:frequencyhigh,conditionaccuracy:frequencynw_high) -0.46 0.17 -0.76 -0.08 1.00 1194 1407
cor(conditionspeed:frequencyhigh,conditionaccuracy:frequencynw_high) -0.46 0.17 -0.75 -0.08 1.00 1243 1467
cor(conditionaccuracy:frequencyhigh,conditionspeed:frequencynw_high) -0.30 0.19 -0.64 0.11 1.00 1210 1328
cor(conditionspeed:frequencyhigh,conditionspeed:frequencynw_high) -0.57 0.17 -0.84 -0.16 1.00 1186 1340
cor(conditionaccuracy:frequencynw_high,conditionspeed:frequencynw_high) 0.63 0.16 0.25 0.87 1.00 1157 1404
cor(conditionaccuracy:frequencyhigh,bs_conditionaccuracy) 0.25 0.19 -0.13 0.59 1.00 1208 1098
cor(conditionspeed:frequencyhigh,bs_conditionaccuracy) -0.04 0.20 -0.40 0.36 1.00 1471 1493
cor(conditionaccuracy:frequencynw_high,bs_conditionaccuracy) 0.15 0.19 -0.23 0.50 1.00 1287 1407
cor(conditionspeed:frequencynw_high,bs_conditionaccuracy) 0.01 0.19 -0.38 0.37 1.00 1026 798
cor(conditionaccuracy:frequencyhigh,bs_conditionspeed) 0.21 0.20 -0.22 0.55 1.00 1052 1195
cor(conditionspeed:frequencyhigh,bs_conditionspeed) -0.18 0.20 -0.54 0.22 1.00 1226 1364
cor(conditionaccuracy:frequencynw_high,bs_conditionspeed) -0.09 0.19 -0.45 0.28 1.00 1241 1234
cor(conditionspeed:frequencynw_high,bs_conditionspeed) 0.12 0.20 -0.28 0.48 1.01 942 1446
cor(bs_conditionaccuracy,bs_conditionspeed) 0.24 0.19 -0.17 0.58 1.00 1324 1549
cor(conditionaccuracy:frequencyhigh,ndt_conditionaccuracy) 0.00 0.19 -0.37 0.38 1.00 1351 1323
cor(conditionspeed:frequencyhigh,ndt_conditionaccuracy) 0.14 0.18 -0.24 0.46 1.00 1409 1539
cor(conditionaccuracy:frequencynw_high,ndt_conditionaccuracy) 0.07 0.18 -0.28 0.42 1.00 1506 1490
cor(conditionspeed:frequencynw_high,ndt_conditionaccuracy) 0.13 0.19 -0.25 0.49 1.00 1149 921
cor(bs_conditionaccuracy,ndt_conditionaccuracy) -0.50 0.16 -0.77 -0.18 1.00 1385 1457
cor(bs_conditionspeed,ndt_conditionaccuracy) 0.05 0.20 -0.32 0.42 1.00 1751 1624
cor(conditionaccuracy:frequencyhigh,ndt_conditionspeed) -0.11 0.20 -0.49 0.28 1.00 1283 1188
cor(conditionspeed:frequencyhigh,ndt_conditionspeed) -0.35 0.18 -0.67 0.05 1.00 1505 1064
cor(conditionaccuracy:frequencynw_high,ndt_conditionspeed) 0.33 0.19 -0.07 0.66 1.00 1226 1135
cor(conditionspeed:frequencynw_high,ndt_conditionspeed) 0.34 0.18 -0.04 0.67 1.00 1480 1419
cor(bs_conditionaccuracy,ndt_conditionspeed) 0.00 0.20 -0.38 0.38 1.00 1379 1319
cor(bs_conditionspeed,ndt_conditionspeed) -0.35 0.18 -0.68 0.01 1.00 1669 1554
cor(ndt_conditionaccuracy,ndt_conditionspeed) 0.08 0.19 -0.29 0.45 1.00 1528 902
cor(conditionaccuracy:frequencyhigh,bias_conditionaccuracy) -0.22 0.21 -0.62 0.21 1.00 1644 1430
cor(conditionspeed:frequencyhigh,bias_conditionaccuracy) 0.15 0.20 -0.26 0.52 1.00 1581 1324
cor(conditionaccuracy:frequencynw_high,bias_conditionaccuracy) -0.36 0.20 -0.70 0.07 1.00 1523 1422
cor(conditionspeed:frequencynw_high,bias_conditionaccuracy) -0.33 0.20 -0.68 0.09 1.00 1635 1793
cor(bs_conditionaccuracy,bias_conditionaccuracy) -0.06 0.21 -0.47 0.36 1.00 1978 1760
cor(bs_conditionspeed,bias_conditionaccuracy) 0.02 0.22 -0.42 0.44 1.00 1788 1587
cor(ndt_conditionaccuracy,bias_conditionaccuracy) -0.34 0.20 -0.67 0.09 1.00 1862 1546
cor(ndt_conditionspeed,bias_conditionaccuracy) -0.25 0.21 -0.63 0.20 1.00 1768 995
cor(conditionaccuracy:frequencyhigh,bias_conditionspeed) 0.17 0.21 -0.26 0.55 1.00 1487 1110
cor(conditionspeed:frequencyhigh,bias_conditionspeed) 0.02 0.21 -0.39 0.42 1.00 1351 1589
cor(conditionaccuracy:frequencynw_high,bias_conditionspeed) -0.30 0.20 -0.66 0.14 1.00 1194 1407
cor(conditionspeed:frequencynw_high,bias_conditionspeed) -0.31 0.20 -0.65 0.11 1.00 1528 1580
cor(bs_conditionaccuracy,bias_conditionspeed) 0.01 0.21 -0.41 0.42 1.00 1712 1402
cor(bs_conditionspeed,bias_conditionspeed) 0.32 0.21 -0.13 0.67 1.00 1752 1693
cor(ndt_conditionaccuracy,bias_conditionspeed) -0.36 0.19 -0.70 0.06 1.00 1820 1552
cor(ndt_conditionspeed,bias_conditionspeed) -0.34 0.20 -0.70 0.08 1.00 1708 1739
cor(bias_conditionaccuracy,bias_conditionspeed) 0.32 0.22 -0.13 0.71 1.00 1329 1357
Population-Level Effects:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
conditionaccuracy:frequencyhigh -2.94 0.20 -3.34 -2.56 1.00 636 834
conditionspeed:frequencyhigh -2.72 0.21 -3.13 -2.30 1.00 710 847
conditionaccuracy:frequencynw_high 2.24 0.14 1.97 2.51 1.00 790 1172
conditionspeed:frequencynw_high 1.99 0.18 1.63 2.33 1.00 849 905
bs_conditionaccuracy 1.90 0.14 1.61 2.19 1.00 933 1198
bs_conditionspeed 1.36 0.08 1.20 1.52 1.00 827 942
ndt_conditionaccuracy 0.32 0.02 0.29 0.36 1.00 1057 924
ndt_conditionspeed 0.26 0.02 0.23 0.29 1.00 1021 1127
bias_conditionaccuracy 0.47 0.01 0.45 0.49 1.00 826 1153
bias_conditionspeed 0.50 0.01 0.47 0.52 1.00 1230 1399
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).
Warning message:
There were 7 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
Here is a minimal example to extract the drift rates based on this:
# extract the fixed effects
df.fix = as.data.frame(fixef(fit_wiener)) %>%
rownames_to_column(., var = "condition") %>%
mutate(
condition = if_else(
substr(condition,1,4) == "cond",
paste0("drift_", condition), condition
)
) %>%
separate(col = "condition", into = c("dpar", "condition", "frequency")) %>%
rename("fix_est" = "Estimate") %>%
filter(dpar == "drift")
# extract the random effects
ran.ddm = ranef(fit_wiener)
# convert to data frame
df.ran = data.frame(matrix(ran.ddm[['id']],
ncol = dim(ran.ddm[['id']])[2],
byrow=TRUE))
# add the colnames
colnames(df.ran) = colnames(ran.ddm[['id']])
# merge everything together
df.sub = df.ran %>%
select(Estimate) %>%
# add the coef and dpar
mutate(
dpar = rep(fit_wiener[["ranef"]][["dpar"]], each = dim(ran.ddm[['id']])[1]),
dpar = if_else(dpar == "", "drift", dpar),
coef = rep(fit_wiener[["ranef"]][["coef"]], each = dim(ran.ddm[['id']])[1]),
subID = rep(unique(fit_wiener[["data"]][["id"]]), times = length(fit_wiener[["ranef"]][["coef"]])),
coef = gsub("phase", "", coef)
) %>%
filter(dpar == "drift") %>%
separate(coef, into = c("condition", "frequency"), sep = ":") %>%
merge(., df.fix %>% select(condition, frequency, dpar, fix_est)) %>%
# combine fixed estimate from the condition with the subject-specific estimate
mutate(
drift = Estimate + fix_est
)
Would this give me the subject-specific estimates for each condition? Or is there an easier way to do this, maybe an existing function?
- Operating System: Ubuntu 22.04 LTS
- brms Version: 2.22.0