Data Files & Libraries

Sys.setlocale("LC_ALL","Greek")
## Warning in Sys.setlocale("LC_ALL", "Greek"): using locale code page other than
## 65001 ("UTF-8") may cause problems
## [1] "LC_COLLATE=Greek_Greece.1253;LC_CTYPE=Greek_Greece.1253;LC_MONETARY=Greek_Greece.1253;LC_NUMERIC=C;LC_TIME=Greek_Greece.1253"
##libraries
library(dplyr)
library(sciplot)
library(ez)
library(rstatix)
## Warning: package 'rstatix' was built under R version 4.3.2
library(psychReport)
library(emmeans)
# Data are in two files, depending on tasks (counting vs. ordering).
################################################################################
# Counting Tasks
# in "cnd" variable, b1 denotes abstract sounds, b2 denotes environmental sounds
################################################################################
datab<-read.table("data_counting.txt", header=T)
#Convert chr to factor
datab <- mutate_if(datab, is.character, as.factor)
str(datab)
## 'data.frame':    5376 obs. of  12 variables:
##  $ blck : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ acc  : int  1 0 1 1 0 0 0 1 0 0 ...
##  $ trial: int  0 1 2 3 4 5 6 7 8 9 ...
##  $ resp : int  4 3 4 4 5 5 3 4 3 5 ...
##  $ rt   : int  8956 1581 4900 1902 3067 2352 1237 3027 2391 4708 ...
##  $ stim1: Factor w/ 8 levels "Bark","bird",..: 7 6 3 8 7 3 8 3 7 3 ...
##  $ stim2: Factor w/ 8 levels "Bark","bird",..: 6 3 7 6 3 8 7 8 8 6 ...
##  $ stim3: Factor w/ 8 levels "Bark","bird",..: 3 8 8 3 8 6 6 7 3 8 ...
##  $ stim4: Factor w/ 8 levels "Bark","bird",..: 8 7 6 7 6 7 3 6 6 7 ...
##  $ sbj  : Factor w/ 28 levels "sbj1","sbj10",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ cnd  : Factor w/ 2 levels "b1","b2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ music: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
#Calculate absolute error (ae), for supplementary analysis
datab$ae<-abs(datab$resp -4)
##################################################################################
# Ordering Tasks
# in "cnd" variable, c1 denotes abstract sounds, c2 denotes environmental sounds
##################################################################################
datac<-read.table("data_ordering.txt", header=T) 
#Convert chr to factor
datac <- mutate_if(datac, is.character, as.factor)
str(datac)
## 'data.frame':    5376 obs. of  32 variables:
##  $ blck                : int  1 1 1 1 1 1 1 1 1 1 ...
##  $ acc1                : int  1 0 0 1 1 0 0 0 0 0 ...
##  $ acc2                : int  1 1 1 0 0 0 0 1 0 0 ...
##  $ acc3                : int  1 0 0 0 0 1 1 0 0 0 ...
##  $ acc4                : int  1 1 0 0 1 1 0 0 0 0 ...
##  $ correct_button_order: Factor w/ 24 levels "['roi1', 'roi2', 'roi3', 'roi4']",..: 15 12 20 10 13 19 7 21 8 4 ...
##  $ trial               : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ resp1_roi           : Factor w/ 4 levels "roi1","roi2",..: 3 3 3 2 3 1 4 3 1 3 ...
##  $ resp2_roi           : Factor w/ 4 levels "roi1","roi2",..: 2 4 1 4 2 4 2 2 2 2 ...
##  $ resp3_roi           : Factor w/ 4 levels "roi1","roi2",..: 1 2 2 1 1 2 3 4 3 1 ...
##  $ resp4_roi           : Factor w/ 4 levels "roi1","roi2",..: 4 1 4 3 4 3 1 1 4 4 ...
##  $ labels              : Factor w/ 1290 levels "['αδιέξοδο', 'επαφή', 'επικοινωνία', 'παράσιτα']",..: 358 353 370 677 673 369 356 355 671 369 ...
##  $ position            : Factor w/ 24 levels "[0, 1, 2, 3]",..: 15 20 12 19 9 10 7 16 8 5 ...
##  $ resp1               : Factor w/ 8 levels "bark","bird",..: 8 8 8 8 8 6 6 6 7 3 ...
##  $ resp2               : Factor w/ 8 levels "bark","bird",..: 7 7 6 6 6 8 8 8 3 7 ...
##  $ resp3               : Factor w/ 8 levels "bark","bird",..: 3 6 7 7 7 7 7 7 8 6 ...
##  $ resp4               : Factor w/ 8 levels "bark","bird",..: 6 3 3 3 3 3 3 3 6 8 ...
##  $ rt1                 : int  11774 5173 10225 12120 13025 5189 8326 11586 4399 6527 ...
##  $ rt2                 : int  511 480 871 856 783 696 560 351 1583 624 ...
##  $ rt3                 : int  703 495 471 1056 464 599 1646 976 1455 783 ...
##  $ rt4                 : int  528 392 840 752 879 336 640 718 423 728 ...
##  $ stim1               : Factor w/ 8 levels "Bark","bird",..: 8 6 3 8 8 8 8 7 3 6 ...
##  $ stim2               : Factor w/ 8 levels "Bark","bird",..: 7 7 6 3 7 6 3 8 7 3 ...
##  $ stim3               : Factor w/ 8 levels "Bark","bird",..: 3 8 8 6 6 7 7 3 6 8 ...
##  $ stim4               : Factor w/ 8 levels "Bark","bird",..: 6 3 7 7 3 3 6 6 8 7 ...
##  $ sbj                 : Factor w/ 28 levels "sbj1","sbj10",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ cnd                 : Factor w/ 2 levels "c1","c2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ acc                 : int  1 0 0 0 0 0 0 0 0 0 ...
##  $ criterion           : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ music               : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ strtg_abs           : Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
##  $ strtg_env           : Factor w/ 2 levels "no","yes": 2 2 2 2 2 2 2 2 2 2 ...

Descriptive Statistics

Counting Tasks

################
# Counting Tasks
################
#calculate average accuracy by participant and condition
datab_sbj<-aggregate(datab$acc, list(datab$sbj, datab$cnd), mean)
colnames(datab_sbj)<-c("sbj","cnd","acc") 
str(datab_sbj)
## 'data.frame':    56 obs. of  3 variables:
##  $ sbj: Factor w/ 28 levels "sbj1","sbj10",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ cnd: Factor w/ 2 levels "b1","b2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ acc: num  0.375 0.2188 0.125 0.0104 0.0312 ...
########abstract sounds
mean(datab_sbj[datab_sbj$cnd=="b1",]$acc)
## [1] 0.3772321
sd(datab_sbj[datab_sbj$cnd=="b1",]$acc)
## [1] 0.2779854
########environmental sounds
mean(datab_sbj[datab_sbj$cnd=="b2",]$acc)
## [1] 0.6796875
sd(datab_sbj[datab_sbj$cnd=="b2",]$acc)
## [1] 0.3154774
#boxplot
col=c("grey75","grey25")
plot(acc~cnd, data=datab_sbj, main= "Counting Tasks, Accuracy", xlab="", ylab="Proportion Correct", frame.plot=F, xaxt="n", boxwex=0.5, las=2, col=col, range=0, lty=1)
axis(side=1, tick=F, at=c(1,2),labels=c("Abstract sounds", "Environmental sounds" ))

#range = 0 means that error bars extend to the full range

Ordering Tasks

#calculate average accuracy by participant and condition
datac_sbj<-aggregate(datac$acc, list(datac$sbj, datac$cnd), mean)
colnames(datac_sbj)<-c("sbj","cnd","acc") 
str(datac_sbj)
## 'data.frame':    56 obs. of  3 variables:
##  $ sbj: Factor w/ 28 levels "sbj1","sbj10",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ cnd: Factor w/ 2 levels "c1","c2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ acc: num  0.302 0.115 0.115 0.323 0.125 ...
########abstract sounds
mean(datac_sbj[datac_sbj$cnd=="c1",]$acc)
## [1] 0.234003
sd(datac_sbj[datac_sbj$cnd=="c1",]$acc)
## [1] 0.1923557
########environmental sounds
mean(datac_sbj[datac_sbj$cnd=="c2",]$acc)
## [1] 0.2689732
sd(datac_sbj[datac_sbj$cnd=="c2",]$acc)
## [1] 0.1063072
#plot accuracy by cnd
col=c("grey75","grey25")
plot(acc~cnd, data=datac_sbj, main="Ordering Tasks, Accuracy", ylim=c(0.0, 1.0), xlab="", ylab="Proportion Correct", frame.plot=F, xaxt="n", boxwex=0.5, las=2, col=col, range=0,lty=5)
axis(side=1, at=c(1,2),tick=F, labels=c("Abstract sounds", "Environmental sounds" ))

Inferential Statistics

Countinig Tasks

Main Accuracy Analyses
#check if performance in each condition is chance performance
#participants were instructed to  press keys from 1 to 9. thus, chance performance is 1/9
t.test(datab_sbj[datab_sbj$cnd=="b1",]$acc, mu=1/9, alternative = 'greater')
## 
##  One Sample t-test
## 
## data:  datab_sbj[datab_sbj$cnd == "b1", ]$acc
## t = 5.0657, df = 27, p-value = 1.277e-05
## alternative hypothesis: true mean is greater than 0.1111111
## 95 percent confidence interval:
##  0.2877511       Inf
## sample estimates:
## mean of x 
## 0.3772321
t.test(datab_sbj[datab_sbj$cnd=="b2",]$acc, mu=1/9, alternative = 'greater')
## 
##  One Sample t-test
## 
## data:  datab_sbj[datab_sbj$cnd == "b2", ]$acc
## t = 9.5367, df = 27, p-value = 1.945e-10
## alternative hypothesis: true mean is greater than 0.1111111
## 95 percent confidence interval:
##  0.5781381       Inf
## sample estimates:
## mean of x 
## 0.6796875
#check if participants' accuracy differed between sound categories
t.test(acc~cnd, data = datab_sbj, paired=T)
## 
##  Paired t-test
## 
## data:  acc by cnd
## t = -4.7449, df = 27, p-value = 6.045e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.4332444 -0.1716663
## sample estimates:
## mean difference 
##      -0.3024554
#there's a difference
#calculate effect size (Cohen's d)
datab_sbj %>% cohens_d(acc~cnd)
## # A tibble: 1 x 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 acc   b1     b2       -1.02    28    28 large
Order of Sounds
########abstract sounds
dat_b1_stim<-droplevels(datab[datab$cnd=="b1",])
#create new factor
dat_b1_stim$stim_code<-as.factor(paste(substring(dat_b1_stim$stim1,1,1),substring(dat_b1_stim$stim2,1,1),substring(dat_b1_stim$stim3,1,1),substring(dat_b1_stim$stim4,1,1), sep=""))
levels(dat_b1_stim$stim_code)
##  [1] "bhHL" "bHhL" "bhLH" "bHLh" "bLhH" "bLHh" "hbHL" "HbhL" "hbLH" "HbLh"
## [11] "hHbL" "HhbL" "hHLb" "HhLb" "hLbH" "HLbh" "hLHb" "HLhb" "LbhH" "LbHh"
## [21] "LhbH" "LHbh" "LhHb" "LHhb"
#h=hiss, b=buzz, L=low, H= high

#calculate average accuracy per sbj and specific order of stimuli
datb1_stim<-aggregate(dat_b1_stim$acc, list(dat_b1_stim$stim_code, dat_b1_stim$sbj), mean)
colnames(datb1_stim)<-c("stim_code", "sbj","acc")
ezANOVA(data=datb1_stim, dv=acc, wid=sbj, within=stim_code, type=3)
## $ANOVA
##      Effect DFn DFd        F          p p<.05        ges
## 2 stim_code  23 621 1.696197 0.02266271     * 0.01996842
## 
## $`Mauchly's Test for Sphericity`
##      Effect            W           p p<.05
## 2 stim_code 4.819655e-09 0.001010348     *
## 
## $`Sphericity Corrections`
##      Effect       GGe      p[GG] p[GG]<.05       HFe      p[HF] p[HF]<.05
## 2 stim_code 0.4782541 0.07344036           0.8315019 0.03277858         *
#Sphericity is not fulfilled, thus we check GG "Sphericity Corrections"
#p=.074, no effect
#because we use the GG correction, we need to multiply dfs by GGe
########environmental sounds
dat_b2_stim<-droplevels(datab[datab$cnd=="b2",])

#create new factor
dat_b2_stim$stim_code<-as.factor(paste(substring(dat_b2_stim$stim1,1,1),substring(dat_b2_stim$stim2,1,1),substring(dat_b2_stim$stim3,1,1),substring(dat_b2_stim$stim4,1,1), sep=""))
levels(dat_b2_stim$stim_code)
##  [1] "bBdg" "Bbdg" "bBgd" "Bbgd" "bdBg" "Bdbg" "bdgB" "Bdgb" "bgBd" "Bgbd"
## [11] "bgdB" "Bgdb" "dbBg" "dBbg" "dbgB" "dBgb" "dgbB" "dgBb" "gbBd" "gBbd"
## [21] "gbdB" "gBdb" "gdbB" "gdBb"
#B=bark, b=bird, d=dog, g=glass breaking

#calculate average accuracy per sbj and specific order of stimuli
datb2_stim<-aggregate(dat_b2_stim$acc, list(dat_b2_stim$stim_code, dat_b2_stim$sbj), mean)
colnames(datb2_stim)<-c("stim_code", "sbj","acc")
ezANOVA(data=datb2_stim, dv=acc, wid=sbj, within=stim_code, type=3)
## $ANOVA
##      Effect DFn DFd        F         p p<.05        ges
## 2 stim_code  23 621 1.310722 0.1513678       0.01058807
## 
## $`Mauchly's Test for Sphericity`
##      Effect            W            p p<.05
## 2 stim_code 6.422529e-12 1.694997e-13     *
## 
## $`Sphericity Corrections`
##      Effect       GGe    p[GG] p[GG]<.05       HFe     p[HF] p[HF]<.05
## 2 stim_code 0.4214435 0.226612           0.6768128 0.1884865
#p=.22 no effect
Musical Training
########abstract sounds
#calculate average accuracy per sbj and musical training
datb1_m<-aggregate(dat_b1_stim$acc, list(dat_b1_stim$music, dat_b1_stim$sbj), mean)
colnames(datb1_m)<-c("music", "sbj","acc")
#some musical training
mean(datb1_m[datb1_m$music=="yes",]$acc)
## [1] 0.5464015
sd(datb1_m[datb1_m$music=="yes",]$acc)
## [1] 0.3215483
#no musical training
mean(datb1_m[datb1_m$music=="no",]$acc)
## [1] 0.2677696
sd(datb1_m[datb1_m$music=="no",]$acc)
## [1] 0.1826917
#check if there's a difference between groups
var.test(acc~music, data=datb1_m)
## 
##  F test to compare two variances
## 
## data:  acc by music
## F = 0.32281, num df = 16, denom df = 10, p-value = 0.04299
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.09232941 0.96395940
## sample estimates:
## ratio of variances 
##          0.3228087
# p = .0042, equal variances cannot be assumed, Welch's t-test
t.test(acc~music,paired=F, data=datb1_m)
## 
##  Welch Two Sample t-test
## 
## data:  acc by music
## t = -2.6139, df = 14.226, p-value = 0.02021
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.50691740 -0.05034642
## sample estimates:
##  mean in group no mean in group yes 
##         0.2677696         0.5464015
datb1_m %>% cohens_d(acc~music, paired=F)
## # A tibble: 1 x 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 acc   no     yes      -1.07    17    11 large
#there's a difference
#functions to draw hatched boxplots
#modified from https://stackoverflow.com/questions/17423115/colorfill-boxplot-in-r-cran-with-lines-dots-or-similar
boxpattern_abstract <- function(y, xcenter, boxwidth, angle=NULL, angle.density=10, ...) {
    # draw an individual box
    bstats <- boxplot.stats(y)
    bxmin <- bstats$stats[1]
    bxq2 <- bstats$stats[2]
    bxmedian <- bstats$stats[3]
    bxq4 <- bstats$stats[4]
    bxmax <- bstats$stats[5]
    bleft <- xcenter-(boxwidth/2)
    bright <- xcenter+(boxwidth/2)
    bleft_seg <- xcenter-(boxwidth/4)
    bright_seg <- xcenter+(boxwidth/4)
    # boxplot
    polygon(c(bleft,bright,bright,bleft,bleft),
            c(bxq2,bxq2,bxq4,bxq4,bxq2), angle=angle[1], col="grey75")
    polygon(c(bleft,bright,bright,bleft,bleft),
            c(bxq2,bxq2,bxq4,bxq4,bxq2), angle=angle[2], density=angle.density)
    # lines
    segments(bleft,bxmedian,bright,bxmedian,lwd=3) # median
    segments(bleft_seg,bxmin,bright_seg,bxmin,lwd=1) # min
    segments(xcenter,bxmin,xcenter,bxq2,lwd=1)
    segments(bleft_seg,bxmax,bright_seg,bxmax,lwd=1) # max
    segments(xcenter,bxq4,xcenter,bxmax,lwd=1)
  }

drawboxplots_abstract <- function(y, x, boxwidth=1, angle=NULL, ...){
  # figure out all the boxes and start the plot
  groups <- split(y,as.factor(x))
  len <- length(groups)
  bxylim <- c((min(y)-0.04*abs(min(y))),(max(y)+0.04*max(y)))
  xcenters <- seq(1,max(2,(len*(1.4))),length.out=len)
  if(is.null(angle)){
    angle <- seq(-90,35,length.out=len)
    angle <- lapply(angle,function(x) c(x,x))
  }
  else if(!length(angle)==len)
    stop("angle must be a vector or list of two-element vectors")
  else if(!is.list(angle))
    angle <- lapply(angle,function(x) c(x,x))
  # draw plot area
  plot(10, xlim=c(.97*(min(xcenters)-1), 1.04*(max(xcenters)+1)),
       ylim=bxylim, 
       xlab="",
       ylab="Proportion Correct", 
       las=1, 
       frame.plot=F, 
       xaxt="n", 
       main="Counting Task, Abstract Sounds\nAccuracy by Musical Training")
  
  axis(1, at=xcenters, tick=F, labels=c("None\nN=11", "Some\nN=17"))
  # draw boxplots
  plots <- mapply(boxpattern_abstract, y=groups, xcenter=xcenters,
                  boxwidth=boxwidth, angle=angle, ...)
}

drawboxplots_abstract(datb1_m$acc, datb1_m$music)

########environmental sounds
#calculate average accuracy per sbj and musical training
datb2_m<-aggregate(dat_b2_stim$acc, list(dat_b2_stim$music, dat_b2_stim$sbj), mean)
colnames(datb2_m)<-c("music", "sbj","acc")
#some musical training
mean(datb2_m[datb2_m$music=="yes",]$acc)
## [1] 0.7613636
sd(datb2_m[datb2_m$music=="yes",]$acc)
## [1] 0.2809974
#no musical training
mean(datb2_m[datb2_m$music=="no",]$acc)
## [1] 0.6268382
sd(datb2_m[datb2_m$music=="no",]$acc)
## [1] 0.3332364
#check if there's a difference between groups
var.test(acc~music, data=datb2_m)
## 
##  F test to compare two variances
## 
## data:  acc by music
## F = 1.4064, num df = 16, denom df = 10, p-value = 0.5937
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4022492 4.1996575
## sample estimates:
## ratio of variances 
##           1.406372
# p = 0.594, equal variances may be assumed
t.test(acc~music,paired=F, var.equal=T, data=datb2_m)
## 
##  Two Sample t-test
## 
## data:  acc by music
## t = -1.1066, df = 26, p-value = 0.2786
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.3844173  0.1153665
## sample estimates:
##  mean in group no mean in group yes 
##         0.6268382         0.7613636
#no difference
Supplementary Analyses of Absolute Error
###################################
#Supplementary Analysis.
#Measure: Absolute Error = abs(response-correct)
###################################


#aggregate by sbj
db_sbj<-aggregate(datab$ae, list(datab$sbj, datab$cnd), mean)
colnames(db_sbj)<-c("sbj","cnd", "ae")
#descriptives
#abstract
mean(db_sbj[db_sbj$cnd=="b1",]$ae)
## [1] 0.7321429
sd(db_sbj[db_sbj$cnd=="b1",]$ae)
## [1] 0.3949273
#environmental
mean(db_sbj[db_sbj$cnd=="b2",]$ae)
## [1] 0.3333333
sd(db_sbj[db_sbj$cnd=="b2",]$ae)
## [1] 0.3399942
#plot ae by cnd
col=c("grey75","grey25")
plot(ae~cnd, data=db_sbj, main="Enumeration Task, DV= Absolute Error", xlab="", ylab="Mean Absolute Error", frame.plot=F, xaxt="n", boxwex=0.5, las=2, col=col, range=0, ylim=c(0,2))
axis(side=1, at=c(1,2),tick=F, labels=c("Abstract sounds", "Environmental sounds" ))

#check if participants' ae was differed between sound categories
t.test(ae~cnd, data = db_sbj, paired=T)
## 
##  Paired t-test
## 
## data:  ae by cnd
## t = 4.7381, df = 27, p-value = 6.158e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.2261047 0.5715143
## sample estimates:
## mean difference 
##       0.3988095
#there's a difference
#calculate effect size (Cohen's d)
db_sbj %>% cohens_d(ae~cnd)
## # A tibble: 1 x 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 ae    b1     b2        1.08    28    28 large
#there's a difference
#############################################
#check if sequence of sounds affects ae
#############################################

########abstract sounds

#calculate average ad per sbj and specific order of stimuli
db1_stim<-aggregate(dat_b1_stim$ae, list(dat_b1_stim$stim_code, dat_b1_stim$sbj), mean)
colnames(db1_stim)<-c("stim_code", "sbj","ae")

ezANOVA(data=db1_stim, dv=ae, wid=sbj, within=stim_code, type=3)
## $ANOVA
##      Effect DFn DFd        F            p p<.05        ges
## 2 stim_code  23 621 2.673234 4.352252e-05     * 0.02535864
## 
## $`Mauchly's Test for Sphericity`
##      Effect            W           p p<.05
## 2 stim_code 1.072643e-08 0.006135115     *
## 
## $`Sphericity Corrections`
##      Effect       GGe       p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 stim_code 0.4832871 0.002641338         * 0.8464332 0.0001450472         *
#Sphericity is not fulfilled, thus we check GG "Sphericity Corrections"
#because we use the GG correction, we need to multiply dfs by GGe
#p = 0.003, there is an effect

#there is an effect, the command below reports partial eta squared (in line 2)
aovEffectSize(ezANOVA(data=db1_stim, dv=ae, wid=sbj, within=stim_code, type=3, return_aov = TRUE, detailed = TRUE))
## $ANOVA
##        Effect DFn DFd        SSn       SSd         F            p p<.05
## 1 (Intercept)   1  27 360.214286 101.06696 96.231106 2.140498e-10     *
## 2   stim_code  23 621   3.566964  36.02679  2.673234 4.352252e-05     *
##          pes
## 1 0.78089947
## 2 0.09008907
## 
## $`Mauchly's Test for Sphericity`
##      Effect            W           p p<.05
## 2 stim_code 1.072643e-08 0.006135115     *
## 
## $`Sphericity Corrections`
##      Effect       GGe       p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 stim_code 0.4832871 0.002641338         * 0.8464332 0.0001450472         *
## 
## $aov
## 
## Call:
## aov(formula = formula(aov_formula), data = data)
## 
## Grand Mean: 0.7321429
## 
## Stratum 1: sbj
## 
## Terms:
##                 Residuals
## Sum of Squares    101.067
## Deg. of Freedom        27
## 
## Residual standard error: 1.934741
## 
## Stratum 2: sbj:stim_code
## 
## Terms:
##                 stim_code Residuals
## Sum of Squares    3.56696  36.02679
## Deg. of Freedom        23       621
## 
## Residual standard error: 0.2408613
## Estimated effects may be unbalanced
ezPlot(data=db1_stim, dv=ae, wid=sbj, within=stim_code, type=3, x=stim_code)

#pairwise t-test, bonferroni corrected
pairwise.t.test(db1_stim$ae, db1_stim$stim_code, paired=T, p.adjust.method= "bonferroni")
## 
##  Pairwise comparisons using paired t tests 
## 
## data:  db1_stim$ae and db1_stim$stim_code 
## 
##      bhHL  bHhL  bhLH  bHLh  bLhH  bLHh  hbHL  HbhL  hbLH  HbLh  hHbL  HhbL 
## bHhL 1.000 -     -     -     -     -     -     -     -     -     -     -    
## bhLH 1.000 1.000 -     -     -     -     -     -     -     -     -     -    
## bHLh 1.000 1.000 1.000 -     -     -     -     -     -     -     -     -    
## bLhH 1.000 1.000 1.000 1.000 -     -     -     -     -     -     -     -    
## bLHh 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -     -     -    
## hbHL 0.179 0.048 1.000 1.000 1.000 0.829 -     -     -     -     -     -    
## HbhL 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -    
## hbLH 0.015 0.027 1.000 1.000 1.000 0.541 1.000 1.000 -     -     -     -    
## HbLh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -    
## hHbL 1.000 0.381 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -    
## HhbL 1.000 1.000 1.000 1.000 1.000 1.000 0.466 1.000 0.370 1.000 1.000 -    
## hHLb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## HhLb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## hLbH 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## HLbh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## hLHb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## HLhb 1.000 1.000 1.000 1.000 1.000 1.000 0.215 1.000 0.029 1.000 1.000 1.000
## LbhH 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LbHh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LhbH 1.000 1.000 1.000 1.000 1.000 1.000 0.012 1.000 0.106 1.000 1.000 1.000
## LHbh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LhHb 1.000 0.203 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LHhb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
##      hHLb  HhLb  hLbH  HLbh  hLHb  HLhb  LbhH  LbHh  LhbH  LHbh  LhHb 
## bHhL -     -     -     -     -     -     -     -     -     -     -    
## bhLH -     -     -     -     -     -     -     -     -     -     -    
## bHLh -     -     -     -     -     -     -     -     -     -     -    
## bLhH -     -     -     -     -     -     -     -     -     -     -    
## bLHh -     -     -     -     -     -     -     -     -     -     -    
## hbHL -     -     -     -     -     -     -     -     -     -     -    
## HbhL -     -     -     -     -     -     -     -     -     -     -    
## hbLH -     -     -     -     -     -     -     -     -     -     -    
## HbLh -     -     -     -     -     -     -     -     -     -     -    
## hHbL -     -     -     -     -     -     -     -     -     -     -    
## HhbL -     -     -     -     -     -     -     -     -     -     -    
## hHLb -     -     -     -     -     -     -     -     -     -     -    
## HhLb 1.000 -     -     -     -     -     -     -     -     -     -    
## hLbH 1.000 1.000 -     -     -     -     -     -     -     -     -    
## HLbh 1.000 1.000 1.000 -     -     -     -     -     -     -     -    
## hLHb 1.000 1.000 1.000 1.000 -     -     -     -     -     -     -    
## HLhb 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -     -    
## LbhH 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -    
## LbHh 1.000 1.000 1.000 1.000 1.000 0.381 1.000 -     -     -     -    
## LhbH 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -    
## LHbh 1.000 1.000 1.000 1.000 1.000 0.370 1.000 1.000 1.000 -     -    
## LhHb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -    
## LHhb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## 
## P value adjustment method: bonferroni
#differences between:
#bhHL - hbLH (0.015)
#bHhL - hbHL (0.048)
#bHhL - hbLH (0.027)
#hbHL - LhbH (0.012)
########environmental sounds

#calculate average accuracy per sbj and specific order of stimuli
db2_stim<-aggregate(dat_b2_stim$ae, list(dat_b2_stim$stim_code, dat_b2_stim$sbj), mean)
colnames(db2_stim)<-c("stim_code", "sbj","ae")

ezANOVA(data=db2_stim, dv=ae, wid=sbj, within=stim_code, type=3)
## $ANOVA
##      Effect DFn DFd        F         p p<.05        ges
## 2 stim_code  23 621 1.677874 0.0250519     * 0.01399553
## 
## $`Mauchly's Test for Sphericity`
##      Effect            W            p p<.05
## 2 stim_code 1.749636e-12 5.894438e-16     *
## 
## $`Sphericity Corrections`
##      Effect       GGe     p[GG] p[GG]<.05       HFe      p[HF] p[HF]<.05
## 2 stim_code 0.3934248 0.0944021           0.6088128 0.05789012
#Sphericity is not fulfilled, thus we check GG "Sphericity Corrections"
#because we use the GG correction, we need to multiply dfs by GGe
#p=.094 no effect
###########################################
#check if musical training affects AE
###########################################

########abstract sounds
#calculate average AE per sbj and musical training

#calculate average accuracy per sbj and musical training
db1_m<-aggregate(dat_b1_stim$ae, list(dat_b1_stim$music, dat_b1_stim$sbj), mean)
colnames(db1_m)<-c("music", "sbj","ae")

#some musical training
mean(db1_m[db1_m$music=="yes",]$ae)
## [1] 0.4668561
sd(db1_m[db1_m$music=="yes",]$ae)
## [1] 0.3362121
#no musical training
mean(db1_m[db1_m$music=="no",]$ae)
## [1] 0.903799
sd(db1_m[db1_m$music=="no",]$ae)
## [1] 0.3359381
#check if there's a difference between groups
var.test(ae~music, data=db1_m)
## 
##  F test to compare two variances
## 
## data:  ae by music
## F = 0.99837, num df = 16, denom df = 10, p-value = 0.9624
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.2855529 2.9812975
## sample estimates:
## ratio of variances 
##          0.9983706
# p = 0.96, equal variances can be assumed
t.test(ae~music, paired=F, var.equal=T, data=db1_m)
## 
##  Two Sample t-test
## 
## data:  ae by music
## t = 3.3602, df = 26, p-value = 0.002415
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  0.1696562 0.7042297
## sample estimates:
##  mean in group no mean in group yes 
##         0.9037990         0.4668561
db1_m %>% cohens_d(ae~music, paired=F)
## # A tibble: 1 x 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 ae    no     yes       1.30    17    11 large
#there's a difference
########environmental sounds
#calculate average AE per sbj and musical training

#calculate average accuracy per sbj and musical training
db2_m<-aggregate(dat_b2_stim$ae, list(dat_b2_stim$music, dat_b2_stim$sbj), mean)
colnames(db2_m)<-c("music", "sbj","ae")

#some musical training
mean(db2_m[db2_m$music=="yes",]$ae)
## [1] 0.2395833
sd(db2_m[db2_m$music=="yes",]$ae)
## [1] 0.2816741
#no musical training
mean(db2_m[db2_m$music=="no",]$ae)
## [1] 0.3939951
sd(db2_m[db2_m$music=="no",]$ae)
## [1] 0.3681418
#check if there's a difference between groups
var.test(ae~music, data=db2_m)
## 
##  F test to compare two variances
## 
## data:  ae by music
## F = 1.7082, num df = 16, denom df = 10, p-value = 0.3934
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.4885751 5.1009377
## sample estimates:
## ratio of variances 
##           1.708191
# p = 0.39, equal variances can be assumed
t.test(ae~music, paired=F, var.equal=T, data=db2_m)
## 
##  Two Sample t-test
## 
## data:  ae by music
## t = 1.1823, df = 26, p-value = 0.2478
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.1140468  0.4228703
## sample estimates:
##  mean in group no mean in group yes 
##         0.3939951         0.2395833
db2_m %>% cohens_d(ae~music, paired=F)
## # A tibble: 1 x 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 ae    no     yes      0.471    17    11 small
#there's no difference

Ordering Tasks

Main Accuracy Analyses
#Check if performance is chance performance
#there are 24 permutations of the four keys. thus, chance performance is 1/24
t.test(datac_sbj[datac_sbj$cnd=="c1",]$acc, mu=1/24, alternative = 'greater')
## 
##  One Sample t-test
## 
## data:  datac_sbj[datac_sbj$cnd == "c1", ]$acc
## t = 5.291, df = 27, p-value = 6.974e-06
## alternative hypothesis: true mean is greater than 0.04166667
## 95 percent confidence interval:
##  0.1720853       Inf
## sample estimates:
## mean of x 
##  0.234003
t.test(datac_sbj[datac_sbj$cnd=="c2",]$acc, mu=1/24, alternative = 'greater')
## 
##  One Sample t-test
## 
## data:  datac_sbj[datac_sbj$cnd == "c2", ]$acc
## t = 11.314, df = 27, p-value = 4.711e-12
## alternative hypothesis: true mean is greater than 0.04166667
## 95 percent confidence interval:
##  0.2347538       Inf
## sample estimates:
## mean of x 
## 0.2689732
#check if participants' accuracy differed between sound categories
t.test(acc~cnd, data = datac_sbj, paired=T)
## 
##  Paired t-test
## 
## data:  acc by cnd
## t = -1.5409, df = 27, p-value = 0.135
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.08153544  0.01159497
## sample estimates:
## mean difference 
##     -0.03497024
#there's a difference
#calculate effect size (Cohen's d)
datac_sbj %>% cohens_d(acc~cnd)
## # A tibble: 1 x 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 acc   c1     c2      -0.225    28    28 small
Order of Sounds
########abstract sounds
dat_c1_stim<-droplevels(datac[datac$cnd=="c1",])

#create new factor
dat_c1_stim$stim_code<-as.factor(paste(substring(dat_c1_stim$stim1,1,1),substring(dat_c1_stim$stim2,1,1),substring(dat_c1_stim$stim3,1,1),substring(dat_c1_stim$stim4,1,1), sep=""))
levels(dat_c1_stim$stim_code)
##  [1] "bhHL" "bHhL" "bhLH" "bHLh" "bLhH" "bLHh" "hbHL" "HbhL" "hbLH" "HbLh"
## [11] "hHbL" "HhbL" "hHLb" "HhLb" "hLbH" "HLbh" "hLHb" "HLhb" "LbhH" "LbHh"
## [21] "LhbH" "LHbh" "LhHb" "LHhb"
#h=hiss, b=buzz, L=low, H= high

#aggregate data
datc1_stim<-aggregate(dat_c1_stim$acc, list(dat_c1_stim$stim_code, dat_c1_stim$sbj), mean)
colnames(datc1_stim)<-c("stim_code", "sbj","acc")
#check if there's an effect 
ezANOVA(data=datc1_stim, dv=acc, wid=sbj, within=stim_code, type=3)
## $ANOVA
##      Effect DFn DFd        F           p p<.05        ges
## 2 stim_code  23 621 3.023819 3.47362e-06     * 0.05750872
## 
## $`Mauchly's Test for Sphericity`
##      Effect           W         p p<.05
## 2 stim_code 7.37825e-08 0.1637892      
## 
## $`Sphericity Corrections`
##      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 stim_code 0.4790159 0.0007568019         * 0.8337478 1.899969e-05         *
#there is an effect, the command below reports partial eta squared (in line 2)
aovEffectSize(ezANOVA(data=datc1_stim, dv=acc, wid=sbj, within=stim_code, type=3, return_aov = TRUE, detailed = TRUE))
## $ANOVA
##        Effect DFn DFd       SSn      SSd         F            p p<.05       pes
## 1 (Intercept)   1  27 36.796968 23.97647 41.437216 6.743071e-07     * 0.6054778
## 2   stim_code  23 621  3.214193 28.69987  3.023819 3.473620e-06     * 0.1007140
## 
## $`Mauchly's Test for Sphericity`
##      Effect           W         p p<.05
## 2 stim_code 7.37825e-08 0.1637892      
## 
## $`Sphericity Corrections`
##      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 stim_code 0.4790159 0.0007568019         * 0.8337478 1.899969e-05         *
## 
## $aov
## 
## Call:
## aov(formula = formula(aov_formula), data = data)
## 
## Grand Mean: 0.234003
## 
## Stratum 1: sbj
## 
## Terms:
##                 Residuals
## Sum of Squares   23.97647
## Deg. of Freedom        27
## 
## Residual standard error: 0.9423467
## 
## Stratum 2: sbj:stim_code
## 
## Terms:
##                 stim_code Residuals
## Sum of Squares   3.214193 28.699870
## Deg. of Freedom        23       621
## 
## Residual standard error: 0.2149781
## Estimated effects may be unbalanced
ezPlot(data=datc1_stim, dv=acc, wid=sbj, within=stim_code, type=3, x=stim_code)

#pairwise t-test, bonferroni corrected
pairwise.t.test(datc1_stim$acc, datc1_stim$stim_code, paired=T, p.adjust.method= "bonferroni")
## 
##  Pairwise comparisons using paired t tests 
## 
## data:  datc1_stim$acc and datc1_stim$stim_code 
## 
##      bhHL  bHhL  bhLH  bHLh  bLhH  bLHh  hbHL  HbhL  hbLH  HbLh  hHbL  HhbL 
## bHhL 0.293 -     -     -     -     -     -     -     -     -     -     -    
## bhLH 1.000 1.000 -     -     -     -     -     -     -     -     -     -    
## bHLh 1.000 1.000 1.000 -     -     -     -     -     -     -     -     -    
## bLhH 1.000 1.000 1.000 1.000 -     -     -     -     -     -     -     -    
## bLHh 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -     -     -    
## hbHL 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -     -    
## HbhL 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -    
## hbLH 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -     -    
## HbLh 0.767 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -    
## hHbL 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -    
## HhbL 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -    
## hHLb 1.000 0.071 1.000 1.000 1.000 1.000 1.000 1.000 1.000 0.474 1.000 1.000
## HhLb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## hLbH 0.604 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## HLbh 1.000 0.430 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## hLHb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## HLhb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LbhH 0.582 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LbHh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LhbH 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LHbh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LhHb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## LHhb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
##      hHLb  HhLb  hLbH  HLbh  hLHb  HLhb  LbhH  LbHh  LhbH  LHbh  LhHb 
## bHhL -     -     -     -     -     -     -     -     -     -     -    
## bhLH -     -     -     -     -     -     -     -     -     -     -    
## bHLh -     -     -     -     -     -     -     -     -     -     -    
## bLhH -     -     -     -     -     -     -     -     -     -     -    
## bLHh -     -     -     -     -     -     -     -     -     -     -    
## hbHL -     -     -     -     -     -     -     -     -     -     -    
## HbhL -     -     -     -     -     -     -     -     -     -     -    
## hbLH -     -     -     -     -     -     -     -     -     -     -    
## HbLh -     -     -     -     -     -     -     -     -     -     -    
## hHbL -     -     -     -     -     -     -     -     -     -     -    
## HhbL -     -     -     -     -     -     -     -     -     -     -    
## hHLb -     -     -     -     -     -     -     -     -     -     -    
## HhLb 1.000 -     -     -     -     -     -     -     -     -     -    
## hLbH 0.909 1.000 -     -     -     -     -     -     -     -     -    
## HLbh 1.000 1.000 0.954 -     -     -     -     -     -     -     -    
## hLHb 1.000 1.000 1.000 1.000 -     -     -     -     -     -     -    
## HLhb 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -     -    
## LbhH 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -     -     -    
## LbHh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -     -    
## LhbH 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -     -    
## LHbh 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -     -    
## LhHb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 -    
## LHhb 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## 
## P value adjustment method: bonferroni
########environmental
dat_c2_stim<-droplevels(datac[datac$cnd=="c2",])

#create new factor
dat_c2_stim$stim_code<-as.factor(paste(substring(dat_c2_stim$stim1,1,1),substring(dat_c2_stim$stim2,1,1),substring(dat_c2_stim$stim3,1,1),substring(dat_c2_stim$stim4,1,1), sep=""))
levels(dat_c2_stim$stim_code)
##  [1] "bBdg" "Bbdg" "bBgd" "Bbgd" "bdBg" "Bdbg" "bdgB" "Bdgb" "bgBd" "Bgbd"
## [11] "bgdB" "Bgdb" "dbBg" "dBbg" "dbgB" "dBgb" "dgbB" "dgBb" "gbBd" "gBbd"
## [21] "gbdB" "gBdb" "gdbB" "gdBb"
#B=bark, b=bird, d=dog, g=glass breaking

#aggregate data
datc2_stim<-aggregate(dat_c2_stim$acc, list(dat_c2_stim$stim_code, dat_c2_stim$sbj), mean)
colnames(datc2_stim)<-c("stim_code", "sbj","acc")
ezANOVA(data=datc2_stim, dv=acc, wid=sbj, within=stim_code, type=3)
## $ANOVA
##      Effect DFn DFd        F            p p<.05       ges
## 2 stim_code  23 621 4.276912 2.366157e-10     * 0.1155254
## 
## $`Mauchly's Test for Sphericity`
##      Effect            W         p p<.05
## 2 stim_code 1.417044e-07 0.3434282      
## 
## $`Sphericity Corrections`
##      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 stim_code 0.4932389 4.742022e-06         * 0.8766074 2.596892e-09         *
#τhere is an effect, the command below reports partial eta squared (in line 2)
aovEffectSize(ezANOVA(data=datc1_stim, dv=acc, wid=sbj, within=stim_code, type=3, return_aov = TRUE, detailed = TRUE))
## $ANOVA
##        Effect DFn DFd       SSn      SSd         F            p p<.05       pes
## 1 (Intercept)   1  27 36.796968 23.97647 41.437216 6.743071e-07     * 0.6054778
## 2   stim_code  23 621  3.214193 28.69987  3.023819 3.473620e-06     * 0.1007140
## 
## $`Mauchly's Test for Sphericity`
##      Effect           W         p p<.05
## 2 stim_code 7.37825e-08 0.1637892      
## 
## $`Sphericity Corrections`
##      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2 stim_code 0.4790159 0.0007568019         * 0.8337478 1.899969e-05         *
## 
## $aov
## 
## Call:
## aov(formula = formula(aov_formula), data = data)
## 
## Grand Mean: 0.234003
## 
## Stratum 1: sbj
## 
## Terms:
##                 Residuals
## Sum of Squares   23.97647
## Deg. of Freedom        27
## 
## Residual standard error: 0.9423467
## 
## Stratum 2: sbj:stim_code
## 
## Terms:
##                 stim_code Residuals
## Sum of Squares   3.214193 28.699870
## Deg. of Freedom        23       621
## 
## Residual standard error: 0.2149781
## Estimated effects may be unbalanced
ezPlot(data=datc2_stim, dv=acc, wid=sbj, within=stim_code, type=3, x=stim_code)

#pairwise t-test, bonferroni corrected
pairwise.t.test(datc2_stim$acc, datc2_stim$stim_code, paired=T, p.adjust.method= "bonferroni")
## 
##  Pairwise comparisons using paired t tests 
## 
## data:  datc2_stim$acc and datc2_stim$stim_code 
## 
##      bBdg    Bbdg    bBgd    Bbgd    bdBg    Bdbg    bdgB    Bdgb    bgBd   
## Bbdg 1.00000 -       -       -       -       -       -       -       -      
## bBgd 1.00000 1.00000 -       -       -       -       -       -       -      
## Bbgd 1.00000 1.00000 1.00000 -       -       -       -       -       -      
## bdBg 1.00000 1.00000 1.00000 1.00000 -       -       -       -       -      
## Bdbg 1.00000 1.00000 1.00000 1.00000 1.00000 -       -       -       -      
## bdgB 0.25904 0.31144 1.00000 1.00000 1.00000 1.00000 -       -       -      
## Bdgb 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 -       -      
## bgBd 0.07804 0.10600 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 -      
## Bgbd 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## bgdB 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## Bgdb 0.18035 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## dbBg 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.14249 0.64426 0.46030
## dBbg 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.12574 1.00000 0.01582
## dbgB 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## dBgb 1.00000 1.00000 1.00000 0.94013 1.00000 1.00000 0.01145 0.27650 0.00016
## dgbB 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## dgBb 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.70274 1.00000 1.00000
## gbBd 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## gBbd 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.42907
## gbdB 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## gBdb 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## gdbB 0.02029 0.05129 0.17334 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## gdBb 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
##      Bgbd    bgdB    Bgdb    dbBg    dBbg    dbgB    dBgb    dgbB    dgBb   
## Bbdg -       -       -       -       -       -       -       -       -      
## bBgd -       -       -       -       -       -       -       -       -      
## Bbgd -       -       -       -       -       -       -       -       -      
## bdBg -       -       -       -       -       -       -       -       -      
## Bdbg -       -       -       -       -       -       -       -       -      
## bdgB -       -       -       -       -       -       -       -       -      
## Bdgb -       -       -       -       -       -       -       -       -      
## bgBd -       -       -       -       -       -       -       -       -      
## Bgbd -       -       -       -       -       -       -       -       -      
## bgdB 1.00000 -       -       -       -       -       -       -       -      
## Bgdb 1.00000 1.00000 -       -       -       -       -       -       -      
## dbBg 0.60917 1.00000 0.37039 -       -       -       -       -       -      
## dBbg 1.00000 1.00000 0.17891 1.00000 -       -       -       -       -      
## dbgB 1.00000 1.00000 1.00000 1.00000 1.00000 -       -       -       -      
## dBgb 0.29690 0.18183 0.15576 1.00000 1.00000 0.72274 -       -       -      
## dgbB 1.00000 1.00000 1.00000 1.00000 0.82889 1.00000 0.14249 -       -      
## dgBb 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 -      
## gbBd 1.00000 1.00000 1.00000 1.00000 0.95550 1.00000 0.17334 1.00000 1.00000
## gBbd 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## gbdB 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## gBdb 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000
## gdbB 1.00000 1.00000 1.00000 0.00613 0.00841 1.00000 0.00117 1.00000 0.18568
## gdBb 1.00000 1.00000 1.00000 1.00000 1.00000 1.00000 0.60404 1.00000 1.00000
##      gbBd    gBbd    gbdB    gBdb    gdbB   
## Bbdg -       -       -       -       -      
## bBgd -       -       -       -       -      
## Bbgd -       -       -       -       -      
## bdBg -       -       -       -       -      
## Bdbg -       -       -       -       -      
## bdgB -       -       -       -       -      
## Bdgb -       -       -       -       -      
## bgBd -       -       -       -       -      
## Bgbd -       -       -       -       -      
## bgdB -       -       -       -       -      
## Bgdb -       -       -       -       -      
## dbBg -       -       -       -       -      
## dBbg -       -       -       -       -      
## dbgB -       -       -       -       -      
## dBgb -       -       -       -       -      
## dgbB -       -       -       -       -      
## dgBb -       -       -       -       -      
## gbBd -       -       -       -       -      
## gBbd 1.00000 -       -       -       -      
## gbdB 1.00000 1.00000 -       -       -      
## gBdb 1.00000 1.00000 1.00000 -       -      
## gdbB 1.00000 1.00000 1.00000 0.31144 -      
## gdBb 1.00000 1.00000 1.00000 1.00000 1.00000
## 
## P value adjustment method: bonferroni
Musical Training
########abstract sounds
#calculate average accuracy per sbj and musical training
datc1_m<-aggregate(dat_c1_stim$acc, list(dat_c1_stim$music, dat_c1_stim$sbj), mean)
colnames(datc1_m)<-c("music", "sbj","acc")
#some musical training
mean(datc1_m[datc1_m$music=="yes",]$acc)
## [1] 0.3238636
sd(datc1_m[datc1_m$music=="yes",]$acc)
## [1] 0.2329449
#no musical training
mean(datc1_m[datc1_m$music=="no",]$acc)
## [1] 0.1758578
sd(datc1_m[datb1_m$music=="no",]$acc)
## [1] 0.1392137
#check if there's a difference between groups
var.test(acc~music, data=datc1_m)
## 
##  F test to compare two variances
## 
## data:  acc by music
## F = 0.35716, num df = 16, denom df = 10, p-value = 0.06465
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.1021533 1.0665251
## sample estimates:
## ratio of variances 
##          0.3571557
# equal variances
t.test(acc~music,paired=F, var.equal=T,data=datc1_m)
## 
##  Two Sample t-test
## 
## data:  acc by music
## t = -2.112, df = 26, p-value = 0.04445
## alternative hypothesis: true difference in means between group no and group yes is not equal to 0
## 95 percent confidence interval:
##  -0.292050942 -0.003960644
## sample estimates:
##  mean in group no mean in group yes 
##         0.1758578         0.3238636
#there's a difference 
datc1_m %>% cohens_d(acc~music, paired=F)
## # A tibble: 1 x 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 acc   no     yes     -0.771    17    11 moderate
#functions to draw hatched boxplots
#modified from https://stackoverflow.com/questions/17423115/colorfill-boxplot-in-r-cran-with-lines-dots-or-similar
boxpattern_2 <- function(y, xcenter, boxwidth, angle=NULL, angle.density=10, ...) {
  # draw an individual box
  bstats <- boxplot.stats(y)
  bxmin <- bstats$stats[1]
  bxq2 <- bstats$stats[2]
  bxmedian <- bstats$stats[3]
  bxq4 <- bstats$stats[4]
  bxmax <- bstats$stats[5]
  bleft <- xcenter-(boxwidth/2)
  bright <- xcenter+(boxwidth/2)
  bleft_seg <- xcenter-(boxwidth/4)
  bright_seg <- xcenter+(boxwidth/4)
  # boxplot
  polygon(c(bleft,bright,bright,bleft,bleft),
          c(bxq2,bxq2,bxq4,bxq4,bxq2), angle=angle[1], col="grey75")
  polygon(c(bleft,bright,bright,bleft,bleft),
          c(bxq2,bxq2,bxq4,bxq4,bxq2), angle=angle[2], density=angle.density)
  # lines
  segments(bleft,bxmedian,bright,bxmedian,lwd=3, lty=5) # median
  segments(bleft_seg,bxmin,bright_seg,bxmin,lwd=1, lty=5) # min
  segments(xcenter,bxmin,xcenter,bxq2,lwd=1, lty=5)
  segments(bleft_seg,bxmax,bright_seg,bxmax,lwd=1, lty=5) # max
  segments(xcenter,bxq4,xcenter,bxmax,lwd=1, lty=5)
}

drawboxplots_2 <- function(y, x, boxwidth=1, angle=NULL, ...){
  # figure out all the boxes and start the plot
  groups <- split(y,as.factor(x))
  len <- length(groups)
  bxylim <- c((min(y)-0.04*abs(min(y))),(max(y)+0.04*max(y)))
  xcenters <- seq(1,max(2,(len*(1.4))),length.out=len)
  if(is.null(angle)){
    angle <- seq(-90,35,length.out=len)
    angle <- lapply(angle,function(x) c(x,x))
  }
  else if(!length(angle)==len)
    stop("angle must be a vector or list of two-element vectors")
  else if(!is.list(angle))
    angle <- lapply(angle,function(x) c(x,x))
  # draw plot area
  plot(10, xlim=c(.97*(min(xcenters)-1), 1.04*(max(xcenters)+1)),
       #ylim=bxylim, 
       ylim=c(0.0, 1.0),
       xlab="",
       ylab="Proportion Correct", 
       las=1, 
       frame.plot=F, 
       xaxt="n", 
       main="Ordering Task, Abstract Sounds\nAccuracy by Musical Training")
  
  axis(1, at=xcenters, tick=F, labels=c("None\nN=11", "Some\nN=17"))
  # draw boxplots
  plots <- mapply(boxpattern_2, y=groups, xcenter=xcenters,
                  boxwidth=boxwidth, angle=angle, ...)
}

drawboxplots_2(datc1_m$acc, datc1_m$music)

Accuracy per Sound Position
#******************
#create data frame
#******************

sbj<-rep("S",28*192*4)
acc<-rep(1,28*192*4)
trial<-rep(1,28*192*4)
sound<-rep("a",28*192*4)
stim<-rep("stim",28*192*4)
cnd<-rep("C",28*192*4)


warren<-data.frame(sbj, cnd, trial, sound, stim, acc)
#Convert all chr columns to factor
warren <- mutate_if(warren, is.character, as.factor)

levels(warren$sbj)<-levels(datac$sbj)
levels(warren$cnd)<-levels(datac$cnd)
levels(warren$stim)<-levels(datac$stim1)
levels(warren$sound)<-c("1st", "2nd", "3rd","4th")

str(warren)
## 'data.frame':    21504 obs. of  6 variables:
##  $ sbj  : Factor w/ 28 levels "sbj1","sbj10",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ cnd  : Factor w/ 2 levels "c1","c2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ trial: num  1 1 1 1 1 1 1 1 1 1 ...
##  $ sound: Factor w/ 4 levels "1st","2nd","3rd",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ stim : Factor w/ 8 levels "Bark","bird",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ acc  : num  1 1 1 1 1 1 1 1 1 1 ...
for (i in 1:length(datac$sbj))
  {

    warren[1+4*(i-1),]$sound<-"1st"
    warren[2+4*(i-1),]$sound<-"2nd"
    warren[3+4*(i-1),]$sound<-"3rd"
    warren[4+4*(i-1),]$sound<-"4th"

    warren[1+4*(i-1),]$trial<-datac[i,]$trial+1
    warren[2+4*(i-1),]$trial<-datac[i,]$trial+1
    warren[3+4*(i-1),]$trial<-datac[i,]$trial+1
    warren[4+4*(i-1),]$trial<-datac[i,]$trial+1

    warren[1+4*(i-1),]$acc<-datac[i,]$acc1
    warren[2+4*(i-1),]$acc<-datac[i,]$acc2
    warren[3+4*(i-1),]$acc<-datac[i,]$acc3
    warren[4+4*(i-1),]$acc<-datac[i,]$acc4

    warren[1+4*(i-1),]$sbj<-datac[i,]$sbj
    warren[2+4*(i-1),]$sbj<-datac[i,]$sbj
    warren[3+4*(i-1),]$sbj<-datac[i,]$sbj
    warren[4+4*(i-1),]$sbj<-datac[i,]$sbj

    warren[1+4*(i-1),]$cnd<-datac[i,]$cnd
    warren[2+4*(i-1),]$cnd<-datac[i,]$cnd
    warren[3+4*(i-1),]$cnd<-datac[i,]$cnd
    warren[4+4*(i-1),]$cnd<-datac[i,]$cnd

    warren[1+4*(i-1),]$stim<-datac[i,]$stim1
    warren[2+4*(i-1),]$stim<-datac[i,]$stim2
    warren[3+4*(i-1),]$stim<-datac[i,]$stim3
    warren[4+4*(i-1),]$stim<-datac[i,]$stim4

    } 

d<-warren


#Analysis of accuracy per sound position
d_sbj_order<-aggregate(d$acc,list(d$sbj, d$cnd, d$sound), mean)
colnames(d_sbj_order)<-c("sbj", "cnd", "sound", "acc")
#draw boxplot
col=c("grey75","grey25")
boxplot(acc~cnd*sound, at=c(1,2,4,5,7,8,10,11), data=d_sbj_order, ylim=c(0.0, 1.0),main="Ordering Tasks", xlab="Sound Position", ylab="Proportion Correct", frame.plot=F, las=2, col=col, range=0, xaxt="n", cex.axis=1, lwd=1, lty=5)
axis(side=1, at=c(1.5,4.5,7.5,10.5),labels=c("1st","2nd","3rd","4th"))
legend(x=7.5, y=1.05, legend=c("Abstract Sounds", "Environmental Sounds"), bty= "n", fill=col)

#split data per sound category (cnd)
########abstract
d_sbj_order_c1<-droplevels(d_sbj_order[d_sbj_order$cnd=="c1",])
str(d_sbj_order_c1)
## 'data.frame':    112 obs. of  4 variables:
##  $ sbj  : Factor w/ 28 levels "sbj1","sbj10",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ cnd  : Factor w/ 1 level "c1": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sound: Factor w/ 4 levels "1st","2nd","3rd",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ acc  : num  0.76 0.375 0.594 0.781 0.365 ...
#check if there's an effect of sound position
ezANOVA(data=d_sbj_order_c1, dv=acc, wid=sbj, within=sound, type=3)
## $ANOVA
##   Effect DFn DFd        F            p p<.05       ges
## 2  sound   3  81 59.14061 2.372417e-20     * 0.2047266
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W           p p<.05
## 2  sound 0.5554823 0.009906582     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05      HFe        p[HF] p[HF]<.05
## 2  sound 0.7967395 1.105746e-16         * 0.879408 3.553499e-18         *
#there is an effect
#report corrects dfs
0.7967395*3
## [1] 2.390219
0.7967395*81
## [1] 64.5359
#report partial eta squared
aovEffectSize(ezANOVA(data=d_sbj_order_c1, dv=acc, wid=sbj, within=sound, type=3, return_aov = TRUE, detailed = TRUE))
## $ANOVA
##        Effect DFn DFd        SSn       SSd         F            p p<.05
## 1 (Intercept)   1  27 25.3175223 2.7598431 247.68549 4.001432e-15     *
## 2       sound   3  81  0.8050828 0.3675518  59.14061 2.372417e-20     *
##         pes
## 1 0.9017058
## 2 0.6865590
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W           p p<.05
## 2  sound 0.5554823 0.009906582     *
## 
## $`Sphericity Corrections`
##   Effect       GGe        p[GG] p[GG]<.05      HFe        p[HF] p[HF]<.05
## 2  sound 0.7967395 1.105746e-16         * 0.879408 3.553499e-18         *
## 
## $aov
## 
## Call:
## aov(formula = formula(aov_formula), data = data)
## 
## Grand Mean: 0.4754464
## 
## Stratum 1: sbj
## 
## Terms:
##                 Residuals
## Sum of Squares   2.759843
## Deg. of Freedom        27
## 
## Residual standard error: 0.319713
## 
## Stratum 2: sbj:sound
## 
## Terms:
##                     sound Residuals
## Sum of Squares  0.8050828 0.3675518
## Deg. of Freedom         3        81
## 
## Residual standard error: 0.06736228
## Estimated effects may be unbalanced
#pairwise comparisons
pairwise.t.test(d_sbj_order_c1$acc, d_sbj_order_c1$sound, paired=T, p.adjust.method= "bonferroni")
## 
##  Pairwise comparisons using paired t tests 
## 
## data:  d_sbj_order_c1$acc and d_sbj_order_c1$sound 
## 
##     1st     2nd    3rd    
## 2nd 2.0e-10 -      -      
## 3rd 9.6e-11 0.7841 -      
## 4th 2.3e-05 0.0058 3.6e-06
## 
## P value adjustment method: bonferroni
#1st sound > 2nd, 3rd, 4th. 
#4th > 2nd, 3rd
########environmental 
d_sbj_order_c2<-droplevels(d_sbj_order[d_sbj_order$cnd=="c2",])
str(d_sbj_order_c2)
## 'data.frame':    112 obs. of  4 variables:
##  $ sbj  : Factor w/ 28 levels "sbj1","sbj10",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ cnd  : Factor w/ 1 level "c2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sound: Factor w/ 4 levels "1st","2nd","3rd",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ acc  : num  0.729 0.781 0.75 0.927 0.708 ...
#check if there's an effect of sound position
ezANOVA(data=d_sbj_order_c2, dv=acc, wid=sbj, within=sound, type=3)
## $ANOVA
##   Effect DFn DFd        F           p p<.05       ges
## 2  sound   3  81 181.7379 7.19857e-36     * 0.6917832
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W            p p<.05
## 2  sound 0.3019402 1.046479e-05     *
## 
## $`Sphericity Corrections`
##   Effect       GGe       p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2  sound 0.6463761 4.98774e-24         * 0.6955809 1.117794e-25         *
#there is an effect
#report corrects dfs
0.6463761*3
## [1] 1.939128
0.6463761*81
## [1] 52.35646
#report partial eta squared
aovEffectSize(ezANOVA(data=d_sbj_order_c2, dv=acc, wid=sbj, within=sound, type=3, return_aov = TRUE, detailed = TRUE))
## $ANOVA
##        Effect DFn DFd       SSn       SSd         F            p p<.05
## 1 (Intercept)   1  27 33.155169 0.8269818 1082.4779 2.520714e-23     *
## 2       sound   3  81  2.784695 0.4137099  181.7379 7.198570e-36     *
##         pes
## 1 0.9756642
## 2 0.8706512
## 
## $`Mauchly's Test for Sphericity`
##   Effect         W            p p<.05
## 2  sound 0.3019402 1.046479e-05     *
## 
## $`Sphericity Corrections`
##   Effect       GGe       p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 2  sound 0.6463761 4.98774e-24         * 0.6955809 1.117794e-25         *
## 
## $aov
## 
## Call:
## aov(formula = formula(aov_formula), data = data)
## 
## Grand Mean: 0.5440848
## 
## Stratum 1: sbj
## 
## Terms:
##                 Residuals
## Sum of Squares  0.8269818
## Deg. of Freedom        27
## 
## Residual standard error: 0.1750113
## 
## Stratum 2: sbj:sound
## 
## Terms:
##                     sound Residuals
## Sum of Squares  2.7846951 0.4137099
## Deg. of Freedom         3        81
## 
## Residual standard error: 0.07146698
## Estimated effects may be unbalanced
#pairwise comparisons
pairwise.t.test(d_sbj_order_c2$acc, d_sbj_order_c2$sound, paired=T, p.adjust.method= "bonferroni")
## 
##  Pairwise comparisons using paired t tests 
## 
## data:  d_sbj_order_c2$acc and d_sbj_order_c2$sound 
## 
##     1st     2nd    3rd    
## 2nd < 2e-16 -      -      
## 3rd < 2e-16 0.0045 -      
## 4th 4.3e-10 0.0010 1.5e-09
## 
## P value adjustment method: bonferroni
#1st sound > 2nd, 3rd, 4th.
#2nd > 3rd
#4th > 2nd, 3rd
#################
#two-way ANOVA
#################
#check if there is an interaction of cnd by sound
ezANOVA(data=d_sbj_order, dv=acc, wid=sbj, within=.(cnd,sound), type=3)
## $ANOVA
##      Effect DFn DFd         F            p p<.05       ges
## 2       cnd   1  27  13.51414 1.035588e-03     * 0.0569589
## 3     sound   3  81 174.23835 3.159680e-35     * 0.4290272
## 4 cnd:sound   3  81  30.46142 2.769595e-13     * 0.0657894
## 
## $`Mauchly's Test for Sphericity`
##      Effect         W            p p<.05
## 3     sound 0.3650997 0.0000941661     *
## 4 cnd:sound 0.3752103 0.0001286879     *
## 
## $`Sphericity Corrections`
##      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     sound 0.7004986 2.158972e-25         * 0.7609812 2.220561e-27         *
## 4 cnd:sound 0.7152290 4.007455e-10         * 0.7789291 7.842104e-11         *
#there is an interaction
#report corrected dfs
#for simple effect of sound
0.7004986*3
## [1] 2.101496
0.7004986*81
## [1] 56.74039
#for interaction
0.7152290*3
## [1] 2.145687
0.7152290*81
## [1] 57.93355
#report eta-squared
aovEffectSize(ezANOVA(data=d_sbj_order, dv=acc, wid=sbj, within=.(cnd,sound), type=3, return_aov = TRUE, detailed = TRUE))
## $ANOVA
##        Effect DFn DFd        SSn       SSd         F            p p<.05
## 1 (Intercept)   1  27 58.2088623 3.0597195 513.65469 4.183786e-19     *
## 2         cnd   1  27  0.2638288 0.5271054  13.51414 1.035588e-03     *
## 3       sound   3  81  3.2821665 0.5086050 174.23835 3.159680e-35     *
## 4   cnd:sound   3  81  0.3076114 0.2726566  30.46142 2.769595e-13     *
##         pes
## 1 0.9500605
## 2 0.3335661
## 3 0.8658307
## 4 0.5301195
## 
## $`Mauchly's Test for Sphericity`
##      Effect         W            p p<.05
## 3     sound 0.3650997 0.0000941661     *
## 4 cnd:sound 0.3752103 0.0001286879     *
## 
## $`Sphericity Corrections`
##      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     sound 0.7004986 2.158972e-25         * 0.7609812 2.220561e-27         *
## 4 cnd:sound 0.7152290 4.007455e-10         * 0.7789291 7.842104e-11         *
## 
## $aov
## 
## Call:
## aov(formula = formula(aov_formula), data = data)
## 
## Grand Mean: 0.5097656
## 
## Stratum 1: sbj
## 
## Terms:
##                 Residuals
## Sum of Squares    3.05972
## Deg. of Freedom        27
## 
## Residual standard error: 0.3366347
## 
## Stratum 2: sbj:cnd
## 
## Terms:
##                       cnd Residuals
## Sum of Squares  0.2638288 0.5271054
## Deg. of Freedom         1        27
## 
## Residual standard error: 0.1397227
## 3 out of 4 effects not estimable
## Estimated effects are balanced
## 
## Stratum 3: sbj:sound
## 
## Terms:
##                    sound Residuals
## Sum of Squares  3.282166  0.508605
## Deg. of Freedom        3        81
## 
## Residual standard error: 0.07924061
## 3 out of 6 effects not estimable
## Estimated effects may be unbalanced
## 
## Stratum 4: sbj:cnd:sound
## 
## Terms:
##                 cnd:sound Residuals
## Sum of Squares  0.3076114 0.2726566
## Deg. of Freedom         3        81
## 
## Residual standard error: 0.05801837
## Estimated effects may be unbalanced
#there is an interaction
#report differences of previously done pairwise comparisons

#abstract
#mean diff between 2 and 3
mean(d_sbj_order_c1[d_sbj_order_c1$sound=="2nd",]$acc - d_sbj_order_c1[d_sbj_order_c1$sound=="3rd",]$acc)
## [1] 0.01971726
sd(d_sbj_order_c1[d_sbj_order_c1$sound=="2nd",]$acc - d_sbj_order_c1[d_sbj_order_c1$sound=="3rd",]$acc)
## [1] 0.06693014
# p = .78

#environmental
#diff between 2 and 3
mean(d_sbj_order_c2[d_sbj_order_c2$sound=="2nd",]$acc - d_sbj_order_c2[d_sbj_order_c2$sound=="3rd",]$acc)
## [1] 0.04427083
sd(d_sbj_order_c2[d_sbj_order_c2$sound=="2nd",]$acc - d_sbj_order_c2[d_sbj_order_c2$sound=="3rd",]$acc)
## [1] 0.06169915
# p = .005
####################################
#analysis of consecutive differences
####################################

ezANOVA(data=d_sbj_order, dv=acc, wid=sbj, within=.(cnd,sound), type=3)
## $ANOVA
##      Effect DFn DFd         F            p p<.05       ges
## 2       cnd   1  27  13.51414 1.035588e-03     * 0.0569589
## 3     sound   3  81 174.23835 3.159680e-35     * 0.4290272
## 4 cnd:sound   3  81  30.46142 2.769595e-13     * 0.0657894
## 
## $`Mauchly's Test for Sphericity`
##      Effect         W            p p<.05
## 3     sound 0.3650997 0.0000941661     *
## 4 cnd:sound 0.3752103 0.0001286879     *
## 
## $`Sphericity Corrections`
##      Effect       GGe        p[GG] p[GG]<.05       HFe        p[HF] p[HF]<.05
## 3     sound 0.7004986 2.158972e-25         * 0.7609812 2.220561e-27         *
## 4 cnd:sound 0.7152290 4.007455e-10         * 0.7789291 7.842104e-11         *
#there is an interaction:
#Pairwise comparisons
a<-aov(acc~sound*cnd + Error(sbj/(sound + cnd)), data=d_sbj_order)
m2<-emmeans(a, ~sound|cnd)
## Note: re-fitting model with sum-to-zero contrasts
con2<-contrast(m2, "consec")

#all comparisons, this is what we need.
pairs(con2, by=NULL)
##  contrast                        estimate     SE  df t.ratio p.value
##  (2nd - 1st c1) - (3rd - 2nd c1)  -0.1767 0.0321 148  -5.497  <.0001
##  (2nd - 1st c1) - (4th - 3rd c1)  -0.2909 0.0262 148 -11.084  <.0001
##  (2nd - 1st c1) - (2nd - 1st c2)   0.1656 0.0219  81   7.549  <.0001
##  (2nd - 1st c1) - (3rd - 2nd c2)  -0.1522 0.0282 109  -5.403  <.0001
##  (2nd - 1st c1) - (4th - 3rd c2)  -0.3344 0.0262 148 -12.742  <.0001
##  (3rd - 2nd c1) - (4th - 3rd c1)  -0.1142 0.0321 148  -3.553  0.0067
##  (3rd - 2nd c1) - (2nd - 1st c2)   0.3423 0.0282 109  12.154  <.0001
##  (3rd - 2nd c1) - (3rd - 2nd c2)   0.0246 0.0219  81   1.120  0.8719
##  (3rd - 2nd c1) - (4th - 3rd c2)  -0.1577 0.0282 109  -5.602  <.0001
##  (4th - 3rd c1) - (2nd - 1st c2)   0.4565 0.0262 148  17.391  <.0001
##  (4th - 3rd c1) - (3rd - 2nd c2)   0.1388 0.0282 109   4.928  <.0001
##  (4th - 3rd c1) - (4th - 3rd c2)  -0.0435 0.0219  81  -1.985  0.3599
##  (2nd - 1st c2) - (3rd - 2nd c2)  -0.3177 0.0321 148  -9.883  <.0001
##  (2nd - 1st c2) - (4th - 3rd c2)  -0.5000 0.0262 148 -19.049  <.0001
##  (3rd - 2nd c2) - (4th - 3rd c2)  -0.1823 0.0321 148  -5.671  <.0001
## 
## P value adjustment: tukey method for comparing a family of 6 estimates
#Additional comparisons of sound category, separated by sound position 
#1st sound
d_sbj_order_1st<-droplevels(d_sbj_order[d_sbj_order$sound=="1st",])
t.test(d_sbj_order_1st$acc~ d_sbj_order_1st$cnd, paired=T)
## 
##  Paired t-test
## 
## data:  d_sbj_order_1st$acc by d_sbj_order_1st$cnd
## t = -6.4213, df = 27, p-value = 7.027e-07
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.2562493 -0.1321436
## sample estimates:
## mean difference 
##      -0.1941964
#difference

mean(d_sbj_order_1st[d_sbj_order_1st$cnd=="c1",]$acc) 
## [1] 0.609003
sd(d_sbj_order_1st[d_sbj_order_1st$cnd=="c1",]$acc) 
## [1] 0.1947023
mean(d_sbj_order_1st[d_sbj_order_1st$cnd=="c2",]$acc) 
## [1] 0.8031994
sd(d_sbj_order_1st[d_sbj_order_1st$cnd=="c2",]$acc) 
## [1] 0.1238482
#2nd sound
d_sbj_order_2nd<-droplevels(d_sbj_order[d_sbj_order$sound=="2nd",])
t.test(d_sbj_order_2nd$acc~ d_sbj_order_2nd$cnd, paired=T)
## 
##  Paired t-test
## 
## data:  d_sbj_order_2nd$acc by d_sbj_order_2nd$cnd
## t = -1.2306, df = 27, p-value = 0.2291
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.07640780  0.01911613
## sample estimates:
## mean difference 
##     -0.02864583
#3rd sound
d_sbj_order_3rd<-droplevels(d_sbj_order[d_sbj_order$sound=="3rd",])
t.test(d_sbj_order_3rd$acc~ d_sbj_order_3rd$cnd, paired=T)
## 
##  Paired t-test
## 
## data:  d_sbj_order_3rd$acc by d_sbj_order_3rd$cnd
## t = -0.2432, df = 27, p-value = 0.8097
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.03861757  0.03043304
## sample estimates:
## mean difference 
##    -0.004092262
#4th sound
d_sbj_order_4th<-droplevels(d_sbj_order[d_sbj_order$sound=="4th",])
t.test(d_sbj_order_4th$acc~ d_sbj_order_4th$cnd, paired=T)
## 
##  Paired t-test
## 
## data:  d_sbj_order_4th$acc by d_sbj_order_4th$cnd
## t = -2.4552, df = 27, p-value = 0.02081
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  -0.087414278 -0.007823817
## sample estimates:
## mean difference 
##     -0.04761905
#doesn't survive Bonferroni corrections for four tests

Session Information

sessionInfo()
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=Greek_Greece.1253  LC_CTYPE=Greek_Greece.1253   
## [3] LC_MONETARY=Greek_Greece.1253 LC_NUMERIC=C                 
## [5] LC_TIME=Greek_Greece.1253    
## system code page: 65001
## 
## time zone: Europe/Athens
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] emmeans_1.8.8     psychReport_3.0.2 rstatix_0.7.2     ez_4.4-0         
## [5] sciplot_1.2-0     dplyr_1.1.3      
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.4       xfun_0.40          bslib_0.5.1        ggplot2_3.4.3     
##  [5] lattice_0.21-8     vctrs_0.6.3        tools_4.3.1        generics_0.1.3    
##  [9] sandwich_3.0-2     tibble_3.2.1       fansi_1.0.4        pkgconfig_2.0.3   
## [13] Matrix_1.6-1       lifecycle_1.0.3    farver_2.1.1       compiler_4.3.1    
## [17] stringr_1.5.0      munsell_0.5.0      codetools_0.2-19   carData_3.0-5     
## [21] htmltools_0.5.6    sass_0.4.7         yaml_2.3.7         pillar_1.9.0      
## [25] car_3.1-2          nloptr_2.0.3       jquerylib_0.1.4    tidyr_1.3.0       
## [29] MASS_7.3-60        cachem_1.0.8       boot_1.3-28.1      abind_1.4-5       
## [33] multcomp_1.4-25    nlme_3.1-162       tidyselect_1.2.0   digest_0.6.33     
## [37] mvtnorm_1.2-3      stringi_1.7.12     reshape2_1.4.4     purrr_1.0.2       
## [41] labeling_0.4.3     splines_4.3.1      fastmap_1.1.1      grid_4.3.1        
## [45] colorspace_2.1-0   cli_3.6.1          magrittr_2.0.3     survival_3.5-5    
## [49] utf8_1.2.3         broom_1.0.5        TH.data_1.1-2      withr_2.5.0       
## [53] scales_1.2.1       backports_1.4.1    estimability_1.4.1 rmarkdown_2.24    
## [57] lme4_1.1-34        zoo_1.8-12         evaluate_0.21      knitr_1.44        
## [61] mgcv_1.8-42        rlang_1.1.1        Rcpp_1.0.11        xtable_1.8-4      
## [65] glue_1.6.2         rstudioapi_0.15.0  minqa_1.2.5        jsonlite_1.8.7    
## [69] R6_2.5.1           plyr_1.8.8