# R code for Bigham et al. (2009) bysubj LOT/THOUGHT study, from Gorman and
# Johnson in press

source('ky.R')
library(nlme)
library(plyr)
library(ggplot2)
options(contrasts = c('contr.sum', 'contr.poly'))

austin <- read.csv('austin.csv')

# average over subjects
bysubj <- ddply(austin, .(Speaker, Vowel), summarize, F1=median(F1.MID), F2=median(F2.MID))

# compute direction
mean(subset(bysubj, Vowel=='LOT')$F1)
mean(subset(bysubj, Vowel=='THOUGHT')$F1)
mean(subset(bysubj, Vowel=='LOT')$F2)
mean(subset(bysubj, Vowel=='THOUGHT')$F2)
median(subset(bysubj, Vowel=='LOT')$F1)
median(subset(bysubj, Vowel=='THOUGHT')$F1)
median(subset(bysubj, Vowel=='LOT')$F2)
median(subset(bysubj, Vowel=='THOUGHT')$F2)

# run paired tests 
t.test(F1 ~ Vowel, paired=TRUE, var.equals=FALSE, data=bysubj)$p.value
t.test(F2 ~ Vowel, paired=TRUE, var.equals=FALSE, data=bysubj)$p.value
wilcox.test(F1 ~ Vowel, paired=TRUE, data=bysubj)$p.value
wilcox.test(F2 ~ Vowel, paired=TRUE, data=bysubj)$p.value

# fixed-effects model (sanity check only!)
m0 <- lm(F2.MID ~ Vowel + Voice + Sex + Vowel:Sex + Vowel:Ethnicity + Vowel:Age, data=austin)
summary(m0)

# heteroscedastic mixed-effects model
m1 <- lme(F2.MID ~ Vowel + Ethnicity + Vowel:Ethnicity + Age + Age:Vowel + Voice + Sex, random=~1 + Vowel|Speaker, weights=varIdent(form=~1|Vowel), data=austin)
summary(m1)
anova(m1)

# sanity check
d <- subset(austin, Ethnicity == 'Anglo')
mean(subset(d, Vowel=='LOT')$F2.MID) - mean(subset(d, Vowel=='THOUGHT')$F2.MID)
d <- subset(austin, Ethnicity == 'Hispanic')
mean(subset(d, Vowel=='LOT')$F2.MID) - mean(subset(d, Vowel=='THOUGHT')$F2.MID)
d <- subset(austin, Ethnicity == 'AfrAme')
mean(subset(d, Vowel=='LOT')$F2.MID) - mean(subset(d, Vowel=='THOUGHT')$F2.MID)

# Compute Pillai score
bysubj <- subset(ddply(austin, .(Speaker), transform, count=length(F1.MID), Pillai=summary(manova(matrix(c(F1.MID, F2.MID), ncol=2) ~ Vowel))$stats[1, 2]), count > 7)

# plot the three most extreme speakers
png('bysubj_pil.png', width=7, height=6.5, units='in', res=600)
#special <- drop(subset(bysubj, Speaker %in% c('ZW52', 'FA06', 'MW40', 'QH54', 'SE43', 'KU78')))
special <- drop(subset(bysubj, Speaker %in% c('ZW52', 'MW40', 'QH54', 'KU78')))
suecial$Speaker <- paste('Pillai = ', round(special$Pillai, 3), sep='')
special.pillai <- ddply(special, .(Speaker, Vowel), summarise, F1.MID=median(F1.MID), F2.MID=median(F2.MID))
qplot(F2.MID, F1.MID, shape=Vowel, point='geom', alpha=I(2/3), main='Low back vowels, Austin', data=special, xlab='F2 (Hz)', ylab='F1 (Hz)') + facet_wrap(~Speaker) + scale_shape_manual(value=c(1, 4)) + theme_bw() + scale_x_reverse() + scale_y_reverse() 
#+ geom_line(data=special.pillai, aes(group=3), colour='black') + geom_point(data=special.pillai, aes(group=3), colour='black', size=2)
dev.off()