Confusion about spread_draws

Load data

Code

library(brms)

Loading required package: Rcpp

Loading 'brms' package (version 2.20.15). Useful instructions
can be found by typing help('brms'). A more detailed introduction
to the package is available through vignette('brms_overview').


Attaching package: 'brms'

The following object is masked from 'package:stats':

    ar

Code

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.0     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2

── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Code

library(tidybayes)


Attaching package: 'tidybayes'

The following objects are masked from 'package:brms':

    dstudent_t, pstudent_t, qstudent_t, rstudent_t

Code

data_long_summarised = readRDS(file = file.path("data","osf_hedge_cwzds","data_long_summarised.Rds"))

data_long_summarised

# A tibble: 428 × 7
     pps session condition     n n_correct  n_go condition_contrast
   <dbl>   <int>     <int> <int>     <int> <dbl>              <dbl>
 1     1       1         0   450       450   449                0.5
 2     1       1         2   150        88    62               -0.5
 3     1       2         0   450       448   448                0.5
 4     1       2         2   150        59    91               -0.5
 5     2       1         0   450       450   449                0.5
 6     2       1         2   150       109    41               -0.5
 7     2       2         0   450       450   449                0.5
 8     2       2         2   150        87    63               -0.5
 9     3       1         0   450       450   449                0.5
10     3       1         2   150       118    32               -0.5
# ℹ 418 more rows

Fit model

Code

model_time_1 = brms::brm(
  n_go | trials(n) ~  1 + condition_contrast + (1 + condition_contrast | pps),
  # family = binomial(link = "probit"),
  family = binomial(link = "logit"),
  data = filter(data_long_summarised, session == 1 ),
  cores   = 2,
  chains  = 2,
  iter    = 10000,
  # init    = .01,
  adapt_delta = .95,
  backend = "cmdstanr",
  threads = threading(4),
  save_warmup = FALSE,
  refresh = 2000
)

Start sampling

Running MCMC with 2 parallel chains, with 4 thread(s) per chain...

Chain 1 Iteration:    1 / 10000 [  0%]  (Warmup) 
Chain 2 Iteration:    1 / 10000 [  0%]  (Warmup) 
Chain 2 Iteration: 2000 / 10000 [ 20%]  (Warmup) 
Chain 1 Iteration: 2000 / 10000 [ 20%]  (Warmup) 
Chain 1 Iteration: 4000 / 10000 [ 40%]  (Warmup) 
Chain 2 Iteration: 4000 / 10000 [ 40%]  (Warmup) 
Chain 1 Iteration: 5001 / 10000 [ 50%]  (Sampling) 
Chain 2 Iteration: 5001 / 10000 [ 50%]  (Sampling) 
Chain 2 Iteration: 7000 / 10000 [ 70%]  (Sampling) 
Chain 1 Iteration: 7000 / 10000 [ 70%]  (Sampling) 
Chain 2 Iteration: 9000 / 10000 [ 90%]  (Sampling) 
Chain 1 Iteration: 9000 / 10000 [ 90%]  (Sampling) 
Chain 2 Iteration: 10000 / 10000 [100%]  (Sampling) 
Chain 2 finished in 16.7 seconds.
Chain 1 Iteration: 10000 / 10000 [100%]  (Sampling) 
Chain 1 finished in 17.9 seconds.

Both chains finished successfully.
Mean chain execution time: 17.3 seconds.
Total execution time: 18.0 seconds.

Model Output

Code

model_time_1

 Family: binomial 
  Links: mu = logit 
Formula: n_go | trials(n) ~ 1 + condition_contrast + (1 + condition_contrast | pps) 
   Data: filter(data_long_summarised, session == 1) (Number of observations: 214) 
  Draws: 2 chains, each with iter = 10000; warmup = 5000; thin = 1;
         total post-warmup draws = 10000

Multilevel Hyperparameters:
~pps (Number of levels: 107) 
                                  Estimate Est.Error l-95% CI u-95% CI Rhat
sd(Intercept)                         0.58      0.05     0.49     0.68 1.00
sd(condition_contrast)                1.73      0.14     1.48     2.03 1.00
cor(Intercept,condition_contrast)     0.46      0.09     0.27     0.63 1.00
                                  Bulk_ESS Tail_ESS
sd(Intercept)                         2247     3730
sd(condition_contrast)                1231     2756
cor(Intercept,condition_contrast)      891     2334

Regression Coefficients:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept              1.76      0.07     1.63     1.89 1.00     2220     3740
condition_contrast     6.58      0.18     6.21     6.94 1.00     1517     2842

Draws were sampled using sample(hmc). 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).

Code

tidybayes::get_variables(model_time_1)

  [1] "b_Intercept"                           
  [2] "b_condition_contrast"                  
  [3] "sd_pps__Intercept"                     
  [4] "sd_pps__condition_contrast"            
  [5] "cor_pps__Intercept__condition_contrast"
  [6] "Intercept"                             
  [7] "r_pps[1,Intercept]"                    
  [8] "r_pps[2,Intercept]"                    
  [9] "r_pps[3,Intercept]"                    
 [10] "r_pps[4,Intercept]"                    
 [11] "r_pps[5,Intercept]"                    
 [12] "r_pps[7,Intercept]"                    
 [13] "r_pps[8,Intercept]"                    
 [14] "r_pps[9,Intercept]"                    
 [15] "r_pps[10,Intercept]"                   
 [16] "r_pps[11,Intercept]"                   
 [17] "r_pps[12,Intercept]"                   
 [18] "r_pps[13,Intercept]"                   
 [19] "r_pps[14,Intercept]"                   
 [20] "r_pps[15,Intercept]"                   
 [21] "r_pps[16,Intercept]"                   
 [22] "r_pps[18,Intercept]"                   
 [23] "r_pps[19,Intercept]"                   
 [24] "r_pps[20,Intercept]"                   
 [25] "r_pps[21,Intercept]"                   
 [26] "r_pps[22,Intercept]"                   
 [27] "r_pps[23,Intercept]"                   
 [28] "r_pps[24,Intercept]"                   
 [29] "r_pps[25,Intercept]"                   
 [30] "r_pps[26,Intercept]"                   
 [31] "r_pps[27,Intercept]"                   
 [32] "r_pps[28,Intercept]"                   
 [33] "r_pps[29,Intercept]"                   
 [34] "r_pps[30,Intercept]"                   
 [35] "r_pps[31,Intercept]"                   
 [36] "r_pps[32,Intercept]"                   
 [37] "r_pps[33,Intercept]"                   
 [38] "r_pps[34,Intercept]"                   
 [39] "r_pps[35,Intercept]"                   
 [40] "r_pps[36,Intercept]"                   
 [41] "r_pps[38,Intercept]"                   
 [42] "r_pps[39,Intercept]"                   
 [43] "r_pps[40,Intercept]"                   
 [44] "r_pps[41,Intercept]"                   
 [45] "r_pps[42,Intercept]"                   
 [46] "r_pps[43,Intercept]"                   
 [47] "r_pps[44,Intercept]"                   
 [48] "r_pps[45,Intercept]"                   
 [49] "r_pps[46,Intercept]"                   
 [50] "r_pps[47,Intercept]"                   
 [51] "r_pps[48,Intercept]"                   
 [52] "r_pps[49,Intercept]"                   
 [53] "r_pps[50,Intercept]"                   
 [54] "r_pps[1001,Intercept]"                 
 [55] "r_pps[1002,Intercept]"                 
 [56] "r_pps[1003,Intercept]"                 
 [57] "r_pps[1004,Intercept]"                 
 [58] "r_pps[1005,Intercept]"                 
 [59] "r_pps[1006,Intercept]"                 
 [60] "r_pps[1007,Intercept]"                 
 [61] "r_pps[1008,Intercept]"                 
 [62] "r_pps[1009,Intercept]"                 
 [63] "r_pps[1010,Intercept]"                 
 [64] "r_pps[1011,Intercept]"                 
 [65] "r_pps[1012,Intercept]"                 
 [66] "r_pps[1013,Intercept]"                 
 [67] "r_pps[1014,Intercept]"                 
 [68] "r_pps[1015,Intercept]"                 
 [69] "r_pps[1016,Intercept]"                 
 [70] "r_pps[1017,Intercept]"                 
 [71] "r_pps[1018,Intercept]"                 
 [72] "r_pps[1019,Intercept]"                 
 [73] "r_pps[1020,Intercept]"                 
 [74] "r_pps[1021,Intercept]"                 
 [75] "r_pps[1022,Intercept]"                 
 [76] "r_pps[1023,Intercept]"                 
 [77] "r_pps[1024,Intercept]"                 
 [78] "r_pps[1025,Intercept]"                 
 [79] "r_pps[1026,Intercept]"                 
 [80] "r_pps[1027,Intercept]"                 
 [81] "r_pps[1029,Intercept]"                 
 [82] "r_pps[1030,Intercept]"                 
 [83] "r_pps[1031,Intercept]"                 
 [84] "r_pps[1032,Intercept]"                 
 [85] "r_pps[1033,Intercept]"                 
 [86] "r_pps[1034,Intercept]"                 
 [87] "r_pps[1035,Intercept]"                 
 [88] "r_pps[1036,Intercept]"                 
 [89] "r_pps[1037,Intercept]"                 
 [90] "r_pps[1038,Intercept]"                 
 [91] "r_pps[1039,Intercept]"                 
 [92] "r_pps[1040,Intercept]"                 
 [93] "r_pps[1041,Intercept]"                 
 [94] "r_pps[1042,Intercept]"                 
 [95] "r_pps[1043,Intercept]"                 
 [96] "r_pps[1044,Intercept]"                 
 [97] "r_pps[1045,Intercept]"                 
 [98] "r_pps[1046,Intercept]"                 
 [99] "r_pps[1047,Intercept]"                 
[100] "r_pps[1048,Intercept]"                 
[101] "r_pps[1049,Intercept]"                 
[102] "r_pps[1050,Intercept]"                 
[103] "r_pps[1051,Intercept]"                 
[104] "r_pps[1052,Intercept]"                 
[105] "r_pps[1053,Intercept]"                 
[106] "r_pps[1054,Intercept]"                 
[107] "r_pps[1055,Intercept]"                 
[108] "r_pps[1057,Intercept]"                 
[109] "r_pps[1058,Intercept]"                 
[110] "r_pps[1059,Intercept]"                 
[111] "r_pps[1060,Intercept]"                 
[112] "r_pps[1061,Intercept]"                 
[113] "r_pps[1062,Intercept]"                 
[114] "r_pps[1,condition_contrast]"           
[115] "r_pps[2,condition_contrast]"           
[116] "r_pps[3,condition_contrast]"           
[117] "r_pps[4,condition_contrast]"           
[118] "r_pps[5,condition_contrast]"           
[119] "r_pps[7,condition_contrast]"           
[120] "r_pps[8,condition_contrast]"           
[121] "r_pps[9,condition_contrast]"           
[122] "r_pps[10,condition_contrast]"          
[123] "r_pps[11,condition_contrast]"          
[124] "r_pps[12,condition_contrast]"          
[125] "r_pps[13,condition_contrast]"          
[126] "r_pps[14,condition_contrast]"          
[127] "r_pps[15,condition_contrast]"          
[128] "r_pps[16,condition_contrast]"          
[129] "r_pps[18,condition_contrast]"          
[130] "r_pps[19,condition_contrast]"          
[131] "r_pps[20,condition_contrast]"          
[132] "r_pps[21,condition_contrast]"          
[133] "r_pps[22,condition_contrast]"          
[134] "r_pps[23,condition_contrast]"          
[135] "r_pps[24,condition_contrast]"          
[136] "r_pps[25,condition_contrast]"          
[137] "r_pps[26,condition_contrast]"          
[138] "r_pps[27,condition_contrast]"          
[139] "r_pps[28,condition_contrast]"          
[140] "r_pps[29,condition_contrast]"          
[141] "r_pps[30,condition_contrast]"          
[142] "r_pps[31,condition_contrast]"          
[143] "r_pps[32,condition_contrast]"          
[144] "r_pps[33,condition_contrast]"          
[145] "r_pps[34,condition_contrast]"          
[146] "r_pps[35,condition_contrast]"          
[147] "r_pps[36,condition_contrast]"          
[148] "r_pps[38,condition_contrast]"          
[149] "r_pps[39,condition_contrast]"          
[150] "r_pps[40,condition_contrast]"          
[151] "r_pps[41,condition_contrast]"          
[152] "r_pps[42,condition_contrast]"          
[153] "r_pps[43,condition_contrast]"          
[154] "r_pps[44,condition_contrast]"          
[155] "r_pps[45,condition_contrast]"          
[156] "r_pps[46,condition_contrast]"          
[157] "r_pps[47,condition_contrast]"          
[158] "r_pps[48,condition_contrast]"          
[159] "r_pps[49,condition_contrast]"          
[160] "r_pps[50,condition_contrast]"          
[161] "r_pps[1001,condition_contrast]"        
[162] "r_pps[1002,condition_contrast]"        
[163] "r_pps[1003,condition_contrast]"        
[164] "r_pps[1004,condition_contrast]"        
[165] "r_pps[1005,condition_contrast]"        
[166] "r_pps[1006,condition_contrast]"        
[167] "r_pps[1007,condition_contrast]"        
[168] "r_pps[1008,condition_contrast]"        
[169] "r_pps[1009,condition_contrast]"        
[170] "r_pps[1010,condition_contrast]"        
[171] "r_pps[1011,condition_contrast]"        
[172] "r_pps[1012,condition_contrast]"        
[173] "r_pps[1013,condition_contrast]"        
[174] "r_pps[1014,condition_contrast]"        
[175] "r_pps[1015,condition_contrast]"        
[176] "r_pps[1016,condition_contrast]"        
[177] "r_pps[1017,condition_contrast]"        
[178] "r_pps[1018,condition_contrast]"        
[179] "r_pps[1019,condition_contrast]"        
[180] "r_pps[1020,condition_contrast]"        
[181] "r_pps[1021,condition_contrast]"        
[182] "r_pps[1022,condition_contrast]"        
[183] "r_pps[1023,condition_contrast]"        
[184] "r_pps[1024,condition_contrast]"        
[185] "r_pps[1025,condition_contrast]"        
[186] "r_pps[1026,condition_contrast]"        
[187] "r_pps[1027,condition_contrast]"        
[188] "r_pps[1029,condition_contrast]"        
[189] "r_pps[1030,condition_contrast]"        
[190] "r_pps[1031,condition_contrast]"        
[191] "r_pps[1032,condition_contrast]"        
[192] "r_pps[1033,condition_contrast]"        
[193] "r_pps[1034,condition_contrast]"        
[194] "r_pps[1035,condition_contrast]"        
[195] "r_pps[1036,condition_contrast]"        
[196] "r_pps[1037,condition_contrast]"        
[197] "r_pps[1038,condition_contrast]"        
[198] "r_pps[1039,condition_contrast]"        
[199] "r_pps[1040,condition_contrast]"        
[200] "r_pps[1041,condition_contrast]"        
[201] "r_pps[1042,condition_contrast]"        
[202] "r_pps[1043,condition_contrast]"        
[203] "r_pps[1044,condition_contrast]"        
[204] "r_pps[1045,condition_contrast]"        
[205] "r_pps[1046,condition_contrast]"        
[206] "r_pps[1047,condition_contrast]"        
[207] "r_pps[1048,condition_contrast]"        
[208] "r_pps[1049,condition_contrast]"        
[209] "r_pps[1050,condition_contrast]"        
[210] "r_pps[1051,condition_contrast]"        
[211] "r_pps[1052,condition_contrast]"        
[212] "r_pps[1053,condition_contrast]"        
[213] "r_pps[1054,condition_contrast]"        
[214] "r_pps[1055,condition_contrast]"        
[215] "r_pps[1057,condition_contrast]"        
[216] "r_pps[1058,condition_contrast]"        
[217] "r_pps[1059,condition_contrast]"        
[218] "r_pps[1060,condition_contrast]"        
[219] "r_pps[1061,condition_contrast]"        
[220] "r_pps[1062,condition_contrast]"        
[221] "lprior"                                
[222] "lp__"                                  
[223] "accept_stat__"                         
[224] "treedepth__"                           
[225] "stepsize__"                            
[226] "divergent__"                           
[227] "n_leapfrog__"                          
[228] "energy__"

Extract draws using spread_draws

spread_draws to get posterior draws for participant 1 on the condition_contrast coefficient.

Weirdly, its giving me 20,000 draws for each participant. I don’t understand this as

Code

subj_draws_time_1 = model_time_1 %>%
  tidybayes::spread_draws(r_pps[pps,condition_contrast]) 

subj_draws_time_1$pps %>% table()

.
    1     2     3     4     5     7     8     9    10    11    12    13    14 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
   15    16    18    19    20    21    22    23    24    25    26    27    28 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
   29    30    31    32    33    34    35    36    38    39    40    41    42 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
   43    44    45    46    47    48    49    50  1001  1002  1003  1004  1005 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
 1006  1007  1008  1009  1010  1011  1012  1013  1014  1015  1016  1017  1018 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
 1019  1020  1021  1022  1023  1024  1025  1026  1027  1029  1030  1031  1032 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
 1033  1034  1035  1036  1037  1038  1039  1040  1041  1042  1043  1044  1045 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
 1046  1047  1048  1049  1050  1051  1052  1053  1054  1055  1057  1058  1059 
20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 20000 
 1060  1061  1062 
20000 20000 20000

Code

pps_1_draws = model_time_1 %>%
  tidybayes::spread_draws(r_pps[pps,condition_contrast]) %>%
  filter(pps == 1)

brms to get same posterior draws

Code

subj_draws_time_1_brms = as_draws_df(model_time_1) %>%
                    select(ends_with("condition_contrast]"))

Warning: Dropping 'draws_df' class as required metadata was removed.

Code

pps_1_draws_brms = subj_draws_time_1_brms$`r_pps[1,condition_contrast]`

Comparison

In both cases i want to analyse the posterior distribution for participant 1.

The first 10,000 draws from pps_1_draws (spread_draws output) match up to the brms output, but the second 10,000 draws look very different.

Code

hist(pps_1_draws_brms)

Code

hist(pps_1_draws$r_pps[1:10000])

Code

hist(pps_1_draws$r_pps[10001:20000])

Code

mean(pps_1_draws_brms)

[1] -0.5908408

Code

mean(pps_1_draws$r_pps[1:10000])

[1] -0.5908408

Code

mean(pps_1_draws$r_pps[10001:20000])

[1] 0.8137379

I’ve checked and its not the case that the first 10,000 draws are from chain 1:

Code

table(pps_1_draws$.chain[1:10000])


   1    2 
5000 5000

Anyway, I’m a little confused if there’s something i’ve done wrong here or if there’s been a strange interaction with brms here.

My setup:

Code

sessionInfo()

R version 4.3.2 (2023-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.3 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so;  LAPACK version 3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] tidybayes_3.0.6 lubridate_1.9.3 forcats_1.0.0   stringr_1.5.1  
 [5] dplyr_1.1.4     purrr_1.0.2     readr_2.1.5     tidyr_1.3.1    
 [9] tibble_3.2.1    ggplot2_3.5.0   tidyverse_2.0.0 brms_2.20.15   
[13] Rcpp_1.0.12    

loaded via a namespace (and not attached):
 [1] svUnit_1.0.6         tidyselect_1.2.0     loo_2.7.0           
 [4] fastmap_1.1.1        tensorA_0.36.2.1     shinystan_2.6.0     
 [7] promises_1.2.1       shinyjs_2.1.0        digest_0.6.34       
[10] timechange_0.3.0     mime_0.12            lifecycle_1.0.4     
[13] StanHeaders_2.32.6   ellipsis_0.3.2       magrittr_2.0.3      
[16] posterior_1.5.0      compiler_4.3.2       rlang_1.1.3         
[19] tools_4.3.2          igraph_2.0.2         utf8_1.2.4          
[22] yaml_2.3.8           knitr_1.45           bridgesampling_1.1-2
[25] htmlwidgets_1.6.4    pkgbuild_1.4.3       curl_5.2.1          
[28] plyr_1.8.9           dygraphs_1.1.1.6     abind_1.4-5         
[31] miniUI_0.1.1.1       withr_3.0.0          grid_4.3.2          
[34] stats4_4.3.2         fansi_1.0.6          xts_0.13.2          
[37] xtable_1.8-4         colorspace_2.1-0     inline_0.3.19       
[40] scales_1.3.0         gtools_3.9.5         cli_3.6.2           
[43] mvtnorm_1.2-4        rmarkdown_2.25       generics_0.1.3      
[46] RcppParallel_5.1.7   rstudioapi_0.15.0    tzdb_0.4.0          
[49] reshape2_1.4.4       rstan_2.32.6         shinythemes_1.2.0   
[52] bayesplot_1.11.1     parallel_4.3.2       matrixStats_1.2.0   
[55] base64enc_0.1-3      vctrs_0.6.5          V8_4.4.2            
[58] Matrix_1.6-5         jsonlite_1.8.8       arrayhelpers_1.1-0  
[61] hms_1.1.3            crosstalk_1.2.1      ggdist_3.3.1        
[64] glue_1.7.0           codetools_0.2-19     DT_0.32             
[67] distributional_0.4.0 stringi_1.8.3        gtable_0.3.4        
[70] later_1.3.2          QuickJSR_1.1.3       munsell_0.5.0       
[73] colourpicker_1.3.0   pillar_1.9.0         htmltools_0.5.7     
[76] Brobdingnag_1.2-9    R6_2.5.1             evaluate_0.23       
[79] shiny_1.8.0          lattice_0.21-9       markdown_1.12       
[82] backports_1.4.1      threejs_0.3.3        httpuv_1.6.14       
[85] rstantools_2.4.0     coda_0.19-4.1        gridExtra_2.3       
[88] nlme_3.1-163         checkmate_2.3.1      xfun_0.42           
[91] zoo_1.8-12           pkgconfig_2.0.3