# Working on making fancy graphs with R / fixed versus random babble

I have been working on learning R for several months now, and continue to get better at it and enjoy it more all the time. I am currently working on a spoken word recognition project at work. The task we are using is quite simple. Participants listen to words that have been mixed with multi-talker babble (kind of like background conversation at a cocktail party), and type in what they hear. We are analyzing the errors they make to try to discover what sorts of words are activated in the brain, and how people organize the words they know in their brains.

We started running subjects in the summer of 2007, and after running 48 subjects, we made a few methodological changes. One change is that we switched from using the same segment of babble as background noise to using a random segment of babble. We also re-leveled the sound files after adding the babble so that all the words were presented at about the same volume. I recently analyzed the data from these two groups. It turns out the the results from the first 48 subjects performed significantly better than subjects 49-96. However, I realized that these are different subjects, so it could be due to the different subjects, and not the change in methods. Therefore I split each group into 2 subgroups. Now I had 4 groups — 1-24, 24-48, 49-72, and 73-96. The subgroups did not differ statistically from one another, but the main groups did. This seemed to indicate that listeners were becoming habituated to the noise, and could therefore tune it out a bit better, and thus hear the words better.

I originally displayed this in a table, but it was pretty cumbersome, so I decided to make a bar chart instead. I have seen some people produce bar charts such as these where they indicate which groups are statistically different using brackets to group bars, and an asterisk to mark them as statistically different. I also wanted to add error bars in. Adam Krawitz (another post-doc at IU) gave me some code to get started on the error bars, and I found some code from a post to the R-help list about making brackets. First here is the figure:

And here is the code:

# make a nice bar plot

# first we find the mean of each group

bardata <- c(mean(items1to24$pc,na.rm=T),
mean(items25to48$pc,na.rm=T),
mean(items49to72$pc,na.rm=T),
mean(items73to96$pc,na.rm=T)
)
#Then we do t-tests to calculate the 95% confidence intervals for our error bars
items1to24T=t.test(items1to24$pc)
items25to48T=t.test(items25to48$pc)
items49to72T=t.test(items49to72$pc)
items73to96T=t.test(items73to96$pc)
lowerConf=c(items1to24T$conf.int[1],
items25to48T$conf.int[1],
items49to72T$conf.int[1],
items73to96T$conf.int[1]
)
upperConf=c(items1to24T$conf.int[2],
items25to48T$conf.int[2],
items49to72T$conf.int[2],
items73to96T$conf.int[2]
)
# we set up some nice colors
barcolors=c('#FF3333','#FF3333','#4444FF','#4444FF')
bar <- barplot(bardata, ylim = c(0,1), col=barcolors, names.arg=c("1-24","25-48","49-72", "73-96"),ylab='proportion correct', xlab='subject group',xlim=c(0,1), width=.24,tck=.04)
#we use the arrows function to do the error bars
arrows(bar, upperConf, bar, lowerConf,
length = 0.05, # width of the arrowhead
angle = 90,
code = 3
)
# function to draw curly braces in red
# x1...y2 are the ends of the brace
# for upside down braces, x1 > x2 and y1 > y2

# taken from R help

# http://finzi.psych.upenn.edu/R/Rhelp02a/archive/112010.html

library(grid) #requires the grid library

Brack <- function(x1,y1,x2,y2,h)
{
x2 <- x2-x1; y2 <- y2-y1
v1 <- viewport(x=x1,y=y1,width=sqrt(x2^2+y2^2),
height=h,angle=180*atan2(y2,x2)/pi,
just=c("left","bottom"),gp=gpar(col="black"))
pushViewport(v1)
grid.curve(x2=0,y2=0,x1=.125,y1=.5,curvature=.5)
grid.move.to(.125,.5)
grid.line.to(.375,.5)
grid.curve(x1=.375,y1=.5,x2=.5,y2=1,curvature=.5)
grid.curve(x2=1,y2=0,x1=.875,y1=.5,curvature=-.5)
grid.move.to(.875,.5)
grid.line.to(.625,.5)
grid.curve(x2=.625,y2=.5,x1=.5,y1=1,curvature=.5)
popViewport()
}
# Now we use the brack function
# brackY1 and brack Y2 determine the y-value for the placement of the brackets.
# I had to fudge this a bit to get it to look right
brackY1=1.05*max(upperConf[1:2])
brackY2=1.21*max(upperConf[3:4])
Brack(bar[1]+.07,brackY1,bar[2]-.01,brackY1,.07)
Brack(bar[3]-.09,brackY2,bar[4]-.17,brackY2,.11)
# And now we also make some straight line groupings
arrows(mean(bar[1:2]),1.4*upperConf[1],mean(bar[3:4]), 1.4*upperConf[1],length=0)
arrows(mean(bar[1:2]),1.4*upperConf[1],mean(bar[1:2]), 1.35*upperConf[1],length=0)
arrows(mean(bar[3:4]),1.4*upperConf[1],mean(bar[3:4]), 1.35*upperConf[1],length=0)
arrows(mean(bar[2:3]),1.4*upperConf[1],mean(bar[2:3]), 1.45*upperConf[1],length=0)
# we label the groupings
text(mean(bar[1:2]),1.25*upperConf[1],"fixed babble",cex=.9,col=barcolors[1])
text(mean(bar[3:4]),1.25*upperConf[1],"random babble",cex=.9,col=barcolors[3])
# we add the stats in
text(mean(bar[2:3]),1.55*upperConf[1],paste("t =", formatC(itemBetweenT$statistic,digits=3), ", p", prettyPval(itemBetweenT$p.value)),cex=.9)
# and we add the values on the bars as well
text(bar[1],.7*bardata[1],formatC(bardata[1],digits=3),col='white')
text(bar[2],.7*bardata[2],formatC(bardata[2],digits=3),col='white')
text(bar[3],.7*bardata[3],formatC(bardata[3],digits=3),col='white')
text(bar[4],.7*bardata[4],formatC(bardata[4],digits=3),col='white')
# and then I print it to a pdf file
dev.print(pdf,"fig/48v96.pdf",height=2.8,width=4,pointsize=11)

While these results were suggestive of a difference, it is also possible that the difference was due to the re-leveling, and not the fixed vs. random babble. If the listeners were really habituating to the fixed babble, this should show up as an improvement over the course of the experiment. Therefore I also compared the learning rates of these two groups. It is normal for participants to get better as they get more familiar with a task, so we will always expect some learning. However if the participants are becoming habituated to a particular aspect of the stimuli, (in this case the fact that the same segment of babble is being used over and over again), then we would expect more learning in this group. This figure displays the percent correct over a moving 30 trial window. That is, the first point represents trials 1-30, the second point 2-31 and so on. It looks like the fixed-babble group is learning more. To test this statistically, I subtracted the fixed babble values from the random babble values, and performed a correlation on these values against the trial window. If the learning rates are the same, then there should be no correlation. If the fixed babble group is learning more rapidly however, we should see the difference between the two groups increase over the course of the experiment, which is what we found. Here is the nice figure:

And here is the code:

#compute the percent correct for each window

numPoints=min(length(unique(Trials1$Trial)),length(unique(Trials1$Trial)))-trialWindow

pc1=c(1:numPoints)

pc2=c(1:numPoints)

xvalues=c(1:numPoints)

for (i in 0:numPoints) {

upperTrials1=subset(Trials1,Trials1$Trial>i)

theseTrials1=subset(upperTrials1,upperTrials1$Trial<=i+trialWindow)
corr1 = (nrow(subset(theseTrials1, theseTrials1$correct=="CORRECT")))
wrong1 = (nrow(subset(theseTrials1,theseTrials1$correct=="WRONG")))
missing1 = (nrow(subset(theseTrials1,theseTrials1$correct=="MISSING")))
nonword1 = (nrow(subset(theseTrials1,theseTrials1$correct=="NONWORD")))
pc1[i] = corr1/(corr1+wrong1+nonword1+missing1)
upperTrials2=subset(Trials2,Trials2$Trial>i)

theseTrials2=subset(upperTrials2,upperTrials2$Trial<=i+trialWindow)
corr2 = (nrow(subset(theseTrials2, theseTrials2$correct=="CORRECT")))
wrong2 = (nrow(subset(theseTrials2,theseTrials2$correct=="WRONG")))
missing2 = (nrow(subset(theseTrials2,theseTrials2$correct=="MISSING")))
nonword2 = (nrow(subset(theseTrials2,theseTrials2$correct=="NONWORD")))
pc2[i] = corr2/(corr2+wrong2+nonword2+missing2)
xvalues[i]=i;
}
# compute the differences between the two groups
diffs=pc1-pc2
#set up some graphical paramets
par(mar=c(3,3,0.3,3),
fin=c(4.1,2.9),
fig=c(0,1,0,1),
mgp=c(2,1,0),
#usr=c(min(xdata),max(xdata),0,1),
lab=c(10,10,3),
#ps=11,
family='serif'
)
# plot the first group
plot(xvalues,pc1,col='#EE0000',pch="+",ylim=c(.4,.65),axes=F,xlab='trial window', ylab='proportion correct',cex=.8)
# add the second group
points(xvalues,pc2,col='#0000EE',pch="x",ylim=c(.4,.65),cex=.8)
# how much to add to the fitline y values to make it fit on the plot
fudge=.42
# plot the slope of the line we fitted to the differences
fitline=lsfit(xvalues,diffs+fudge)
abline(fitline,col='#000F00')
# do a correlation test
slope=cor.test(xvalues,diffs);
# the stats to put on the graph
stats=paste('r= ', formatC(slope$estimate,digits=3), ", p ", prettyPval(slope$p.value),sep='')
# put the values on the axes
axis(1) #bottom
axis2values=c(.4,.45,.5,.55,.6,.65)
axis(2,at=axis2values) #left
axis(4,at=axis2values,lab=formatC(axis2values-fudge,digits=2)) #right
box() #draw a box around the figure
# put the stats on the graph
text(280,.53,stats,cex=.75)
#label the right axis
mtext("difference between groups", side=4,line=2)
# label the points
text(80,.64,"fixed babble",col='#EE0000', cex=.8)
text(275,.44,"random babble",col='#0000EE', cex=.8)
# print as pdf
dev.print(pdf,"fig/learning.pdf",height=2.8,width=4,pointsize=11)