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
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 "alekws23","anadim22",..: 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 "alekws23","anadim22",..: 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.802 0.542 0.573 0.667 ...
#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
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 "alekws23","anadim22",..: 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.823 0.812 0.844 0.844 ...
#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
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
#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
####################################
#analysis of consecutive differences
####################################
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
