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

dfAcc3<-read.table("Acc_Task3_1.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)

##USELESS Free learning curves (ROUGHT AND NOT ACCURATE)

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

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

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


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

#Proper

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


## FREE LEARNING CURVES ##

# Free learning curve 1
# per subject per sequence
dfAcc3 %>% 
  ggplot(aes(x = Time,
             y = DomHandX,
             color=Subject)) +
  geom_point() +
  geom_smooth(se = F) +
  facet_wrap(~Subject, 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

dfAcc3 %>% 
  ggplot(aes(x = Time,
             y = DomHandX,
             group = Subject,
             color= Subject)) +
  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
dfAcc3 %>% 
  ggplot(aes(x = TimeRepX,
             y = DomHandX,
             group = Subject,
             color=Subject)) +
  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(~Subject)
## `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(DomHandX ~ asym + ampl * exp(-rate * TimeRepX))

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

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 <- 
  dfAcc3 %>% 
  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("Subject"))
## 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[249] is inf, but must be finite!  (in 'modeld2dc59074fd6_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[92] is inf, but must be finite!  (in 'modeld2dc59074fd6_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 1: 
## Chain 1: Gradient evaluation took 0.268 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 2680 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: 3.266 seconds (Warm-up)
## Chain 1:                0.645 seconds (Sampling)
## Chain 1:                3.911 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL 'ce11eef883174df5125e8f267d0f1a59' NOW (CHAIN 2).
## Chain 2: Rejecting initial value:
## Chain 2:   Error evaluating the log probability at the initial value.
## Chain 2: Exception: exp_mod_normal_lpdf: Location parameter[8635] is inf, but must be finite!  (in 'modeld2dc59074fd6_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 2: Rejecting initial value:
## Chain 2:   Error evaluating the log probability at the initial value.
## Chain 2: Exception: exp_mod_normal_lpdf: Location parameter[16914] is inf, but must be finite!  (in 'modeld2dc59074fd6_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 2: 
## Chain 2: Gradient evaluation took 0.191 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 1910 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: 2.892 seconds (Warm-up)
## Chain 2:                0.698 seconds (Sampling)
## Chain 2:                3.59 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[31280] is inf, but must be finite!  (in 'modeld2dc59074fd6_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[1887] is inf, but must be finite!  (in 'modeld2dc59074fd6_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## 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:   Error evaluating the log probability at the initial value.
## Chain 3: Exception: exp_mod_normal_lpdf: Location parameter[66790] is inf, but must be finite!  (in 'modeld2dc59074fd6_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 3: 
## Chain 3: Gradient evaluation took 0.183 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 1830 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: 3.851 seconds (Warm-up)
## Chain 3:                0.806 seconds (Sampling)
## Chain 3:                4.657 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[181] is inf, but must be finite!  (in 'modeld2dc59074fd6_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[168654] is inf, but must be finite!  (in 'modeld2dc59074fd6_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## 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[2254] is inf, but must be finite!  (in 'modeld2dc59074fd6_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[11906] is inf, but must be finite!  (in 'modeld2dc59074fd6_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[7934] is inf, but must be finite!  (in 'modeld2dc59074fd6_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[41042] is inf, but must be finite!  (in 'modeld2dc59074fd6_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[15930] is inf, but must be finite!  (in 'modeld2dc59074fd6_ce11eef883174df5125e8f267d0f1a59' at line 85)
## 
## Chain 4: 
## Chain 4: Gradient evaluation took 0.206 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 2060 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.816 seconds (Warm-up)
## Chain 4:                0.67 seconds (Sampling)
## Chain 4:                3.486 seconds (Total)
## Chain 4:
## Warning: There were 8 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: 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
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 1.211874 0.1563928 1.515678
rate 2.219202 0.4245065 4.791077
asym 0.450876 0.3334083 1.690092 
bayr::grpef(P_1)
Coefficient estimates with 95% credibility limits
nonlin center lower upper
ampl 3.3174678 1.4809413 5.680950
rate 0.5223732 0.1536373 2.802884
asym 0.2368261 0.1477081 1.077902 
P_1 %>% re_scores() %>% bayr::ranef()
Coefficient estimates with 95% credibility limits
nonlin re_entity center lower upper
ampl 1 2.6358219 -9.0047866 8.3553644
ampl 2 5.7410657 2.3027403 8.6012503
ampl 3 5.7518058 -0.6253278 8.9939622
ampl 4 0.5794748 -5.8492460 5.1089596
ampl 5 -1.4035362 -1.9675785 3.0359894
ampl 6 -0.2396824 -9.2771927 2.1248745
ampl 7 -3.5691472 -9.0038300 0.8769815
ampl 8 4.7961620 -2.1315149 8.3950812
ampl 9 1.0247942 -3.5138951 5.8079898
ampl 10 1.1590246 0.4990026 9.1232928
rate 1 2.4964408 0.2027952 8.2132935
rate 2 2.8186105 0.1555028 5.6094428
rate 3 3.1861053 0.2673572 8.1641143
rate 4 2.9556873 0.4510080 6.3170517
rate 5 2.0896352 0.2029208 4.1969567
rate 6 1.9627975 0.2403868 5.0708120
rate 7 2.6208997 0.2548388 6.7182022
rate 8 1.3382137 0.5980535 4.2853503
rate 9 1.8667536 0.6335630 3.1384891
rate 10 0.8189017 0.2007836 4.8217236
asym 1 0.6512649 0.3417985 1.3790956
asym 2 0.5548112 0.1540614 2.3304532
asym 3 0.5453482 0.1498880 2.9254197
asym 4 0.4316749 -0.1887018 0.5468500
asym 5 0.5666214 0.2247566 1.9885948
asym 6 0.4927933 0.1754465 0.6485812
asym 7 0.5775389 0.1763669 3.1296315
asym 8 0.4868023 -0.1882663 0.7325478
asym 9 0.4719411 -0.3081539 0.8010359
asym 10 0.4247143 0.2343402 0.6106832 
## LEARNING CURVES BASED ON MODEL ESTIMATES ##
dfAcc3$Subject<-dfAcc3$Subject
# learning curves per Subject
dfAcc3 %>% 
  mutate(M_1 = predict(PP_1)$center) %>% 
  ggplot(aes(x = TimeRepX,
             y = M_1,
             color=Subject)) +
  facet_wrap(~Subject, 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 106260 rows containing non-finite values (stat_smooth).
## Warning: Removed 106260 rows containing missing values (geom_point).

# crossbar plots for each parameter and Subject
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()