Setting my environment- adding packhages-wd

knitr::opts_chunk$set(echo = TRUE)
library(lme4); library(ggplot2) ; library(lattice) ; library(tidyverse) ; library(readxl) ;library(readr) ; library(ggpubr) ; library(sjPlot)
## Caricamento del pacchetto richiesto: Matrix
## -- Attaching packages --------------------------------------- tidyverse 1.3.1 --
## v tibble  3.1.4     v dplyr   1.0.7
## v tidyr   1.1.4     v stringr 1.4.0
## v readr   2.0.2     v forcats 0.5.1
## v purrr   0.3.4
## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x tidyr::expand() masks Matrix::expand()
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
## x tidyr::pack()   masks Matrix::pack()
## x tidyr::unpack() masks Matrix::unpack()

Importing Data that are already cleaned using inter-quartile range of boxplots- in total only 8% of trials were deleted while one participant was excluded due to accuracy

data <- read_excel("Final_table_last.xlsx", 
    sheet = "clean")
accuracy_final <- read_delim("accuracy.csv", delim = ";", 
    escape_double = FALSE, trim_ws = TRUE)
## Rows: 12 Columns: 3
## -- Column specification --------------------------------------------------------
## Delimiter: ";"
## dbl (3): Participant, Accuracy, Sd
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
questions = read_excel("Questions.xlsx")
recognition = read_delim("Questions.csv", 
    delim = ";", escape_double = FALSE, trim_ws = TRUE)
## Rows: 2 Columns: 2
## -- Column specification --------------------------------------------------------
## Delimiter: ";"
## chr (1): Answer
## dbl (1): Responses
## 
## i Use `spec()` to retrieve the full column specification for this data.
## i Specify the column types or set `show_col_types = FALSE` to quiet this message.
sbjs = unique(data$Participant)
n = length(sbjs)
obs = rep(NA,n); obs = mapply(function(i)NROW(data[data$Participant==unique(data$Participant)[i],]),1:n) #observations for each subject

#----------------------------------------Descriptive stats--------------------------------------------#
tapply(data$RTs, data$Condition, summary)
## $Random
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   284.0   537.8   719.2   877.7  1146.8  2284.0 
## 
## $Sequence
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   263.2   462.3   638.9   763.1   987.9  1962.2
tapply(data$ACCURACY,data$Participant, mean)
##       S01       S02       S03       S05       S06       S07       S08       S09 
## 0.9635036 0.9928058 0.9571429 0.9694656 0.9635036 0.9705882 0.9858156 0.9411765 
##       S10       S11       S12 
## 1.0000000 0.9338235 0.9191176
data_agg = with(data,aggregate(RTs,list(Participant,Condition,Category),mean)); names(data_agg) = c("Participant","Condition","Category","RTS")

Making some plots

ggplot(data = data,aes(y=RTs, fill=Condition)) + geom_boxplot() + theme_bw() +facet_wrap(~Participant,scales = "free") + scale_fill_brewer(palette="Purples")

ggplot(data = data,aes(x=RTs,fill=Condition)) + geom_density(alpha=1) + theme_bw() + facet_wrap(~Participant,scales = "free") + scale_fill_brewer(palette="Purples")

ggplot(data = data,aes(x=Id,y=RTs,col=Condition)) + geom_point() + theme_bw() + facet_wrap(~Participant,scales = "free") + scale_color_brewer(palette="Purples") 

ggplot(accuracy_final, aes(y=Accuracy, x=Participant))  + 
  geom_line(aes(group = 1), data = accuracy_final) +
  geom_errorbar(aes(ymin = Accuracy-Sd, ymax = Accuracy+Sd),
    data = accuracy_final, width = 0.2)

p = ggline(data, x = "Participant", y = "ACCURACY", add = "mean_sd")
ggpar(p, ylim = c(0,1.2)) 

#---------------------aggregated plots-----------------------------------------------#
ggplot(data = data_agg,aes(x=RTS,fill=Condition)) + geom_density(alpha=0.5) + theme_bw() + scale_fill_brewer(palette="Purples")

ggplot(data = data_agg,aes(x=Condition,y=RTS)) + geom_boxplot(alpha=1) + theme_bw() + scale_color_brewer(palette="Purples")

#questions--------------------------
 ggbarplot(questions, "Answers", "Responses",
   fill = "Answers", color = "Answers",
   palette = c("#64A8D1", "#0969A2", "#9B59B6"))

ggbarplot(recognition, "Answer","Responses",
   fill = "Answer", color = "Answer",
   palette = c("#03436A", "#8E44AD"))

Running some models

data$Condition =as.factor(data$Condition)
contrasts(data$Condition) <- contr.sum(2)
mod0 = glmer(formula = I(RTs/1000)~1+(1|Participant), data = data,family = Gamma)
mod1 = glmer(formula = I(RTs/1000)~ Condition+(1+Condition|Participant),data = data,family = Gamma)

#model.matrix(~Condition,data=data)
#-----model comparison------------
anova(mod0,mod1)
## Data: data
## Models:
## mod0: I(RTs/1000) ~ 1 + (1 | Participant)
## mod1: I(RTs/1000) ~ Condition + (1 + Condition | Participant)
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)    
## mod0    3 413.48 429.44 -203.74   407.48                         
## mod1    6 286.20 318.12 -137.10   274.20 133.27  3  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
aics = AIC(mod0,mod1)$AIC
l_aics = exp(-0.5*(aics-min(aics)))
w_aics = l_aics/sum(l_aics) # probability of one model being in favor over the other
AIC_table = data.frame(paste0("mod",0:1),aics,w_aics); names(AIC_table) = c("model","AIC","relative AIC")
print(AIC_table)
##   model      AIC relative AIC
## 1  mod0 413.4758 2.305424e-28
## 2  mod1 286.2016 1.000000e+00
summary(mod1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: Gamma  ( inverse )
## Formula: I(RTs/1000) ~ Condition + (1 + Condition | Participant)
##    Data: data
## 
##      AIC      BIC   logLik deviance df.resid 
##    286.2    318.1   -137.1    274.2     1505 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.7873 -0.6878 -0.1971  0.4664  5.2510 
## 
## Random effects:
##  Groups      Name        Variance Std.Dev. Corr 
##  Participant (Intercept) 0.026543 0.16292       
##              Condition1  0.008276 0.09097  -0.16
##  Residual                0.145779 0.38181       
## Number of obs: 1511, groups:  Participant, 11
## 
## Fixed effects:
##             Estimate Std. Error t value Pr(>|z|)    
## (Intercept)  1.39998    0.10955  12.779   <2e-16 ***
## Condition1  -0.09326    0.05278  -1.767   0.0772 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##            (Intr)
## Condition1 -0.197
plot(effects::allEffects(mod1))

qqmath(ranef(mod1, condVar=TRUE))
## $Participant

plot_model(mod1, type = "eff", terms = c("Condition"))

data$residuals=residuals(mod1)
ggplot(data, aes(x=Condition, y=residuals,fill=Participant)) +   geom_boxplot() + theme_bw()

Trying to analyze over time

step <-  read_excel("Final_table_last.xlsx", 
    sheet = "Sheet3")
ggplot(step, aes(x=Step, y=Mean)) + 
    geom_errorbar(aes(ymin=Mean-Sd, ymax=Mean+Sd), width=.1) +
    geom_line() +
   scale_color_brewer(palette="Paired")+theme_bw() + facet_wrap(~Condition)

mod4= lm(formula = I(Mean/1000)~Condition*Step, data = step)
plot(effects::allEffects.default(mod4))

summary(mod4)
## 
## Call:
## lm(formula = I(Mean/1000) ~ Condition * Step, data = step)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.070350 -0.012808 -0.003688  0.027708  0.061552 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             0.98492    0.04227  23.303 2.35e-09 ***
## ConditionSequence      -0.30572    0.05709  -5.355 0.000459 ***
## Step                   -0.05083    0.01085  -4.684 0.001146 ** 
## ConditionSequence:Step  0.06036    0.01383   4.363 0.001815 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0454 on 9 degrees of freedom
## Multiple R-squared:  0.799,  Adjusted R-squared:  0.732 
## F-statistic: 11.93 on 3 and 9 DF,  p-value: 0.00173