##Set working directory

setwd("C:/Users/ChanRW/OneDrive - Universiteit Twente/2020_MSCA_IF/2_Bachelor_thesis/2022_Joel_Broncho_X/Data")

##Packages that are required for the analysis

library(tidyverse)
library(broom.mixed)
library(devtools)
library(brms)
library(mascutils)
library(bayr)
library(readxl)
library(ggplot2)
library(effects)

##Import dataset

#dfAcc1<-read.table("Acc_Task1_220124_2.csv", sep = ",", header = T, stringsAsFactors = F)
dfAcc2<-read.table("Acc_Task1_220124_4.csv", sep = ",", header = T, stringsAsFactors = F)
#df3<-read.table("Task3Behav.csv", sep = ",", header = T, stringsAsFactors = F)

##No use Free learning curves

#dfAcc1 %>% 
  #filter(RT < 10) %>% 
  #ggplot(aes(x = Time,
             #y = DomHand_X)) +
  #geom_point() +
  #geom_smooth(se = F) +
  #facet_wrap(~Participant, scales = "free_y")

#dfAcc1 %>% 
  #filter(RT < 10) %>% 
  #ggplot(aes(x = Time,
             #y = DomHand_X,
             #group = Participant)) +
  #facet_grid(~Trial) +
 #geom_smooth(se = F)

#dfAcc1 %>% 
  #filter(RT < 10) %>% 
  #ggplot(aes(x = Time,
             #y = DomHand_X,
             #group = Participant)) +
  #facet_grid(~Trial) +
  #geom_smooth(se = F)


#df %>% 
  #filter(RT < 10) %>% 
  #ggplot(aes(x = repetition,
             #y = RT,
             #group = subject)) +
  #geom_smooth(se = F)

#Real Free learning curves

dfAcc2$Participant <- factor(dfAcc2$Participant)
mycolors=c("#6b5f3c","#ccc627","#54ab8e","#587ed1","#b04366","#de2d26","#d95f0e", "#a1d99b", "#dd1c77", "#a6bddb")


## FREE LEARNING CURVES ##

# Free learning curve 1
# per subject per sequence
dfAcc2 %>% 
  ggplot(aes(x = Time,
             y = DomHand_X,
             color=Participant)) +
  geom_point() +
  geom_smooth(se = F) +
  facet_wrap(~Participant, scales = "free_y")+
  ylab("DomHand Acc")+
  xlab("Time")+
  #ylim(0,1.288)+
  theme_minimal()+
  scale_color_manual(values=mycolors)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

# Free learning curve 2
# Individual learning curves over blocks/reps

dfAcc2 %>% 
  ggplot(aes(x = Time,
             y = DomHand_X,
             group = Participant,
             color=Participant)) +
  geom_smooth(se = F)+
  #scale_x_continuous(limits = c(0,192), expand = c(0, 0))+
  theme_classic()+
  ylab("DomHand_X")+
  xlab("Time")+
  scale_color_manual(values=mycolors)+
  facet_wrap(~Trial)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

#Over Trials / TimeRepX
dfAcc2 %>% 
  ggplot(aes(x = TimeRepX,
             y = DomHand_X,
             group = Participant,
             color=Participant)) +
  geom_smooth(se = F)+
  #scale_x_continuous(limits = c(0,192), expand = c(0, 0))+
  theme_classic()+
  ylab("DomHand_X")+
  xlab("Time")+
  scale_color_manual(values=mycolors)+
  facet_wrap(~Participant)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

  #geom_vline(xintercept = c(24,48,72,96,120,144,168),colour="grey”, show.legend=TRUE)+
  #geom_text(aes(x=12, label="Trial_1", y=0.1), colour="grey") +
  #geom_text(aes(x=36, label="Trial 2", y=0.1), colour="grey") +
  #geom_text(aes(x=60, label="Trial 3", y=0.1), colour="grey") +
  #geom_text(aes(x=84, label="Trial 4", y=0.1), colour="grey") +
  #geom_text(aes(x=108, label="Trial 5", y=0.1), colour="grey")
  
  
  #geom_vline(xintercept = c(24,48,72,96,120,144,168),colour="grey”, show.legend=TRUE)+
  #geom_text(aes(x=12, label="Block 1", y=0.1), colour="grey") +
  #geom_text(aes(x=36, label="Block 2", y=0.1), colour="grey") +
  #geom_text(aes(x=60, label="Block 3", y=0.1), colour="grey") +
  #geom_text(aes(x=84, label="Block 4", y=0.1), colour="grey") +
  #geom_text(aes(x=108, label="Block 5", y=0.1), colour="grey") +
  #geom_text(aes(x=132, label="Block 6", y=0.1), colour="grey") +
  #geom_text(aes(x=156, label="Block 7", y=0.1), colour="grey") +
  #geom_text(aes(x=180, label="Block 8", y=0.1), colour="grey") 
## NON-LINEAR MULTILEVEL REGRESSION ##

# specify formula, variables and weakly informative priors
F_ary <- formula(DomHand_X ~ asym + ampl * exp(-rate * TimeRepX))

F_ary_ef_1 <- list(formula(ampl ~ 1|Participant),
                   formula(rate ~ 1|Participant),
                   formula(asym ~ 1|Participant))

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))

# create model including MCMC sampling
M_1 <- 
  dfAcc2 %>% 
  brm(bf(F_ary,
         flist = F_ary_ef_1,
         nl = T), 
      prior = F_ary_prior,
      family = "exgaussian",
      data = .,
      iter = 10, 
      warmup = 8,
      save_pars=save_pars("Participant"))
## Compiling Stan program...
## Start sampling
## 
## SAMPLING FOR MODEL 'ce11eef883174df5125e8f267d0f1a59' NOW (CHAIN 1).
## Chain 1: Rejecting initial value:
## Chain 1:   Error evaluating the log probability at the initial value.
## Chain 1: Exception: exp_mod_normal_lpdf: Location parameter[2310] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 1: Rejecting initial value:
## Chain 1:   Error evaluating the log probability at the initial value.
## Chain 1: Exception: exp_mod_normal_lpdf: Location parameter[4714] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 1: Rejecting initial value:
## Chain 1:   Error evaluating the log probability at the initial value.
## Chain 1: Exception: exp_mod_normal_lpdf: Location parameter[10666] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 1: 
## Chain 1: Gradient evaluation took 0.077 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 770 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: WARNING: No variance estimation is
## Chain 1:          performed for num_warmup < 20
## Chain 1: 
## Chain 1: Iteration: 1 / 10 [ 10%]  (Warmup)
## Chain 1: Iteration: 2 / 10 [ 20%]  (Warmup)
## Chain 1: Iteration: 3 / 10 [ 30%]  (Warmup)
## Chain 1: Iteration: 4 / 10 [ 40%]  (Warmup)
## Chain 1: Iteration: 5 / 10 [ 50%]  (Warmup)
## Chain 1: Iteration: 6 / 10 [ 60%]  (Warmup)
## Chain 1: Iteration: 7 / 10 [ 70%]  (Warmup)
## Chain 1: Iteration: 8 / 10 [ 80%]  (Warmup)
## Chain 1: Iteration: 9 / 10 [ 90%]  (Sampling)
## Chain 1: Iteration: 10 / 10 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 2.092 seconds (Warm-up)
## Chain 1:                0.313 seconds (Sampling)
## Chain 1:                2.405 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'ce11eef883174df5125e8f267d0f1a59' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 0.052 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 520 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: WARNING: No variance estimation is
## Chain 2:          performed for num_warmup < 20
## Chain 2: 
## Chain 2: Iteration: 1 / 10 [ 10%]  (Warmup)
## Chain 2: Iteration: 2 / 10 [ 20%]  (Warmup)
## Chain 2: Iteration: 3 / 10 [ 30%]  (Warmup)
## Chain 2: Iteration: 4 / 10 [ 40%]  (Warmup)
## Chain 2: Iteration: 5 / 10 [ 50%]  (Warmup)
## Chain 2: Iteration: 6 / 10 [ 60%]  (Warmup)
## Chain 2: Iteration: 7 / 10 [ 70%]  (Warmup)
## Chain 2: Iteration: 8 / 10 [ 80%]  (Warmup)
## Chain 2: Iteration: 9 / 10 [ 90%]  (Sampling)
## Chain 2: Iteration: 10 / 10 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 1.315 seconds (Warm-up)
## Chain 2:                0.378 seconds (Sampling)
## Chain 2:                1.693 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL 'ce11eef883174df5125e8f267d0f1a59' NOW (CHAIN 3).
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: exp_mod_normal_lpdf: Location parameter[40095] is -inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: exp_mod_normal_lpdf: Location parameter[4112] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 3: Rejecting initial value:
## Chain 3:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: exp_mod_normal_lpdf: Location parameter[1278] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 3: 
## Chain 3: Gradient evaluation took 0.063 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 630 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: WARNING: No variance estimation is
## Chain 3:          performed for num_warmup < 20
## Chain 3: 
## Chain 3: Iteration: 1 / 10 [ 10%]  (Warmup)
## Chain 3: Iteration: 2 / 10 [ 20%]  (Warmup)
## Chain 3: Iteration: 3 / 10 [ 30%]  (Warmup)
## Chain 3: Iteration: 4 / 10 [ 40%]  (Warmup)
## Chain 3: Iteration: 5 / 10 [ 50%]  (Warmup)
## Chain 3: Iteration: 6 / 10 [ 60%]  (Warmup)
## Chain 3: Iteration: 7 / 10 [ 70%]  (Warmup)
## Chain 3: Iteration: 8 / 10 [ 80%]  (Warmup)
## Chain 3: Iteration: 9 / 10 [ 90%]  (Sampling)
## Chain 3: Iteration: 10 / 10 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.902 seconds (Warm-up)
## Chain 3:                0.199 seconds (Sampling)
## Chain 3:                1.101 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL 'ce11eef883174df5125e8f267d0f1a59' NOW (CHAIN 4).
## 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[28300] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## 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[1] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## 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[10183] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## 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[13694] is inf, but must be finite!  (in 'modela0847d7558a_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 4: 
## Chain 4: Gradient evaluation took 0.045 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 450 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: WARNING: No variance estimation is
## Chain 4:          performed for num_warmup < 20
## Chain 4: 
## Chain 4: Iteration: 1 / 10 [ 10%]  (Warmup)
## Chain 4: Iteration: 2 / 10 [ 20%]  (Warmup)
## Chain 4: Iteration: 3 / 10 [ 30%]  (Warmup)
## Chain 4: Iteration: 4 / 10 [ 40%]  (Warmup)
## Chain 4: Iteration: 5 / 10 [ 50%]  (Warmup)
## Chain 4: Iteration: 6 / 10 [ 60%]  (Warmup)
## Chain 4: Iteration: 7 / 10 [ 70%]  (Warmup)
## Chain 4: Iteration: 8 / 10 [ 80%]  (Warmup)
## Chain 4: Iteration: 9 / 10 [ 90%]  (Sampling)
## Chain 4: Iteration: 10 / 10 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 2.724 seconds (Warm-up)
## Chain 4:                0.354 seconds (Sampling)
## Chain 4:                3.078 seconds (Total)
## Chain 4:
## Warning: There were 5 divergent transitions after warmup. See
## http://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 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
## http://mc-stan.org/misc/warnings.html#bfmi-low
## 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
## http://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
## http://mc-stan.org/misc/warnings.html#tail-ess
P_1 <- posterior(M_1) 
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
PP_1 <- post_pred(M_1, thin = 10)
## Warning: 'nsamples.brmsfit' is deprecated. Please use 'ndraws' instead.
## Warning: Argument 'nsamples' is deprecated. Please use argument 'ndraws'
## instead.
# save model estimates
save(M_1, P_1, PP_1, df, file = "model_estimates.Rda")

# get parameters for fixed effects, random factor variation and random effects
bayr::fixef(P_1)
Coefficient estimates with 95% credibility limits
nonlin center lower upper
ampl 0.6357903 0.1668794 1.0310393
rate 1.7375299 0.4383288 3.8141671
asym 0.1420251 0.0505451 0.2827376 
bayr::grpef(P_1)
Coefficient estimates with 95% credibility limits
nonlin center lower upper
ampl 0.4981550 0.2486024 5.3796286
rate 0.2305420 0.1577480 0.6702555
asym 0.2505574 0.0274151 0.7625544 
P_1 %>% re_scores() %>% bayr::ranef()
Coefficient estimates with 95% credibility limits
nonlin re_entity center lower upper
ampl 1 0.9688968 0.1967065 1.2081402
ampl 2 0.3633486 -4.9269568 1.8590059
ampl 3 0.5535167 -0.1226836 8.0832378
ampl 4 1.0191972 0.5492328 8.6672361
ampl 5 0.3898570 -0.9982284 1.3397911
ampl 6 0.7838995 -3.0463937 1.6617784
ampl 7 -0.0845413 -1.6363487 0.6072287
ampl 8 0.3973571 -0.6370416 1.2374719
ampl 9 0.6634273 -3.2027133 0.8949979
ampl 10 -0.0176411 -1.8148987 1.2172768
rate 1 2.0748281 0.4417415 4.2679959
rate 2 1.6082447 0.3873175 3.2998544
rate 3 1.6276396 0.5214304 3.9119443
rate 4 1.6921749 0.6667222 3.5293919
rate 5 1.7416284 0.1741848 3.8666446
rate 6 1.6943163 0.4248442 3.4381555
rate 7 2.0284465 0.3386768 3.3150175
rate 8 1.2449642 0.5940424 4.1280756
rate 9 1.2655970 0.3138959 3.9372651
rate 10 2.2099059 0.6199482 4.0142250
asym 1 0.0181369 -0.2099792 0.1180537
asym 2 -0.0828012 -0.6863989 0.0134162
asym 3 0.0978848 -0.0523248 0.2242465
asym 4 0.1523575 -0.1210832 0.6341334
asym 5 -0.0275523 -0.3151141 0.0956343
asym 6 0.2297999 0.0586725 0.8026807
asym 7 0.0315837 -0.4508254 0.1161643
asym 8 0.0572859 -0.4507247 0.1584670
asym 9 0.0122410 -0.4304944 0.7584419
asym 10 0.1914040 0.0787843 0.5680915 
## LEARNING CURVES BASED ON MODEL ESTIMATES ##
dfAcc2$Participant<-dfAcc2$Participant
# learning curves per participant
dfAcc2 %>% 
  mutate(M_1 = predict(PP_1)$center) %>% 
  ggplot(aes(x = TimeRepX,
             y = M_1,
             color=Participant)) +
  facet_wrap(~Participant, scales = "free_y") +
  geom_smooth(se = F) +
  geom_point(alpha=0.2, size=1)+
  ylim(0,1.5)+
  labs(x="Trial",y="Model estimates")+
  theme_classic()+
  scale_color_manual(values=mycolors)
## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
## Warning: Removed 33221 rows containing non-finite values (stat_smooth).
## Warning: Removed 33221 rows containing missing values (geom_point).

# crossbar plots for each parameter and participant
P_1 %>% 
  re_scores() %>% 
  bayr::ranef() %>% 
  ggplot(aes(x = re_entity, 
             y = center, 
             ymin = lower, 
             ymax = upper)) +
  facet_grid(nonlin~1, scales = "free_y") +
  geom_crossbar(width = .2) +
  labs(x = "Subject", y = "Model estimates") +
  theme_classic()