Learning Curve new

library(tidyverse)

Warning: package 'tidyr' was built under R version 4.4.1

── 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.1     ✔ 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

library(broom.mixed)

Warning: package 'broom.mixed' was built under R version 4.4.1

library(devtools)

Loading required package: usethis

library(brms)

Loading required package: Rcpp
Loading 'brms' package (version 2.21.0). 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

library(mascutils)

Attaching package: 'mascutils'

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

    expand_grid

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

    mode

library(bayr)

Registered S3 methods overwritten by 'bayr':
  method             from     
  coef.brmsfit       brms     
  knit_print.tbl_obs mascutils
  predict.brmsfit    brms     
  print.tbl_obs      mascutils

Attaching package: 'bayr'

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

    as_tbl_obs, discard_all_na, discard_redundant, expand_grid,
    go_arrange, go_first, left_union, reorder_levels, rescale_centered,
    rescale_unit, rescale_zero_one, update_by, z_score, z_trans

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

    fixef, ranef

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

    expand_grid

library(readxl)

Warning: package 'readxl' was built under R version 4.4.1

library(ggplot2)
library(effects)

Loading required package: carData

Warning: package 'carData' was built under R version 4.4.1

lattice theme set by effectsTheme()
See ?effectsTheme for details.

library(lme4)

Warning: package 'lme4' was built under R version 4.4.1

Loading required package: Matrix

Warning: package 'Matrix' was built under R version 4.4.1

Attaching package: 'Matrix'

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

    expand, pack, unpack


Attaching package: 'lme4'

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

    fixef, ranef

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

    ngrps

library(haven)
library(lattice)

Warning: package 'lattice' was built under R version 4.4.1

library(car)

Warning: package 'car' was built under R version 4.4.1

Attaching package: 'car'

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

    logit

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

    recode

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

    some

library(knitr)

Warning: package 'knitr' was built under R version 4.4.1

library(reshape2)

Attaching package: 'reshape2'

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

    smiths

library(dplyr)
library(forcats)
library(DHARMa)

Warning: package 'DHARMa' was built under R version 4.4.1

This is DHARMa 0.4.7. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

library(Hmisc)

Warning: package 'Hmisc' was built under R version 4.4.1

Attaching package: 'Hmisc'

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

    src, summarize

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

    format.pval, units

library(phia)
library(lsmeans)

Loading required package: emmeans

Warning: package 'emmeans' was built under R version 4.4.1

Welcome to emmeans.
Caution: You lose important information if you filter this package's results.
See '? untidy'

Attaching package: 'emmeans'

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

    test

The 'lsmeans' package is now basically a front end for 'emmeans'.
Users are encouraged to switch the rest of the way.
See help('transition') for more information, including how to
convert old 'lsmeans' objects and scripts to work with 'emmeans'.

library(emmeans)
library(multcomp)

Warning: package 'multcomp' was built under R version 4.4.1

Loading required package: mvtnorm
Loading required package: survival

Attaching package: 'survival'

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

    kidney

Loading required package: TH.data

Warning: package 'TH.data' was built under R version 4.4.1

Loading required package: MASS

Warning: package 'MASS' was built under R version 4.4.1

Attaching package: 'MASS'

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

    select


Attaching package: 'TH.data'

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

    geyser

library(plotly)

Attaching package: 'plotly'

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

    select

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

    subplot

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

    last_plot

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

    filter

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

    layout

library(lmerTest)

Attaching package: 'lmerTest'

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

    lmer

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

    step

library(tinytex)
library(tidyverse)
library(devtools)
library(ggthemes)
library(janitor)

Attaching package: 'janitor'

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

    chisq.test, fisher.test

library(magrittr)

Attaching package: 'magrittr'

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

    set_names

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

    extract

library(tidyverse)
library(mascutils)
library(ggthemes)
library(afex)

Warning: package 'afex' was built under R version 4.4.1

************
Welcome to afex. For support visit: http://afex.singmann.science/
- Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
- Methods for calculating p-values with mixed(): 'S', 'KR', 'LRT', and 'PB'
- 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
- Get and set global package options with: afex_options()
- Set sum-to-zero contrasts globally: set_sum_contrasts()
- For example analyses see: browseVignettes("afex")
************

Attaching package: 'afex'

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

    lmer

setwd("/Users/can/Documents/Uni/Thesis/Data/E-Prime")

df <- read_csv("df.csv")

Rows: 31968 Columns: 14
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): procedure, feedback.CRESP, feedback.RESP, cue.OnsetTime, cue.OnsetD...
dbl (9): subject, session, sub.trial.number, feedback.ACC, feedback.RT, h, t...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

# Read the dataset
df <- read.csv("df.csv")

# Convert trial.RT from milliseconds to seconds
df <- df %>%
  mutate(trial.RTS = trial.RT / 1000)

# Convert trial.ACC into a dummy variable
df <- df %>%
  mutate(trial.acc = ifelse(trial.acc == 1, 1, 0))  # 1 if correct, 0 otherwise

# Aggregate data to have one row per trial
df_lc <- df %>%
  group_by(subject, session, trial) %>%
  summarise(
    trial.RT = mean(trial.RT, na.rm = TRUE),   # Mean RT per trial
    trial.RTS = mean(trial.RTS, na.rm = TRUE), # Mean RT in seconds
    trial.acc = mean(trial.acc, na.rm = TRUE)  # Mean accuracy per trial (should remain binary)
  ) %>%
  ungroup()

`summarise()` has grouped output by 'subject', 'session'. You can override
using the `.groups` argument.

# Adjust trial numbers so that they continue across sessions
df_lc <- df_lc %>%
  mutate(trial = case_when(
    session == 1 ~ trial,
    session == 2 ~ trial + 48,
    session == 3 ~ trial + 96
  ))

# Save the modified dataset
write.csv(df_lc, "df_lc.csv", row.names = FALSE)

# Load your dataset
df_lc1 <- read.csv("df_lc.csv")

#Raw data plots

df_new <- df_lc1 %>%
  ggplot(aes(x = trial, y = trial.RTS, colour = subject)) +
  geom_point() +
  geom_smooth(se = FALSE) +
  facet_wrap(~subject) +
  ylab("RT (s)") +
  xlab("Trial No.") + 
  scale_y_continuous()

plot(df_new)

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

df_lc1 %>% 
  filter(trial.RTS < 10) %>% 
  ggplot(aes(x = trial,
             y = trial.RTS)) +
  geom_point() +
  geom_smooth(se = F) +
  ylim(0,2.0) +
  facet_wrap(~session, scales = "free_y")

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Warning: Removed 10 rows containing non-finite outside the scale range
(`stat_smooth()`).

Warning: Removed 10 rows containing missing values or values outside the scale range
(`geom_point()`).

df_lc1 %>% 
  filter(trial.RTS < 10) %>% 
  ggplot(aes(x = trial,
             y = trial.RTS,
             group = subject)) +
  facet_grid(~session) +
  geom_smooth(se = F)

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

df_lc1 %>% 
  filter(trial.RTS < 10) %>% 
  ggplot(aes(x = trial,
             y = trial.RTS,
             group = subject)) +
  geom_smooth(se = F)

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

#Classic Learning curve model all 3 sessions

# ampl = possible amount of learning before diminishing returns
# rate = how fast a person learns
# asymptote = maximal learning (lowest RT)

F_ary <- formula(trial.RTS ~ asym + ampl * exp(-rate * trial))

# Base model only accounts for subject-level model
F_ary_ef_1 <- list(formula(ampl ~ 1|subject),
                   formula(rate ~ 1|subject),
                   formula(asym ~ 1|subject))

# Priors (Weakly informative)
F_ary_prior <- c(set_prior("normal(5, 100)", nlpar = "ampl", lb = 0),
                 set_prior("normal(.5, 3)", nlpar = "rate", lb = 0),
                 set_prior("normal(3, 20)", nlpar = "asym", lb = 0))

# Modeling at the subject level   
M_1 <- 
  df_lc1 %>% 
  brm(bf(F_ary,
         flist = F_ary_ef_1,
         nl = T), 
      prior = F_ary_prior,
      family = "exgaussian",
      data = .,
      iter = 5000, #2000 - #5000
      warmup = 4000, #1000 - #4000
      save_pars = save_pars(all=TRUE))

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.6)’
using SDK: ‘’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.0005 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 1: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1681.92 seconds (Warm-up)
Chain 1:                281.4 seconds (Sampling)
Chain 1:                1963.32 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000298 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 2.98 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 2: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1118.65 seconds (Warm-up)
Chain 2:                285.292 seconds (Sampling)
Chain 2:                1403.94 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 3:   Stan can't start sampling from this initial value.
Chain 3: 
Chain 3: Gradient evaluation took 0.000295 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 2.95 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 3: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1061.67 seconds (Warm-up)
Chain 3:                296.707 seconds (Sampling)
Chain 3:                1358.38 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4:   Error evaluating the log probability at the initial value.
Chain 4: Exception: exp_mod_normal_lpdf: Location parameter[1193] is -inf, but must be finite! (in 'anon_model', line 104, column 4 to column 67)
Chain 4: Rejecting initial value:
Chain 4:   Error evaluating the log probability at the initial value.
Chain 4: Exception: exp_mod_normal_lpdf: Location parameter[110] is -inf, but must be finite! (in 'anon_model', line 104, column 4 to column 67)
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4:   Error evaluating the log probability at the initial value.
Chain 4: Exception: exp_mod_normal_lpdf: Location parameter[60] is inf, but must be finite! (in 'anon_model', line 104, column 4 to column 67)
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: 
Chain 4: Gradient evaluation took 0.000299 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2.99 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 4: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1119.07 seconds (Warm-up)
Chain 4:                297.726 seconds (Sampling)
Chain 4:                1416.79 seconds (Total)
Chain 4: 

Warning: There were 67 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: There were 3839 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: The largest R-hat is 2.05, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

P_1 <- posterior(M_1) 
PP_1 <- post_pred(M_1, thin = 10)

save(M_1, P_1, PP_1, df_lc1, file = "2025_analyses_long.Rda")

fixef(M_1)

                Estimate  Est.Error       Q2.5     Q97.5
ampl_Intercept 0.6918332 0.35764996 0.16993574 1.4471944
rate_Intercept 0.2348865 0.10112664 0.09865021 0.4816092
asym_Intercept 0.5727385 0.06574976 0.46036416 0.7414013

grpef(P_1)

Coefficient estimates with 95% credibility limits

nonlin	center	lower	upper
ampl	1.1942577	0.7252313	1.9853807
rate	0.2383765	0.1144898	0.4396119
asym	0.2559751	0.1533502	0.4382487

ranef(M_1)

$subject
, , ampl_Intercept

      Estimate Est.Error       Q2.5        Q97.5
2   0.11753494 0.3576811 -0.6400706  0.651299532
3  -0.01488935 0.3650256 -0.7808086  0.555309122
4  -0.96440337 0.3531747 -1.7051928 -0.431180421
5   0.72544185 1.1637706 -1.4039795  2.462198441
7  -0.66918376 0.4017671 -1.5264499 -0.089457400
8   0.11154595 0.3690096 -0.6784986  0.679738721
10 -0.41045452 0.3894940 -1.2533689  0.297294252
11  0.15812291 0.5567299 -0.7787330  1.429960394
12  0.41720795 1.2322087 -2.1092242  2.781110012
13 -0.32215999 0.3620590 -1.0844900  0.219906677
14 -0.85207933 0.3608505 -1.6258601 -0.314372191
15 -0.22064036 0.3633444 -0.9938895  0.338733392
16 -0.14691836 0.6402021 -1.1933635  1.039570537
17 -0.14495024 0.3698142 -0.9422539  0.447618705
18  0.64152078 0.4518566 -0.3162152  1.466171010
19 -0.61460813 1.4096388 -2.0427359  3.085060266
20  3.73640285 1.0073064  1.9786250  5.735368799
22  0.37476807 0.4018551 -0.4653025  1.090931275
23 -0.88855732 0.4967742 -1.7973498 -0.003986641

, , rate_Intercept

       Estimate  Est.Error       Q2.5        Q97.5
2  -0.193357429 0.10113490 -0.4385640 -0.056248444
3  -0.033473630 0.10574499 -0.2781545  0.134145787
4  -0.220900250 0.10054567 -0.4691718 -0.082866216
5   0.255734833 0.27159609 -0.2351429  0.666735217
7   0.124596643 0.22588369 -0.2112108  0.646977584
8  -0.056187806 0.10071196 -0.2909616  0.088611792
10  0.087746156 0.23928727 -0.2674117  0.615593808
11  0.277966456 0.26812045 -0.1973774  0.815130786
12  0.001678211 0.23944265 -0.3520494  0.599103114
13 -0.142229458 0.10272207 -0.3870061  0.004651431
14 -0.176820115 0.10654255 -0.4351661 -0.012994103
15 -0.058061058 0.11975892 -0.3092521  0.158074946
16  0.072562176 0.23826970 -0.2420494  0.530076125
17 -0.135653424 0.10325893 -0.3754811  0.020376706
18  0.117917981 0.12144793 -0.1559355  0.339752804
19 -0.029083110 0.30958253 -0.3708509  0.730946588
20  0.294755184 0.10265997  0.0729073  0.487625460
22  0.104756087 0.12578310 -0.1500827  0.362138340
23 -0.229938106 0.09897886 -0.4737896 -0.096050631

, , asym_Intercept

      Estimate  Est.Error         Q2.5       Q97.5
2  -0.06501642 0.06743103 -0.230709468  0.05075396
3  -0.08504995 0.06694851 -0.256506878  0.02990957
4   0.18736787 0.09985424  0.021554908  0.42455384
5   0.09557825 0.20185543 -0.109835499  0.63637921
7  -0.13148964 0.06625120 -0.298385193 -0.01852103
8  -0.15321832 0.06674682 -0.327471343 -0.03618161
10 -0.30512767 0.10934491 -0.511301523 -0.04327792
11 -0.24627612 0.09210245 -0.465575779 -0.10906591
12 -0.23754770 0.12175086 -0.535939330 -0.07559719
13 -0.13619972 0.06636590 -0.305034274 -0.02142851
14 -0.09909659 0.06728726 -0.261013124  0.02055028
15  0.05303346 0.06686933 -0.118578932  0.16502826
16  0.03717742 0.12503546 -0.140797844  0.36323165
17 -0.05514775 0.06640160 -0.218951557  0.05549084
18  0.15314047 0.06675378 -0.017199574  0.26711778
19  0.45508251 0.26420711  0.008309245  0.93730757
20  0.51536706 0.06979060  0.334383839  0.63213894
22  0.21826458 0.06741601  0.043325316  0.33366748
23 -0.10722189 0.25591027 -0.631161133  0.39472699

# Extract random effects for the "subject" grouping factor
ranef_df <- ranef(M_1)$subject %>%  
  as.data.frame() %>%
  rownames_to_column(var = "subject")  # Convert row names to a column

# Rename columns for clarity (adjust based on `ranef(M_1)` output)
colnames(ranef_df) <- c("subject", "term", "Estimate", "Q2.5", "Q97.5")

# Convert subject to a factor for clear axis labels
ranef_df$subject <- as.factor(ranef_df$subject)

# Create a single plot with all subjects
ggplot(ranef_df, aes(x = subject,  
                      y = Estimate,  
                      ymin = Q2.5,  
                      ymax = Q97.5)) +  
  geom_crossbar(width = 0.4, color = "blue", fill = "lightblue") +  
  geom_point(size = 3, color = "black") +  # Highlight mean estimates
  theme_minimal() +
  labs(title = "Subject-Level Random Effects",
       x = "Participant",
       y = "Random Effect Estimate") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))  # Rotate x-axis labels

# Generate fitted values (predicted learning curves)
fitted_df <- fitted(M_1, summary = TRUE) %>%  # Extract fitted responses
  as.data.frame()

# Ensure the dataset has subject & trial numbers
fitted_df <- df_lc %>% 
  select(subject, trial) %>%  # Keep only subject & trial
  bind_cols(fitted_df)  # Merge with fitted values

# Plot estimated learning curves
ggplot(fitted_df, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curves
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # 95% Credible Interval
  labs(title = "Estimated Learning Curves",
       x = "Trial Number",
       y = "Predicted RT (s)",
       color = "Subject") +
  theme_minimal()

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

ggplot() +
  geom_line(data = df_lc, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +  
  geom_line(data = fitted_df, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_df, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  labs(title = "Observed vs. Estimated Learning Curves",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curves for each subject in separate panels
ggplot(fitted_df, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curve
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # Confidence band
  facet_wrap(~subject, scales = "free_y") +  # Separate plot per subject
  labs(title = "Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Predicted RT (s)") +
  theme_minimal() +
  theme(legend.position = "none")  

# Plot observed vs. estimated learning curves for each subject
ggplot() +
  geom_line(data = df_lc, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +
  geom_line(data = fitted_df, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_df, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  facet_wrap(~subject, scales = "free_y") +
  labs(title = "Observed vs. Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Reaction Time (s)",
       color = "Subject") +
  theme_minimal() +
  theme(legend.position = "none")

# Ensure fitted values contain trial information
fitted_df <- fitted(M_1, summary = TRUE) %>%
  as.data.frame() %>%
  mutate(trial = df_lc$trial)  # Explicitly add trial column from df_lc

# Aggregate model predictions at the group level
group_fitted <- fitted_df %>%
  group_by(trial) %>%
  summarise(
    mean_Estimate = mean(Estimate, na.rm = TRUE),
    lower_CI = mean(Q2.5, na.rm = TRUE),
    upper_CI = mean(Q97.5, na.rm = TRUE)
  )

# Aggregate observed data at the group level
group_obs <- df_lc %>%
  group_by(trial) %>%
  summarise(
    mean_RT = mean(trial.RTS, na.rm = TRUE),
    se_RT = sd(trial.RTS, na.rm = TRUE) / sqrt(n())
  )




# Plot group-level learning curve
ggplot() +
  # Observed mean RT with standard error
  geom_line(data = group_obs, aes(x = trial, y = mean_RT), color = "black", size = 1, alpha = 0.5) +
  geom_ribbon(data = group_obs, aes(x = trial, ymin = mean_RT - se_RT, ymax = mean_RT + se_RT), alpha = 0.2, fill = "grey") +

  # Estimated mean RT from model
  geom_line(data = group_fitted, aes(x = trial, y = mean_Estimate), color = "blue", size = 1) +
  geom_ribbon(data = group_fitted, aes(x = trial, ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +

  # Labels and theme
  labs(title = "Group-Level Learning Curve: Observed vs. Estimated RT",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curve (group-level)
ggplot(group_fitted, aes(x = trial, y = mean_Estimate)) +
  geom_line(color = "blue", size = 1) +  # Estimated RT line
  geom_ribbon(aes(ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +  # Confidence interval
  labs(title = "Group-Level Estimated Learning Curve",
       x = "Trial Number",
       y = "Predicted Reaction Time (s)") +
  theme_minimal()

# Generate fitted values (predicted learning curves)
fitted_df <- fitted(M_1, summary = TRUE) %>%
  as.data.frame()

# Merge fitted values with the dataset
fitted_df <- df_lc %>%
  select(subject, trial, session) %>%  # Keep subject, trial, and session
  bind_cols(fitted_df)

# Filter trials to match session-specific trial numbering
group_fitted <- fitted_df %>%
  group_by(session, trial) %>%
  summarise(
    mean_Estimate = mean(Estimate, na.rm = TRUE),
    lower_CI = mean(Q2.5, na.rm = TRUE),
    upper_CI = mean(Q97.5, na.rm = TRUE)
  ) %>%
  filter((session == 1 & trial >= 1 & trial <= 48) |  # Only keep valid trials per session
         (session == 2 & trial >= 49 & trial <= 96) |
         (session == 3 & trial >= 97 & trial <= 144))

`summarise()` has grouped output by 'session'. You can override using the
`.groups` argument.

# Plot estimated learning curves separately for each session with correct x-axis
ggplot(group_fitted, aes(x = trial, y = mean_Estimate, color = as.factor(session))) +
  geom_line(size = 1) +  # Estimated RT line per session
  geom_ribbon(aes(ymin = lower_CI, ymax = upper_CI, fill = as.factor(session)), alpha = 0.2) +  # Confidence interval
  facet_wrap(~session, scales = "free_x") +  # Adjust x-axis to match session trials
  scale_x_continuous(breaks = seq(1, 144, by = 5)) +  # Show relevant trial numbers
  labs(title = "Session-Specific Estimated Learning Curves",
       x = "Trial Number",
       y = "Predicted Reaction Time (s)",
       color = "Session",
       fill = "Session") +
  theme_minimal()

#Models for each session individually

# Function to fit models for each session
fit_model_for_session <- function(data, session_number) {
  brm(
    bf(F_ary, flist = F_ary_ef_1, nl = TRUE),
    prior = F_ary_prior,
    family = "exgaussian",
    data = data,
    iter = 5000,
    warmup = 4000,
    save_pars = save_pars(all = TRUE)
  )
}

# Fit models for each session
M_session1 <- fit_model_for_session(df_lc1 %>% filter(session == 1), 1)

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.6)’
using SDK: ‘’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000258 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.58 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 1: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 309.628 seconds (Warm-up)
Chain 1:                94.619 seconds (Sampling)
Chain 1:                404.247 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 2:   Stan can't start sampling from this initial value.
Chain 2: Rejecting initial value:
Chain 2:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 2:   Stan can't start sampling from this initial value.
Chain 2: 
Chain 2: Gradient evaluation took 0.000106 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.06 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 2: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 313.718 seconds (Warm-up)
Chain 2:                80.01 seconds (Sampling)
Chain 2:                393.728 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 3:   Stan can't start sampling from this initial value.
Chain 3: 
Chain 3: Gradient evaluation took 0.000101 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.01 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 3: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 232.596 seconds (Warm-up)
Chain 3:                71.334 seconds (Sampling)
Chain 3:                303.93 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000101 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.01 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 4: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 296.423 seconds (Warm-up)
Chain 4:                94.563 seconds (Sampling)
Chain 4:                390.986 seconds (Total)
Chain 4: 

Warning: There were 261 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: There were 2645 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: The largest R-hat is 1.57, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

M_session2 <- fit_model_for_session(df_lc1 %>% filter(session == 2), 2)

Compiling Stan program...
Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.6)’
using SDK: ‘’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.000234 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.34 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 1: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 107.664 seconds (Warm-up)
Chain 1:                13.117 seconds (Sampling)
Chain 1:                120.781 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 2:   Stan can't start sampling from this initial value.
Chain 2: 
Chain 2: Gradient evaluation took 0.000103 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.03 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 2: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 107.2 seconds (Warm-up)
Chain 2:                12.991 seconds (Sampling)
Chain 2:                120.191 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: Rejecting initial value:
Chain 3:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 3:   Stan can't start sampling from this initial value.
Chain 3: Rejecting initial value:
Chain 3:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 3:   Stan can't start sampling from this initial value.
Chain 3: Rejecting initial value:
Chain 3:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 3:   Stan can't start sampling from this initial value.
Chain 3: 
Chain 3: Gradient evaluation took 0.000101 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.01 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 3: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 101.511 seconds (Warm-up)
Chain 3:                12.793 seconds (Sampling)
Chain 3:                114.304 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000102 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.02 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 4: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 87.909 seconds (Warm-up)
Chain 4:                12.785 seconds (Sampling)
Chain 4:                100.694 seconds (Total)
Chain 4: 

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

M_session3 <- fit_model_for_session(df_lc1 %>% filter(session == 3), 3)

Compiling Stan program...
Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.6)’
using SDK: ‘’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Gradient evaluation took 0.000255 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2.55 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 1: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 196.369 seconds (Warm-up)
Chain 1:                15.127 seconds (Sampling)
Chain 1:                211.496 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000106 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1.06 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 2: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 103.198 seconds (Warm-up)
Chain 2:                21.531 seconds (Sampling)
Chain 2:                124.729 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000108 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1.08 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 3: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 79.381 seconds (Warm-up)
Chain 3:                105.768 seconds (Sampling)
Chain 3:                185.149 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: Rejecting initial value:
Chain 4:   Error evaluating the log probability at the initial value.
Chain 4: Exception: exp_mod_normal_lpdf: Location parameter[236] is -inf, but must be finite! (in 'anon_model', line 104, column 4 to column 67)
Chain 4: Rejecting initial value:
Chain 4:   Error evaluating the log probability at the initial value.
Chain 4: Exception: exp_mod_normal_lpdf: Location parameter[145] is inf, but must be finite! (in 'anon_model', line 104, column 4 to column 67)
Chain 4: 
Chain 4: Gradient evaluation took 0.000107 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 1.07 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 4: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 78.619 seconds (Warm-up)
Chain 4:                13.8 seconds (Sampling)
Chain 4:                92.419 seconds (Total)
Chain 4: 

Warning: There were 6 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: There were 997 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

# Save models
save(M_session1, M_session2, M_session3, file = "tryout1")

PS_1 <- posterior(M_session1) 
PS_2 <- posterior(M_session2)
PS_3 <- posterior(M_session3)

PPS_1 <- post_pred(M_session1, thin = 10)
PPS_2 <- post_pred(M_session2, thin = 10)
PPS_3 <- post_pred(M_session3, thin = 10)

save(M_session1, PS_1, PPS_1, df_lc1, file = "2025_analyses_long_session1.Rda")
save(M_session2, PS_2, PPS_2, df_lc1, file = "2025_analyses_long_session1.Rda")
save(M_session3, PS_3, PPS_3, df_lc1, file = "2025_analyses_long_session1.Rda")

fixef(M_session1)

                Estimate  Est.Error      Q2.5     Q97.5
ampl_Intercept 1.0583820 0.45605673 0.3800999 2.0140477
rate_Intercept 0.3678676 0.13913260 0.1344763 0.6424165
asym_Intercept 0.4640709 0.04593364 0.3713876 0.5535515

fixef(M_session2)

                  Estimate   Est.Error       Q2.5       Q97.5
ampl_Intercept 145.7146485 60.07387043 51.2020254 282.0981433
rate_Intercept   0.1470722  0.01268626  0.1223888   0.1726177
asym_Intercept   0.5267272  0.05720103  0.4102401   0.6336388

fixef(M_session3)

                  Estimate   Est.Error       Q2.5        Q97.5
ampl_Intercept 97.17413926 61.66895053 8.27460397 239.95354800
rate_Intercept  0.07356789  0.01073071 0.04664041   0.09207138
asym_Intercept  0.53267696  0.05618905 0.42144579   0.64917737

grpef(PS_1)

Coefficient estimates with 95% credibility limits

nonlin	center	lower	upper
ampl	1.2736745	0.2937075	2.2058547
rate	0.2897463	0.1084438	0.5860325
asym	0.1772758	0.1115694	0.2858992

grpef(PS_2)

Coefficient estimates with 95% credibility limits

nonlin	center	lower	upper
ampl	1.9131909	0.0916077	11.8147077
rate	0.0255370	0.0150346	0.0504462
asym	0.2456421	0.1791863	0.3685418

grpef(PS_3)

Coefficient estimates with 95% credibility limits

nonlin	center	lower	upper
ampl	2.0751380	0.0906450	11.0793197
rate	0.0121615	0.0058752	0.0279838
asym	0.2184276	0.1598950	0.3209582

ranef(M_session1)

$subject
, , ampl_Intercept

      Estimate Est.Error       Q2.5       Q97.5
2  -0.09138094 0.4523406 -1.0906605  0.61649385
3  -0.32450985 0.4622850 -1.3009494  0.36811885
4  -0.80324212 0.5469425 -1.8357259  0.30111083
5   0.94741595 0.6287180 -0.1334461  2.27515203
7  -0.99992995 0.5038288 -2.0624463 -0.23617447
8  -0.24815417 0.4547141 -1.2079579  0.43796722
10 -0.60785055 0.4784022 -1.6137250  0.14871284
11  0.06377522 0.5477800 -1.0965059  1.14759789
13 -0.62488172 0.4541604 -1.5897053  0.06794245
14 -1.05243932 0.4871263 -2.1284769 -0.31951950
15 -0.55148576 0.4626417 -1.5103273  0.15237535
16  0.15002782 0.4560756 -0.8728660  0.98310422
17  0.56468406 0.9671868 -0.9450025  2.45803611
18  0.15819204 0.4874791 -0.8026650  1.10512714
19  1.25907453 1.2135939 -0.3275939  4.22557799
20  3.08704613 2.0279054 -0.3871582  5.68075779
22 -0.10727256 0.4599813 -1.0710418  0.62294023
23 -0.88795963 0.4857751 -1.9195896 -0.11791136

, , rate_Intercept

      Estimate Est.Error        Q2.5         Q97.5
2  -0.34373623 0.1401350 -0.61832406 -0.1069713572
3  -0.26526592 0.1405848 -0.54694988 -0.0267051934
4   0.28554355 0.3581923 -0.24318106  1.1477085964
5   0.18884290 0.1380480 -0.08841152  0.4774395982
7   0.01089929 0.3304151 -0.48570104  0.8322756401
8  -0.23947817 0.1405256 -0.52153479 -0.0046291056
10  0.05643699 0.2556618 -0.37591004  0.6670312493
11  0.26422210 0.1955397 -0.07719000  0.7109405906
13 -0.31782025 0.1420963 -0.59779504 -0.0790247517
14  0.13480117 0.3143118 -0.37295754  0.8300818225
15 -0.25684856 0.1461837 -0.55339157 -0.0006197833
16  0.08922457 0.1478653 -0.21190759  0.3775766875
17  0.06576447 0.2271024 -0.34351562  0.4869924450
18 -0.06766627 0.1453418 -0.37119687  0.1979454024
19  0.30169813 0.1992253 -0.05347284  0.7150831178
20  0.04640372 0.1888312 -0.30508499  0.3416677760
22 -0.16453999 0.1435160 -0.45186121  0.0737285584
23  0.21343332 0.3083069 -0.30421009  0.9540496809

, , asym_Intercept

       Estimate  Est.Error        Q2.5        Q97.5
2  -0.182615024 0.13789816 -0.46758286  0.052138663
3  -0.100151367 0.05540527 -0.20934927  0.003732132
4   0.076818477 0.05495173 -0.03087858  0.177851936
5   0.003269459 0.04885359 -0.09411154  0.099059212
7  -0.036013046 0.06974700 -0.17461442  0.087960657
8  -0.103571231 0.05221548 -0.20710626 -0.002522979
10 -0.216290499 0.04954845 -0.31658868 -0.119486057
11 -0.141316870 0.04793219 -0.23519764 -0.045296072
13 -0.122015798 0.08750474 -0.33699194  0.017101655
14 -0.044959269 0.05976750 -0.15433243  0.073934458
15  0.081781249 0.06143686 -0.05454576  0.192610931
16 -0.003761176 0.04846442 -0.10033641  0.093825940
17  0.107079678 0.05221210  0.00100723  0.208418998
18  0.231703014 0.04888827  0.13659281  0.329882017
19  0.069667806 0.04808187 -0.02667731  0.164888715
20  0.343022913 0.22573222 -0.17711251  0.586168854
22  0.224279751 0.04948635  0.12849787  0.326691130
23 -0.160105602 0.05111858 -0.25785316 -0.061857310

ranef(M_session2)

$subject
, , ampl_Intercept

       Estimate Est.Error       Q2.5     Q97.5
2   0.530384679  5.764456  -7.227947 10.155942
3   0.393310509  4.824213  -6.957956  9.265926
4   0.053567010  5.649642  -8.534081  9.249861
5   0.320245595  5.564162  -7.971223  9.806939
7   0.017780906  4.665565  -8.406713  8.693344
8   0.042257390  4.874238  -8.576835  8.941472
10 -0.302390423  5.032245  -9.674653  7.542486
11 -0.021419616  4.571853  -8.605071  8.610228
13  0.002802078  4.864545  -8.732492  8.711812
14 -0.021880917  6.655613  -8.227736  8.019480
15 -0.198997793  5.006877 -10.074703  7.969929
16  0.289931868  5.227323  -7.875511  9.506049
17  0.102662736  5.268960  -8.308518  9.251516
18  0.055346012  4.804475  -8.195502  8.813200
19 -0.306985159  5.469714  -9.388313  7.521462
20 -0.268189063  7.418496  -9.192686  9.185767
22 -0.168857205  5.532145  -8.718267  9.060860
23 -0.028053717  4.723481  -8.550613  8.906066

, , rate_Intercept

       Estimate   Est.Error          Q2.5        Q97.5
2  -0.029791761 0.010291866 -0.0520253142 -0.011940585
3  -0.030965757 0.010070163 -0.0539645757 -0.014266168
4  -0.003816855 0.014009323 -0.0287034630  0.025397854
5  -0.036711803 0.009724289 -0.0589698532 -0.021121868
7   0.004332453 0.016792904 -0.0224683125  0.043298715
8  -0.016763618 0.010772982 -0.0395976566  0.001925457
10  0.028633626 0.019833235 -0.0006260964  0.077008897
11 -0.005555724 0.013936328 -0.0296003763  0.024949978
13  0.002840346 0.017438888 -0.0242036805  0.046985286
14  0.025596526 0.018177522 -0.0021071998  0.071210011
15  0.015685594 0.017819459 -0.0130368642  0.059156199
16 -0.014788943 0.012079500 -0.0391146501  0.009647957
17 -0.021116705 0.010560194 -0.0442982061 -0.002616781
18  0.009274116 0.018973321 -0.0191403926  0.054933822
19  0.019387546 0.020133352 -0.0102962820  0.067584107
20  0.033047824 0.019833061  0.0024172538  0.077214895
22  0.025474595 0.020229764 -0.0046436330  0.073324180
23 -0.008053617 0.012607137 -0.0315414440  0.018043923

, , asym_Intercept

      Estimate  Est.Error         Q2.5        Q97.5
2  -0.09011823 0.06039019 -0.207012466  0.031540889
3  -0.02970794 0.05979984 -0.141404555  0.090122602
4   0.07290794 0.05918395 -0.039797706  0.193863000
5  -0.02766204 0.06002121 -0.140812890  0.092334524
7  -0.12028295 0.05850930 -0.229871231 -0.001914643
8  -0.17895952 0.05955478 -0.291709496 -0.056196495
10 -0.30450529 0.05873606 -0.414361247 -0.184790180
11 -0.23429030 0.05843591 -0.345303690 -0.114758257
13 -0.13320618 0.05894598 -0.246179448 -0.012759201
14 -0.05554354 0.05808855 -0.167520370  0.059346554
15  0.11687538 0.05941820  0.005862119  0.236635939
16 -0.01989221 0.05909908 -0.131802964  0.099822282
17 -0.04659973 0.05897489 -0.160588118  0.074956999
18  0.12577866 0.05910122  0.014175863  0.246146651
19  0.11065520 0.05887236 -0.002517795  0.228699555
20  0.72037428 0.06065569  0.606451214  0.842854766
22  0.30709430 0.05992261  0.193875090  0.429436817
23 -0.20605059 0.05845149 -0.316942669 -0.087284300

ranef(M_session3)

$subject
, , ampl_Intercept

      Estimate Est.Error       Q2.5     Q97.5
2   0.91894546  6.446026  -6.210488 12.078280
3  -0.10658768  5.046939  -8.427917  8.390724
4   0.01489996  5.450941  -8.506889  8.541094
5   0.45108907  4.876487  -7.398015  9.985549
7  -0.16240583  4.697411  -8.480970  7.780661
8   0.36606654  4.434718  -7.247009  9.371419
10  0.31412061  4.245442  -7.228767  9.325259
11  0.21204197  4.986138  -7.893876  9.075655
12  0.76396071  5.161050  -6.786306 11.584688
13 -0.50869887  6.102272 -10.515500  8.126229
14 -0.05183704  4.431837  -8.670466  8.345640
15 -0.10821227  5.933654  -8.498698  8.534354
16 -0.18299523  3.923670  -8.621261  7.438441
17  0.34007028  4.344750  -7.649562  9.271899
18 -0.78245074  6.240655 -12.362325  7.523469
19 -0.12148027  4.936748  -8.543388  8.853207
20  0.11834616  4.617213  -7.730508  7.782869
22 -0.84313787  6.968981 -12.033779  6.972541
23 -0.12734026  4.550771  -9.320858  8.544892

, , rate_Intercept

        Estimate   Est.Error         Q2.5         Q97.5
2  -0.0185741911 0.006027828 -0.032340639 -0.0084040340
3   0.0052191559 0.011160067 -0.011239619  0.0308382788
4  -0.0005077607 0.009669861 -0.016265493  0.0216831603
5  -0.0121864103 0.006846709 -0.026707699  0.0004761656
7   0.0057668374 0.010503360 -0.010247932  0.0306843111
8  -0.0066474599 0.008693531 -0.022263216  0.0115302083
10 -0.0079740326 0.007247817 -0.022858970  0.0057911137
11 -0.0026993036 0.010030887 -0.018693479  0.0196133704
12 -0.0152005928 0.006698292 -0.029031492 -0.0026017764
13  0.0117643207 0.010490225 -0.004291721  0.0367086307
14  0.0041271880 0.010084529 -0.011516206  0.0279643114
15  0.0063707890 0.011451092 -0.010415308  0.0340006883
16  0.0048691766 0.010029413 -0.010635790  0.0286684784
17 -0.0095788276 0.006761604 -0.024301207  0.0023229049
18  0.0138953593 0.010638814 -0.002577255  0.0396680974
19  0.0069521340 0.010592571 -0.009194587  0.0318986367
20 -0.0032804218 0.009746202 -0.019663232  0.0203233877
22  0.0140012377 0.010645022 -0.002389689  0.0391306365
23  0.0048623263 0.010342511 -0.011892386  0.0290826858

, , asym_Intercept

      Estimate  Est.Error         Q2.5         Q97.5
2  -0.09349963 0.06401007 -0.227143818  0.0232925018
3  -0.04256937 0.05862423 -0.164120541  0.0729393201
4   0.11154060 0.05943559 -0.007477662  0.2266290476
5   0.08431131 0.06151942 -0.040972128  0.2048795065
7  -0.12632079 0.05831472 -0.242577691 -0.0101739669
8  -0.11078024 0.06052842 -0.234213291  0.0052229352
10 -0.25678950 0.06115490 -0.377104389 -0.1378693348
11 -0.12713146 0.06017918 -0.242734378 -0.0112672790
12 -0.26694318 0.06252497 -0.392276750 -0.1458675600
13 -0.12498908 0.05818583 -0.244551062 -0.0110333120
14 -0.13202029 0.05800613 -0.248364425 -0.0202317748
15  0.02770584 0.05859092 -0.089453486  0.1396342292
16  0.01104030 0.05792008 -0.102155753  0.1208195507
17 -0.11642547 0.05992841 -0.236677647 -0.0002791129
18  0.34067147 0.05926350  0.220890745  0.4582545348
19  0.25832448 0.05931375  0.141996648  0.3739073975
20  0.41965038 0.06267956  0.297269513  0.5440867738
22  0.35445877 0.05830648  0.237805273  0.4730233012
23 -0.16845180 0.05829696 -0.285954026 -0.0557859987

# Extract full fixed-effect summaries from each session model
fixef_s1 <- as.data.frame(fixef(M_session1))
fixef_s2 <- as.data.frame(fixef(M_session2))
fixef_s3 <- as.data.frame(fixef(M_session3))

# Combine estimates into a table
fixef_table <- tibble(
  Parameter = rownames(fixef_s1),  
  Session_1_Estimate = fixef_s1$Estimate,
  Session_1_CI_Lower = fixef_s1$`Q2.5`,
  Session_1_CI_Upper = fixef_s1$`Q97.5`,
  
  Session_2_Estimate = fixef_s2$Estimate,
  Session_2_CI_Lower = fixef_s2$`Q2.5`,
  Session_2_CI_Upper = fixef_s2$`Q97.5`,
  
  Session_3_Estimate = fixef_s3$Estimate,
  Session_3_CI_Lower = fixef_s3$`Q2.5`,
  Session_3_CI_Upper = fixef_s3$`Q97.5`
)

# Print table
print(fixef_table)

# A tibble: 3 × 10
  Parameter      Session_1_Estimate Session_1_CI_Lower Session_1_CI_Upper
  <chr>                       <dbl>              <dbl>              <dbl>
1 ampl_Intercept              1.06               0.380              2.01 
2 rate_Intercept              0.368              0.134              0.642
3 asym_Intercept              0.464              0.371              0.554
# ℹ 6 more variables: Session_2_Estimate <dbl>, Session_2_CI_Lower <dbl>,
#   Session_2_CI_Upper <dbl>, Session_3_Estimate <dbl>,
#   Session_3_CI_Lower <dbl>, Session_3_CI_Upper <dbl>

# Save as CSV if needed
write.csv(fixef_table, "fixef_table_full.csv", row.names = FALSE)

#Random effects for each session

#session1
# Extract random effects for the "subject" grouping factor
ranef_ms1 <- ranef(M_session1)$subject %>%  
  as.data.frame() %>%
  rownames_to_column(var = "subject")  # Convert row names to a column

# Rename columns for clarity (adjust based on `ranef(M_1)` output)
colnames(ranef_ms1) <- c("subject", "term", "Estimate", "Q2.5", "Q97.5")

# Convert subject to a factor for clear axis labels
ranef_ms1$subject <- as.factor(ranef_ms1$subject)

# Create a single plot with all subjects
ggplot(ranef_ms1, aes(x = subject,  
                      y = Estimate,  
                      ymin = Q2.5,  
                      ymax = Q97.5)) +  
  geom_crossbar(width = 0.4, color = "blue", fill = "lightblue") +  
  geom_point(size = 3, color = "black") +  # Highlight mean estimates
  theme_minimal() +
  labs(title = "Subject-Level Random Effects",
       x = "Participant",
       y = "Random Effect Estimate") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))  # Rotate x-axis labels

#Session2
# Extract random effects for the "subject" grouping factor
ranef_ms2 <- ranef(M_session2)$subject %>%  
  as.data.frame() %>%
  rownames_to_column(var = "subject")  # Convert row names to a column

# Rename columns for clarity (adjust based on `ranef(M_1)` output)
colnames(ranef_ms2) <- c("subject", "term", "Estimate", "Q2.5", "Q97.5")

# Convert subject to a factor for clear axis labels
ranef_ms2$subject <- as.factor(ranef_ms2$subject)

# Create a single plot with all subjects
ggplot(ranef_ms2, aes(x = subject,  
                      y = Estimate,  
                      ymin = Q2.5,  
                      ymax = Q97.5)) +  
  geom_crossbar(width = 0.4, color = "blue", fill = "lightblue") +  
  geom_point(size = 3, color = "black") +  # Highlight mean estimates
  theme_minimal() +
  labs(title = "Subject-Level Random Effects",
       x = "Participant",
       y = "Random Effect Estimate") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))  # Rotate x-axis labels

# Session3
# Extract random effects for the "subject" grouping factor
ranef_ms3 <- ranef(M_session3)$subject %>%  
  as.data.frame() %>%
  rownames_to_column(var = "subject")  # Convert row names to a column

# Rename columns for clarity (adjust based on `ranef(M_1)` output)
colnames(ranef_ms3) <- c("subject", "term", "Estimate", "Q2.5", "Q97.5")

# Convert subject to a factor for clear axis labels
ranef_ms3$subject <- as.factor(ranef_ms3$subject)

# Create a single plot with all subjects
ggplot(ranef_ms3, aes(x = subject,  
                      y = Estimate,  
                      ymin = Q2.5,  
                      ymax = Q97.5)) +  
  geom_crossbar(width = 0.4, color = "blue", fill = "lightblue") +  
  geom_point(size = 3, color = "black") +  # Highlight mean estimates
  theme_minimal() +
  labs(title = "Subject-Level Random Effects",
       x = "Participant",
       y = "Random Effect Estimate") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))  # Rotate x-axis labels

# Create separate datasets for each session
df_session1 <- df_lc %>% filter(session == 1)
df_session2 <- df_lc %>% filter(session == 2)
df_session3 <- df_lc %>% filter(session == 3)

# Save them as separate CSV files
write.csv(df_session1, "df_session1.csv", row.names = FALSE)
write.csv(df_session2, "df_session2.csv", row.names = FALSE)
write.csv(df_session3, "df_session3.csv", row.names = FALSE)

#session 1

# Generate fitted values (predicted learning curves)
fitted_ms1 <- fitted(M_session1, summary = TRUE) %>%  # Extract fitted responses
  as.data.frame()

# Ensure the dataset has subject & trial numbers
fitted_ms1 <- df_session1 %>% 
  select(subject, trial) %>%  # Keep only subject & trial
  bind_cols(fitted_ms1)  # Merge with fitted values

# Plot estimated learning curves
ggplot(fitted_ms1, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curves
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # 95% Credible Interval
  labs(title = "Estimated Learning Curves",
       x = "Trial Number",
       y = "Predicted RT (s)",
       color = "Subject") +
  theme_minimal()

ggplot() +
  geom_line(data = df_session1, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +  
  geom_line(data = fitted_ms1, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_ms1, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  labs(title = "Observed vs. Estimated Learning Curves",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curves for each subject in separate panels
ggplot(fitted_ms1, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curve
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # Confidence band
  facet_wrap(~subject, scales = "free_y") +  # Separate plot per subject
  labs(title = "Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Predicted RT (s)") +
  theme_minimal() +
  theme(legend.position = "none")  

# Plot observed vs. estimated learning curves for each subject
ggplot() +
  geom_line(data = df_session1, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +
  geom_line(data = fitted_ms1, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_ms1, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  facet_wrap(~subject, scales = "free_y") +
  labs(title = "Observed vs. Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Reaction Time (s)",
       color = "Subject") +
  theme_minimal() +
  theme(legend.position = "none")

#session 2

# Generate fitted values (predicted learning curves)
fitted_ms2 <- fitted(M_session2, summary = TRUE) %>%  # Extract fitted responses
  as.data.frame()

# Ensure the dataset has subject & trial numbers
fitted_ms2 <- df_session2 %>% 
  select(subject, trial) %>%  # Keep only subject & trial
  bind_cols(fitted_ms2)  # Merge with fitted values

# Plot estimated learning curves
ggplot(fitted_ms2, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curves
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # 95% Credible Interval
  labs(title = "Estimated Learning Curves",
       x = "Trial Number",
       y = "Predicted RT (s)",
       color = "Subject") +
  theme_minimal()

ggplot() +
  geom_line(data = df_session2, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +  
  geom_line(data = fitted_ms2, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_ms2, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  labs(title = "Observed vs. Estimated Learning Curves",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curves for each subject in separate panels
ggplot(fitted_ms2, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curve
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # Confidence band
  facet_wrap(~subject, scales = "free_y") +  # Separate plot per subject
  labs(title = "Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Predicted RT (s)") +
  theme_minimal() +
  theme(legend.position = "none")  

# Plot observed vs. estimated learning curves for each subject
ggplot() +
  geom_line(data = df_session2, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +
  geom_line(data = fitted_ms2, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_ms2, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  facet_wrap(~subject, scales = "free_y") +
  labs(title = "Observed vs. Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Reaction Time (s)",
       color = "Subject") +
  theme_minimal() +
  theme(legend.position = "none")

#session 3

# Generate fitted values (predicted learning curves)
fitted_ms3 <- fitted(M_session3, summary = TRUE) %>%  # Extract fitted responses
  as.data.frame()

# Ensure the dataset has subject & trial numbers
fitted_ms3 <- df_session3 %>% 
  select(subject, trial) %>%  # Keep only subject & trial
  bind_cols(fitted_ms3)  # Merge with fitted values

# Plot estimated learning curves
ggplot(fitted_ms3, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curves
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # 95% Credible Interval
  labs(title = "Estimated Learning Curves",
       x = "Trial Number",
       y = "Predicted RT (s)",
       color = "Subject") +
  theme_minimal()

ggplot() +
  geom_line(data = df_session3, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +  
  geom_line(data = fitted_ms3, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_ms3, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  labs(title = "Observed vs. Estimated Learning Curves",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curves for each subject in separate panels
ggplot(fitted_ms3, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curve
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # Confidence band
  facet_wrap(~subject, scales = "free_y") +  # Separate plot per subject
  labs(title = "Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Predicted RT (s)") +
  theme_minimal() +
  theme(legend.position = "none")  

# Plot observed vs. estimated learning curves for each subject
ggplot() +
  geom_line(data = df_session3, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +
  geom_line(data = fitted_ms3, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_ms3, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  facet_wrap(~subject, scales = "free_y") +
  labs(title = "Observed vs. Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Reaction Time (s)",
       color = "Subject") +
  theme_minimal() +
  theme(legend.position = "none")

#Session 1 group

# Ensure fitted values contain trial information
fitted_ms1 <- fitted(M_session1, summary = TRUE) %>%
  as.data.frame() %>%
  mutate(trial = df_session1$trial)  # Explicitly add trial column from df_lc

# Aggregate model predictions at the group level
group_fitted_ms1 <- fitted_ms1 %>%
  group_by(trial) %>%
  summarise(
    mean_Estimate = mean(Estimate, na.rm = TRUE),
    lower_CI = mean(Q2.5, na.rm = TRUE),
    upper_CI = mean(Q97.5, na.rm = TRUE)
  )


# Aggregate observed data at the group level
group_obs_ms1 <- df_session1 %>%
  group_by(trial) %>%
  summarise(
    mean_RT = mean(trial.RTS, na.rm = TRUE),
    se_RT = sd(trial.RTS, na.rm = TRUE) / sqrt(n())
  )




# Plot group-level learning curve
ggplot() +
  # Observed mean RT with standard error
  geom_line(data = group_obs_ms1, aes(x = trial, y = mean_RT), color = "black", size = 1, alpha = 0.5) +
  geom_ribbon(data = group_obs_ms1, aes(x = trial, ymin = mean_RT - se_RT, ymax = mean_RT + se_RT), alpha = 0.2, fill = "grey") +

  # Estimated mean RT from model
  geom_line(data = group_fitted_ms1, aes(x = trial, y = mean_Estimate), color = "blue", size = 1) +
  geom_ribbon(data = group_fitted_ms1, aes(x = trial, ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +

  # Labels and theme
  labs(title = "Group-Level Learning Curve: Observed vs. Estimated RT",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curve (group-level)
ggplot(group_fitted_ms1, aes(x = trial, y = mean_Estimate)) +
  geom_line(color = "blue", size = 1) +  # Estimated RT line
  geom_ribbon(aes(ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +  # Confidence interval
  labs(title = "Group-Level Estimated Learning Curve",
       x = "Trial Number",
       y = "Predicted Reaction Time (s)") +
  theme_minimal()

#Session 2 group

# Ensure fitted values contain trial information
fitted_ms2 <- fitted(M_session2, summary = TRUE) %>%
  as.data.frame() %>%
  mutate(trial = df_session2$trial)  # Explicitly add trial column from df_lc

# Aggregate model predictions at the group level
group_fitted_ms2 <- fitted_ms2 %>%
  group_by(trial) %>%
  summarise(
    mean_Estimate = mean(Estimate, na.rm = TRUE),
    lower_CI = mean(Q2.5, na.rm = TRUE),
    upper_CI = mean(Q97.5, na.rm = TRUE)
  )


# Aggregate observed data at the group level
group_obs_ms2 <- df_session2 %>%
  group_by(trial) %>%
  summarise(
    mean_RT = mean(trial.RTS, na.rm = TRUE),
    se_RT = sd(trial.RTS, na.rm = TRUE) / sqrt(n())
  )




# Plot group-level learning curve
ggplot() +
  # Observed mean RT with standard error
  geom_line(data = group_obs_ms2, aes(x = trial, y = mean_RT), color = "black", size = 1, alpha = 0.5) +
  geom_ribbon(data = group_obs_ms2, aes(x = trial, ymin = mean_RT - se_RT, ymax = mean_RT + se_RT), alpha = 0.2, fill = "grey") +

  # Estimated mean RT from model
  geom_line(data = group_fitted_ms2, aes(x = trial, y = mean_Estimate), color = "blue", size = 1) +
  geom_ribbon(data = group_fitted_ms2, aes(x = trial, ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +

  # Labels and theme
  labs(title = "Group-Level Learning Curve: Observed vs. Estimated RT",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curve (group-level)
ggplot(group_fitted_ms2, aes(x = trial, y = mean_Estimate)) +
  geom_line(color = "blue", size = 1) +  # Estimated RT line
  geom_ribbon(aes(ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +  # Confidence interval
  labs(title = "Group-Level Estimated Learning Curve",
       x = "Trial Number",
       y = "Predicted Reaction Time (s)") +
  theme_minimal()

#Session 3 group

# Ensure fitted values contain trial information
fitted_ms3 <- fitted(M_session3, summary = TRUE) %>%
  as.data.frame() %>%
  mutate(trial = df_session3$trial)  # Explicitly add trial column from df_lc

# Aggregate model predictions at the group level
group_fitted_ms3 <- fitted_ms3 %>%
  group_by(trial) %>%
  summarise(
    mean_Estimate = mean(Estimate, na.rm = TRUE),
    lower_CI = mean(Q2.5, na.rm = TRUE),
    upper_CI = mean(Q97.5, na.rm = TRUE)
  )


# Aggregate observed data at the group level
group_obs_ms3 <- df_session3 %>%
  group_by(trial) %>%
  summarise(
    mean_RT = mean(trial.RTS, na.rm = TRUE),
    se_RT = sd(trial.RTS, na.rm = TRUE) / sqrt(n())
  )




# Plot group-level learning curve
ggplot() +
  # Observed mean RT with standard error
  geom_line(data = group_obs_ms3, aes(x = trial, y = mean_RT), color = "black", size = 1, alpha = 0.5) +
  geom_ribbon(data = group_obs_ms3, aes(x = trial, ymin = mean_RT - se_RT, ymax = mean_RT + se_RT), alpha = 0.2, fill = "grey") +

  # Estimated mean RT from model
  geom_line(data = group_fitted_ms3, aes(x = trial, y = mean_Estimate), color = "blue", size = 1) +
  geom_ribbon(data = group_fitted_ms3, aes(x = trial, ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +

  # Labels and theme
  labs(title = "Group-Level Learning Curve: Observed vs. Estimated RT",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curve (group-level)
ggplot(group_fitted_ms3, aes(x = trial, y = mean_Estimate)) +
  geom_line(color = "blue", size = 1) +  # Estimated RT line
  geom_ribbon(aes(ymin = lower_CI, ymax = upper_CI), alpha = 0.2, fill = "blue") +  # Confidence interval
  labs(title = "Group-Level Estimated Learning Curve",
       x = "Trial Number",
       y = "Predicted Reaction Time (s)") +
  theme_minimal()

#Inverse decay function tryout

# Define the inverse decay function
F_inv_decay <- formula(trial.RTS ~ baseline + (equilibrium - baseline) * exp(-speed * trial))

# Define hierarchical effects
F_inv_decay_ef_1 <- list(
  formula(baseline ~ 1|subject),
  formula(equilibrium ~ 1|subject),
  formula(speed ~ 1|subject)
)

# Define priors
F_inv_decay_prior <- c(
  set_prior("normal(3, 20)", nlpar = "baseline", lb = 0),
  set_prior("normal(5, 100)", nlpar = "equilibrium", lb = 0),
  set_prior("normal(0.5, 3)", nlpar = "speed", lb = 0)
)

# Fit Bayesian model
M_inv_decay <- 
  df_lc1 %>% 
  brm(
    bf(F_inv_decay, flist = F_inv_decay_ef_1, nl = TRUE), 
    prior = F_inv_decay_prior,
    family = "exgaussian",
    data = .,
    iter = 5000,
    warmup = 4000,
    save_pars = save_pars(all = TRUE)
  )

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.6)’
using SDK: ‘’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: 
Chain 1: Gradient evaluation took 0.00056 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 5.6 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 1: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 2447.55 seconds (Warm-up)
Chain 1:                270.308 seconds (Sampling)
Chain 1:                2717.86 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: Rejecting initial value:
Chain 2:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 2:   Stan can't start sampling from this initial value.
Chain 2: 
Chain 2: Gradient evaluation took 0.000314 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.14 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 2: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 1216.03 seconds (Warm-up)
Chain 2:                762.794 seconds (Sampling)
Chain 2:                1978.82 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000312 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.12 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 3: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 989.054 seconds (Warm-up)
Chain 3:                305.408 seconds (Sampling)
Chain 3:                1294.46 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: Rejecting initial value:
Chain 4:   Log probability evaluates to log(0), i.e. negative infinity.
Chain 4:   Stan can't start sampling from this initial value.
Chain 4: 
Chain 4: Gradient evaluation took 0.000312 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 3.12 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 4: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 1152.28 seconds (Warm-up)
Chain 4:                326.029 seconds (Sampling)
Chain 4:                1478.31 seconds (Total)
Chain 4: 

Warning: There were 132 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: There were 3633 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: The largest R-hat is 1.57, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

# Extract posterior samples
P_inv_decay <- posterior(M_inv_decay) 
PP_inv_decay <- post_pred(M_inv_decay, thin = 10)

# Save results
save(M_inv_decay, P_inv_decay, PP_inv_decay, df_lc1, file = "2025_analyses_long_inv_decay.Rda")

# Summarize fixed effects
fixef(M_inv_decay)

                       Estimate Est.Error       Q2.5     Q97.5
baseline_Intercept    0.6116028 0.1120223 0.47484699 0.9096236
equilibrium_Intercept 1.3799560 0.5630718 0.53115535 2.3879450
speed_Intercept       0.2644908 0.1515227 0.01738646 0.5146940

# Summarize group-level effects
grpef(M_inv_decay)

Coefficient estimates with 95% credibility limits

nonlin	center	lower	upper
baseline	0.2197496	0.1513009	0.5821730
equilibrium	1.3583441	0.2577583	2.4847778
speed	0.2843252	0.0393428	0.5013047

ranef(M_inv_decay)

$subject
, , baseline_Intercept

      Estimate Est.Error        Q2.5        Q97.5
2  -0.10252994 0.1126292 -0.39732458  0.039956429
3  -0.06085510 0.1197214 -0.25388885  0.255605485
4   0.15178613 0.1277062 -0.13768365  0.388856708
5   0.09082584 0.2089260 -0.09518495  0.668845107
7  -0.16370035 0.1127753 -0.46259307 -0.025573624
8  -0.19162561 0.1149513 -0.49448760 -0.048997250
10 -0.27404444 0.1695965 -0.63606359  0.186764798
11 -0.18928816 0.1713227 -0.50717402  0.295643455
12 -0.27062755 0.2026075 -0.78177657 -0.013540384
13 -0.17270274 0.1099302 -0.46433923 -0.036023128
14 -0.13544141 0.1084924 -0.41996269  0.005895949
15  0.01430032 0.1157016 -0.28812604  0.155068188
16 -0.01033380 0.1281181 -0.17159556  0.352119895
17  0.02434837 0.1770176 -0.14485926  0.564004794
18  0.25579552 0.2074167  0.07227594  0.836625751
19  0.19784095 0.2818285 -0.03573493  0.996156879
20  0.47848240 0.1178386  0.17067640  0.625081351
22  0.27932486 0.1332058  0.12562537  0.658845250
23 -0.03290774 0.1676673 -0.31823203  0.357711659

, , equilibrium_Intercept

      Estimate Est.Error       Q2.5       Q97.5
2  -0.06858399 0.5519492 -1.0665251  0.77256051
3  -0.39886908 0.3784910 -1.2051849  0.29905457
4  -0.89014961 0.5715323 -1.8935435 -0.01491107
5   0.91508687 0.8890968 -0.3264513  2.59876564
7  -0.92003313 0.5963420 -2.0148130 -0.02486773
8  -0.17587727 0.5340799 -1.1880046  0.65872788
10 -0.75775052 0.4669414 -1.6971290  0.01141939
11 -0.15337098 0.5355943 -1.1634866  1.02388731
12  0.36720972 1.4968690 -2.8109498  3.42895585
13 -0.57652827 0.5601028 -1.5889147  0.28155035
14 -1.05749744 0.5770725 -2.0826634 -0.16814493
15 -0.30271306 0.5389564 -1.2917037  0.53360102
16  0.12635918 0.4811461 -0.6617160  1.20811641
17 -0.27411902 0.6277614 -1.3134840  1.45742085
18  0.42806815 0.5380078 -0.4298879  1.64846730
19  1.30126186 1.3931243 -0.3289615  4.54098969
20  3.35615104 2.1265922 -0.2085473  6.22415371
22  0.26251022 0.4742963 -0.5828039  1.33995956
23 -1.11447059 0.5699350 -2.1281051 -0.24426488

, , speed_Intercept

      Estimate  Est.Error        Q2.5        Q97.5
2  -0.22322789 0.15113557 -0.47300933  0.023836712
3  -0.10772990 0.09511469 -0.31906200  0.032741673
4  -0.25012747 0.14961346 -0.49810023 -0.004925643
5   0.27289306 0.21665633 -0.04738386  0.656262309
7   0.11243878 0.24599962 -0.28793580  0.698163009
8  -0.09345069 0.14032538 -0.33888138  0.138540421
10  0.16769231 0.25469458 -0.17745835  0.783188637
11  0.26376400 0.25388586 -0.06118711  0.776912981
12 -0.04109912 0.24202261 -0.41479359  0.550123394
13 -0.17477127 0.14864499 -0.42544704  0.069873643
14 -0.17947403 0.20259893 -0.45883063  0.300236146
15 -0.10258419 0.14555093 -0.35962040  0.136941404
16  0.18639561 0.18738576 -0.05448605  0.582773804
17 -0.15081806 0.16625481 -0.41603573  0.304279575
18  0.03308517 0.12111536 -0.17689516  0.303574776
19  0.37389400 0.28647025 -0.04708342  0.863492403
20  0.18793638 0.11927546  0.01885704  0.418388151
22  0.03418489 0.13329001 -0.17654642  0.357450386
23 -0.25630275 0.15017511 -0.50530520 -0.010958965

# Extract random effects for the "subject" grouping factor
ranef_M_inv_decay <- ranef(M_inv_decay)$subject %>%  
  as.data.frame() %>%
  rownames_to_column(var = "subject")  # Convert row names to a column

# Rename columns for clarity (adjust based on `ranef(M_1)` output)
colnames(ranef_M_inv_decay) <- c("subject", "term", "Estimate", "Q2.5", "Q97.5")

# Convert subject to a factor for clear axis labels
ranef_M_inv_decay$subject <- as.factor(ranef_M_inv_decay$subject)

# Create a single plot with all subjects
ggplot(ranef_M_inv_decay, aes(x = subject,  
                      y = Estimate,  
                      ymin = Q2.5,  
                      ymax = Q97.5)) +  
  geom_crossbar(width = 0.4, color = "blue", fill = "lightblue") +  
  geom_point(size = 3, color = "black") +  # Highlight mean estimates
  theme_minimal() +
  labs(title = "Subject-Level Random Effects",
       x = "Participant",
       y = "Random Effect Estimate") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))  # Rotate x-axis labels

# Generate fitted values (predicted learning curves)
fitted_M_inv_decay <- fitted(M_inv_decay, summary = TRUE) %>%  # Extract fitted responses
  as.data.frame()

# Ensure the dataset has subject & trial numbers
fitted_M_inv_decay <- df_lc %>% 
  select(subject, trial) %>%  # Keep only subject & trial
  bind_cols(fitted_M_inv_decay)  # Merge with fitted values

# Plot estimated learning curves
ggplot(fitted_M_inv_decay, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curves
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # 95% Credible Interval
  labs(title = "Estimated Learning Curves",
       x = "Trial Number",
       y = "Predicted RT (s)",
       color = "Subject") +
  theme_minimal()

ggplot() +
  geom_line(data = df_lc, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +  
  geom_line(data = fitted_M_inv_decay, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_M_inv_decay, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  labs(title = "Observed vs. Estimated Learning Curves",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curves for each subject in separate panels
ggplot(fitted_M_inv_decay, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curve
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # Confidence band
  facet_wrap(~subject, scales = "free_y") +  # Separate plot per subject
  labs(title = "Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Predicted RT (s)") +
  theme_minimal() +
  theme(legend.position = "none")  

# Plot observed vs. estimated learning curves for each subject
ggplot() +
  geom_line(data = df_lc, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +
  geom_line(data = fitted_M_inv_decay, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_M_inv_decay, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  facet_wrap(~subject, scales = "free_y") +
  labs(title = "Observed vs. Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Reaction Time (s)",
       color = "Subject") +
  theme_minimal() +
  theme(legend.position = "none")

#logistic growth function tryout

# Define the logistic growth function
F_logistic <- formula(trial.RTS ~ equilibrium / (1 + exp(-speed * (trial - halftime))))

# Define hierarchical effects (subject-level variations)
F_logistic_ef_1 <- list(
  formula(equilibrium ~ 1 | subject),
  formula(speed ~ 1 | subject),
  formula(halftime ~ 1 | subject)
)

# Define priors for Bayesian estimation
F_logistic_prior <- c(
  set_prior("normal(3, 20)", nlpar = "equilibrium", lb = 0),
  set_prior("normal(0.5, 3)", nlpar = "speed", lb = 0),
  set_prior("normal(50, 30)", nlpar = "halftime", lb = 0)
)

# Fit Bayesian model using logistic growth function
M_logistic <- 
  df_lc1 %>% 
  brm(
    bf(F_logistic, flist = F_logistic_ef_1, nl = TRUE), 
    prior = F_logistic_prior,
    family = "exgaussian",
    data = .,
    iter = 5000,
    warmup = 4000,
    save_pars = save_pars(all = TRUE)
  )

Compiling Stan program...

Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c
using C compiler: ‘Apple clang version 16.0.0 (clang-1600.0.26.6)’
using SDK: ‘’
clang -arch arm64 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG   -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/Rcpp/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/unsupported"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/BH/include" -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/src/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppParallel/include/"  -I"/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/rstan/include" -DEIGEN_NO_DEBUG  -DBOOST_DISABLE_ASSERTS  -DBOOST_PENDING_INTEGER_LOG2_HPP  -DSTAN_THREADS  -DUSE_STANC3 -DSTRICT_R_HEADERS  -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION  -D_HAS_AUTO_PTR_ETC=0  -include '/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp'  -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1   -I/opt/R/arm64/include    -fPIC  -falign-functions=64 -Wall -g -O2  -c foo.c -o foo.o
In file included from <built-in>:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1:
In file included from /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/Core:19:
/Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: 'cmath' file not found
  679 | #include <cmath>
      |          ^~~~~~~
1 error generated.
make: *** [foo.o] Error 1

Start sampling

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 1).
Chain 1: Rejecting initial value:
Chain 1:   Gradient evaluated at the initial value is not finite.
Chain 1:   Stan can't start sampling from this initial value.
Chain 1: 
Chain 1: Gradient evaluation took 0.000314 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.14 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1: 
Chain 1: 
Chain 1: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 1: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 1: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 1: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 1: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 1: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 1: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 1: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 1: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 1: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 1: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 1: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 1: 
Chain 1:  Elapsed Time: 1350.38 seconds (Warm-up)
Chain 1:                323.733 seconds (Sampling)
Chain 1:                1674.12 seconds (Total)
Chain 1: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 2).
Chain 2: 
Chain 2: Gradient evaluation took 0.000312 seconds
Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 3.12 seconds.
Chain 2: Adjust your expectations accordingly!
Chain 2: 
Chain 2: 
Chain 2: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 2: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 2: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 2: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 2: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 2: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 2: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 2: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 2: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 2: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 2: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 2: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 2: 
Chain 2:  Elapsed Time: 125.918 seconds (Warm-up)
Chain 2:                13.135 seconds (Sampling)
Chain 2:                139.053 seconds (Total)
Chain 2: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 3).
Chain 3: 
Chain 3: Gradient evaluation took 0.000314 seconds
Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 3.14 seconds.
Chain 3: Adjust your expectations accordingly!
Chain 3: 
Chain 3: 
Chain 3: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 3: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 3: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 3: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 3: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 3: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 3: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 3: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 3: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 3: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 3: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 3: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 3: 
Chain 3:  Elapsed Time: 1213.19 seconds (Warm-up)
Chain 3:                317.247 seconds (Sampling)
Chain 3:                1530.44 seconds (Total)
Chain 3: 

SAMPLING FOR MODEL 'anon_model' NOW (CHAIN 4).
Chain 4: 
Chain 4: Gradient evaluation took 0.000306 seconds
Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 3.06 seconds.
Chain 4: Adjust your expectations accordingly!
Chain 4: 
Chain 4: 
Chain 4: Iteration:    1 / 5000 [  0%]  (Warmup)
Chain 4: Iteration:  500 / 5000 [ 10%]  (Warmup)
Chain 4: Iteration: 1000 / 5000 [ 20%]  (Warmup)
Chain 4: Iteration: 1500 / 5000 [ 30%]  (Warmup)
Chain 4: Iteration: 2000 / 5000 [ 40%]  (Warmup)
Chain 4: Iteration: 2500 / 5000 [ 50%]  (Warmup)
Chain 4: Iteration: 3000 / 5000 [ 60%]  (Warmup)
Chain 4: Iteration: 3500 / 5000 [ 70%]  (Warmup)
Chain 4: Iteration: 4000 / 5000 [ 80%]  (Warmup)
Chain 4: Iteration: 4001 / 5000 [ 80%]  (Sampling)
Chain 4: Iteration: 4500 / 5000 [ 90%]  (Sampling)
Chain 4: Iteration: 5000 / 5000 [100%]  (Sampling)
Chain 4: 
Chain 4:  Elapsed Time: 109.217 seconds (Warm-up)
Chain 4:                26.726 seconds (Sampling)
Chain 4:                135.943 seconds (Total)
Chain 4: 

Warning: There were 619 divergent transitions after warmup. See
https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
to find out why this is a problem and how to eliminate them.

Warning: There were 1997 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded

Warning: There were 1 chains where the estimated Bayesian Fraction of Missing Information was low. See
https://mc-stan.org/misc/warnings.html#bfmi-low

Warning: Examine the pairs() plot to diagnose sampling problems

Warning: The largest R-hat is 2.48, indicating chains have not mixed.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#r-hat

Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#bulk-ess

Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
Running the chains for more iterations may help. See
https://mc-stan.org/misc/warnings.html#tail-ess

# Extract posterior samples
P_logistic <- posterior(M_logistic) 
PP_logistic <- post_pred(M_logistic, thin = 10)

# Save results
save(M_logistic, P_logistic, PP_logistic, df_lc1, file = "2025_analyses_long_logistic.Rda")

# Summarize fixed effects
fixef(M_logistic)

                       Estimate  Est.Error       Q2.5     Q97.5
equilibrium_Intercept 0.5990902 0.06970355 0.49847265 0.7795943
speed_Intercept       4.6284772 2.58915166 0.02852926 8.7349549
halftime_Intercept    1.0578903 2.63367949 0.01472050 9.2354361

# Summarize group-level effects
grpef(M_logistic)

Coefficient estimates with 95% credibility limits

nonlin	center	lower	upper
equilibrium	0.2195314	0.1519384	0.353103
speed	1.7712805	0.0278271	4.688802
halftime	0.3896337	0.0241850	74.586228

ranef(M_logistic)

$subject
, , equilibrium_Intercept

      Estimate  Est.Error        Q2.5       Q97.5
2  -0.02143697 0.07320593 -0.21080624  0.08373278
3  -0.06118754 0.05575555 -0.14430956  0.07336024
4   0.08089777 0.05954229 -0.05661963  0.18261790
5   0.04449573 0.11187305 -0.04863988  0.41070696
7  -0.13306849 0.06947234 -0.30637274 -0.01812806
8  -0.09260144 0.08892375 -0.17893826  0.16828401
10 -0.26480198 0.06520325 -0.34211716 -0.08937996
11 -0.18416302 0.07154941 -0.27684129  0.00793487
12 -0.20138198 0.07518132 -0.38251513 -0.08973074
13 -0.12871244 0.07079976 -0.31417018 -0.02152292
14 -0.11759721 0.06281862 -0.28230920 -0.02273443
15  0.04930514 0.07407181 -0.14171209  0.15413620
16  0.06941437 0.15608412 -0.10950544  0.37978437
17 -0.04585003 0.07146326 -0.23150490  0.05652511
18  0.21116241 0.10546286  0.10627379  0.51740685
19  0.14600030 0.16169862  0.01466776  0.53631332
20  0.49585084 0.07919193  0.30064576  0.61452467
22  0.26888397 0.09418105  0.18208753  0.57467567
23 -0.12116911 0.11020520 -0.26461903  0.08645757

, , speed_Intercept

      Estimate Est.Error      Q2.5    Q97.5
2  -0.31557725  1.974430 -4.473506 2.600323
3  -0.54582466  1.800025 -3.396401 2.974891
4  -2.51333417  2.980010 -6.509792 1.734907
5   0.19797270  1.384046 -2.774114 3.534131
7   0.54725211  2.192486 -3.815936 4.067291
8  -0.41623020  1.982819 -4.159023 4.396444
10  0.52280961  2.023222 -2.554069 6.181857
11  0.29336723  1.751274 -3.024692 3.194810
12  1.47158976  2.219032 -1.612658 5.189496
13  0.26320476  1.514965 -2.589887 2.506728
14 -1.71426380  2.452831 -6.469921 1.987494
15  0.19809219  1.764490 -3.397079 4.565384
16 -1.25527983  2.828419 -6.560174 2.768247
17 -0.46314723  1.646631 -4.981210 2.255215
18  0.02572509  1.710541 -2.586214 4.048392
19 -0.07943453  1.359379 -2.884963 3.578209
20 -0.17191876  1.515689 -3.637640 2.609714
22  0.28452085  1.894519 -3.890838 4.603135
23 -1.31715476  2.865878 -6.565804 1.781547

, , halftime_Intercept

       Estimate Est.Error         Q2.5      Q97.5
2  -14.13058677 35.124564 -118.5024539  0.3047669
3   -6.65289987 18.233094  -62.1986754  0.6078313
4  -10.62237228 22.542290  -73.6518598  0.5359371
5    2.41666616 14.548758  -18.5859955 53.1483308
7  -13.04379502 32.792716 -113.5953858  0.3147363
8   -1.04391171 19.527187  -47.5578349 47.3505591
10   3.85983554 23.926260  -32.1475387 77.4551578
11  -3.18752318 19.686998  -63.7463242 34.3215870
12  -1.78585796 20.365420  -62.6529571 38.2629947
13 -11.72394925 27.599628  -95.9856852  0.4811219
14  -6.24637515 15.085942  -52.8994791  0.3300736
15 -14.26041599 34.370688 -120.3663730  0.2676684
16  -6.32192568 17.426417  -61.5366400  1.6392452
17 -14.33011551 36.692886 -124.0181524  0.3148622
18  -6.50873759 20.159435  -71.2474463  3.0583192
19   5.79439281 13.906268   -0.5542792 50.4151845
20  -2.10344320  6.490356  -21.1250729  0.6517338
22  -5.58138779 12.839403  -43.4501258  0.3633626
23   0.03964254 15.823954  -37.1484366 39.4762393

# Extract random effects for the "subject" grouping factor
ranef_M_logistic <- ranef(M_logistic)$subject %>%  
  as.data.frame() %>%
  rownames_to_column(var = "subject")  # Convert row names to a column

# Rename columns for clarity (adjust based on `ranef(M_1)` output)
colnames(ranef_M_logistic) <- c("subject", "term", "Estimate", "Q2.5", "Q97.5")

# Convert subject to a factor for clear axis labels
ranef_M_logistic$subject <- as.factor(ranef_M_logistic$subject)

# Create a single plot with all subjects
ggplot(ranef_M_logistic, aes(x = subject,  
                      y = Estimate,  
                      ymin = Q2.5,  
                      ymax = Q97.5)) +  
  geom_crossbar(width = 0.4, color = "blue", fill = "lightblue") +  
  geom_point(size = 3, color = "black") +  # Highlight mean estimates
  theme_minimal() +
  labs(title = "Subject-Level Random Effects",
       x = "Participant",
       y = "Random Effect Estimate") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))  # Rotate x-axis labels

# Generate fitted values (predicted learning curves)
fitted_M_logistic <- fitted(M_logistic, summary = TRUE) %>%  # Extract fitted responses
  as.data.frame()

# Ensure the dataset has subject & trial numbers
fitted_M_logistic <- df_lc %>% 
  select(subject, trial) %>%  # Keep only subject & trial
  bind_cols(fitted_M_logistic)  # Merge with fitted values

# Plot estimated learning curves
ggplot(fitted_M_logistic, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curves
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # 95% Credible Interval
  labs(title = "Estimated Learning Curves",
       x = "Trial Number",
       y = "Predicted RT (s)",
       color = "Subject") +
  theme_minimal()

ggplot() +
  geom_line(data = df_lc, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +  
  geom_line(data = fitted_M_logistic, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_M_logistic, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  labs(title = "Observed vs. Estimated Learning Curves",
       x = "Trial Number",
       y = "Reaction Time (s)") +
  theme_minimal()

# Plot estimated learning curves for each subject in separate panels
ggplot(fitted_M_logistic, aes(x = trial, y = Estimate, color = as.factor(subject))) +
  geom_line(size = 1) +  # Estimated learning curve
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5), alpha = 0.2) +  # Confidence band
  facet_wrap(~subject, scales = "free_y") +  # Separate plot per subject
  labs(title = "Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Predicted RT (s)") +
  theme_minimal() +
  theme(legend.position = "none")  

# Plot observed vs. estimated learning curves for each subject
ggplot() +
  geom_line(data = df_lc, aes(x = trial, y = trial.RTS, color = as.factor(subject)), alpha = 0.3) +
  geom_line(data = fitted_M_logistic, aes(x = trial, y = Estimate, color = as.factor(subject)), size = 1) +
  geom_ribbon(data = fitted_M_logistic, aes(x = trial, ymin = Q2.5, ymax = Q97.5, fill = as.factor(subject)), alpha = 0.2) +
  facet_wrap(~subject, scales = "free_y") +
  labs(title = "Observed vs. Estimated Learning Curves per Subject",
       x = "Trial Number",
       y = "Reaction Time (s)",
       color = "Subject") +
  theme_minimal() +
  theme(legend.position = "none")