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 ...
################
# 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
#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" ))
#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
########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
########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 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
#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
########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
########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)
#******************
#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
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