annotate latex/scripts/results_and_figures.R @ 4:3e666d0329b5 tip

updated code

cat("\n____________________________\n")
cat("\n____ MAIN PAPER RESULTS ____\n")
cat("\n____________________________\n\n")
cat("This script assumes that the script data_import.R has been run.\n")
cat("This script generates statistics and figures used in the JASA submission.\n\n")

##############################################################################
# METHOD SECTION OUTPUT
##############################################################################

cat("\n____ METHOD SECTION ____\n\n")
cat("- data for self-assessment table:\n")
print(with(meta[meta$instudy,], summary(data.frame(age, musical.background, sing.ability, sing.experience, choir.experience))))

cat("\n____ ACCURACY SECTION ____\n\n")

cat("- Interval error significance tests:\n\n")
significance.level <- 0.0001
res <- with(data[data$ispoor == F, ],
    t.test(intervalerror[ expectedinterval == 8 ]))
cat(        "  - Upward minor sixth is flat.\n")
cat(sprintf("    Is the difference significant at significance level %0.4f? %s (p = %g, df=%g, t = %g)\n\n", significance.level, res$p.value < significance.level, res$p.value, res$parameter, res$statistic))

res <- with(data[data$ispoor == F, ],
    t.test(intervalerror[ expectedinterval == 12 ]))
cat(        "  - Upward octave is flat.\n")
cat(sprintf("    Is the difference significant at significance level %0.4f? %s (p = %g, df=%g, t = %g)\n\n", significance.level, res$p.value < significance.level, res$p.value, res$parameter, res$statistic))

res <- with(data[data$ispoor == F, ],
    t.test(intervalerror[ expectedinterval == 0 ]))
cat(sprintf("  - Prime (0 interval) is sharp by %0.2f.\n", res$estimate))
cat(sprintf("    Is the difference significant at significance level %0.4f? %s (p = %g, df=%g, t = %g)\n\n", significance.level, res$p.value < significance.level, res$p.value, res$parameter, res$statistic))

cat("- explaining pitch error by following interval or beginning of phrase\n\n")
is.beginning.of.phrase <- I(data$noteid[!data$ispoor] %in% c(1,7,13,20))
res <- with(data[!data$ispoor, ],
    summary(lm(pitcherror ~ is.beginning.of.phrase + expectedfutureinterval)))
cat(sprintf("  - p: %0.3g, %0.3g\n",res$coefficients[3,4],res$coefficients[2,4]))
cat(sprintf("  - F = %g\n",res$fstatistic[1]))
cat(sprintf("  - DF: %0.0f\n",res$fstatistic[3]))
cat(sprintf("  - Variance explained: %0.1f%s\n",100*res$r.squared,"%"))
cat(sprintf("  - Sharpening in cents per note in interval: %0.1f\n\n",100*res$coefficients[3,1]))

res <- with(data[!data$ispoor, ],
    summary(lm(pitcherror ~ is.beginning.of.phrase)))
cat(sprintf("  - Variance explained only by intercept and following interval: %0.1f%s\n",100*res$r.squared,"%"))
cat(sprintf("  - p: %0.3g\n",res$coefficients[2,4]))
cat(sprintf("  - F = %g\n",res$fstatistic[1]))
cat(sprintf("  - DF: %0.0f\n\n",res$fstatistic[3]))
# cat(sprintf("  - Variance explained: %0.1f%s\n\n",100*res$r.squared,"%"))

res <- with(data[!data$ispoor, ],
    summary(lm(pitcherror ~ expectedfutureinterval)))
cat(sprintf("  - Variance explained only by intercept and beginning of phrase: %0.1f%s\n",100*res$r.squared,"%"))
cat(sprintf("  - p: %0.3g\n",res$coefficients[2,4]))
cat(sprintf("  - F = %g\n",res$fstatistic[1]))
cat(sprintf("  - DF: %0.0f\n\n",res$fstatistic[3]))
# cat(sprintf("  - Variance explained: %0.1f%s\n\n",100*res$r.squared,"%"))

et.ji.test <- with(data[!data$ispoor,], t.test(abs(pitcherror), abs(jipitcherror), paired=T))
etMAPE <- with(data[!data$ispoor,], mean(abs(pitcherror), na.rm=T))
jiMAPE <- with(data[!data$ispoor,], mean(abs(jipitcherror), na.rm=T))
cat("- Error difference between ET and JI error:\n")
cat("  Is it significant? Paired t-test says:")
cat(sprintf("%s (p = %g, df: %0.0f, t = %g)\n\n", et.ji.test$p.value < 0.01, et.ji.test$p.value, et.ji.test$parameter, et.ji.test$statistic), "\n")
cat(sprintf("ET (MAPE: %g) fits better than JI (MAPE: %g), by %g\n\n", etMAPE, jiMAPE, et.ji.test$estimate))

##############################################################################
# RESULTS SECTION OUTPUT
##############################################################################
cat("\n____ RESULTS SECTION ____\n\n")
cat("Basic metrics statistics.\n")
cat(sprintf("  - MAPE. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.mape$mape), median(agg.mape$mape), sd(agg.mape$mape)))
cat(sprintf("  - MAIE. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.maie$maie), median(agg.maie$maie), sd(agg.maie$maie)))
# cat(sprintf("  - IB. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.ib$ib), median(agg.ib$ib), sd(agg.ib$ib)))
cat(sprintf("  - D. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", mean(agg.drifts$diffs), median(agg.drifts$diffs), sd(agg.drifts$diffs)))
cat(sprintf("  - |D|. Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n\n", mean(abs(agg.drifts$diffs)), median(abs(agg.drifts$diffs)), sd(abs(agg.drifts$diffs))))
cat(sprintf("  - lin D. in cents: Mean: %0.3f;  median: %0.3f;  sd: %0.3f\n", 100*mean(linear.drift$lineardrift), 100*median(linear.drift$lineardrift), 100*sd(linear.drift$lineardrift)))

cat(sprintf("- Significance of drift:"))
cat(sprintf("  - significant drifts: %i of %i (%0.0f%s)\n", sum(agg.drifts.p$p<0.01), length(agg.drifts.p$p), 100*mean(agg.drifts.p$p<0.01), "%"))
with(agg.drifts.p[agg.drifts$diffs>0,], cat(sprintf("  - number of upward drifts: %i, of which %i significantly so.\n", length(p), sum(p<0.01))))
with(agg.drifts.p[agg.drifts$diffs<0,], cat(sprintf("  - number of downward drifts: %i, of which %i significantly so.\n\n", length(p), sum(p<0.01))))

cat(sprintf("- Significance of linear drift"))
cat(sprintf("  - significant drifts: %i of %i (%0.0f%%) at conf. level 0.01, %i of %i (%0.0f%%) at conf. level 0.05 \n", sum(linear.drift$p.value<0.01), length(linear.drift$p.value), 100*mean(linear.drift$p.value<0.01), sum(linear.drift$p.value<0.05), length(linear.drift$p.value), 100*mean(linear.drift$p.value<0.05)))
with(linear.drift[linear.drift$lineardrift>0,], cat(sprintf("  - number of upward drifts: %i, of which %i significantly so at 0.01.\n", length(p.value), sum(p.value<0.01))))
with(linear.drift[linear.drift$lineardrift<0,], cat(sprintf("  - number of downward drifts: %i, of which %i significantly so at 0.01.\n\n", length(p.value), sum(p.value<0.01))))
with(linear.drift[linear.drift$lineardrift>0,], cat(sprintf("  - number of upward drifts: %i, of which %i significantly so at 0.05.\n", length(p.value), sum(p.value<0.05))))
with(linear.drift[linear.drift$lineardrift<0,], cat(sprintf("  - number of downward drifts: %i, of which %i significantly so at 0.05.\n\n", length(p.value), sum(p.value<0.05))))


cat("- Correlation table (LaTeX syntax)\n")
# cor.df <- data.frame(sg.abl=as.numeric(merged.agg$sing.ability), sg.exp=as.numeric(merged.agg$sing.experience), mus.bg=as.numeric(merged.agg$musical.background), ch.exp=as.numeric(merged.agg$choir.experience), MAIE=merged.agg$maie, MAPE=merged.agg$mape, IB=merged.agg$ib, D=merged.agg$diffs, abs.D=abs(merged.agg$diffs))
cor.df <- data.frame(sg.abl=as.numeric(merged.agg$sing.ability), sg.exp=as.numeric(merged.agg$sing.experience), mus.bg=as.numeric(merged.agg$musical.background), ch.exp=as.numeric(merged.agg$choir.experience), MAIE=merged.agg$maie, MAPE=merged.agg$mape, lin.D=merged.agg$lineardrift, abs.D=abs(merged.agg$diffs), D=merged.agg$diffs)
# cortable(data.frame(sg.abl=as.numeric(merged.agg$sing.ability), mus.bg=as.numeric(merged.agg$musical.background), ch.exp=as.numeric(merged.agg$choir.experience), MAIE=merged.agg$maie, ib=merged.agg$ib, MAPE=merged.agg$mape, D=merged.agg$diffs, abs.D=abs(merged.agg$diffs)))
cortable(cor.df)
cat("\n\n")

# testing singer predicting MAPE, MAIE, and drift
cat("- Variance explained by singer for different measures\n\n")
res <- summary(lm(mape~singer, data=agg.mape))
pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
df <- res$fstatistic[2]
cat(sprintf("  - MAPE: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared*100, pval, df, res$fstatistic[1]))
res <- summary(lm(maie~singer, data=agg.maie))
pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
df <- res$fstatistic[2]
cat(sprintf("  - MAIE: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared * 100, pval, df, res$fstatistic[1]))
# res <- summary(lm(ib~singer, data=agg.ib))
# pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
# cat(sprintf("  - IB: %0.2f%%, p-value: %0.5g\n", res$r.squared * 100, pval))
res <- summary(lm(lineardrift~singer, data=linear.drift))
pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
df <- res$fstatistic[2]
cat(sprintf("  - lin.D: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared * 100, pval, df, res$fstatistic[1]))
res <- summary(lm(diffs~singer, data=agg.drifts))
pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
df <- res$fstatistic[2]
cat(sprintf("  - D: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n", res$r.squared * 100, pval, df, res$fstatistic[1]))
res <- summary(lm(abs(diffs)~singer, data=agg.drifts))
pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
df <- res$fstatistic[2]
cat(sprintf("  - |D|: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n\n", res$r.squared * 100, pval, df, res$fstatistic[1]))


# testing influence of condition
cat("- Testing influence of singing condition (Normal/Masked/Imagined) via Kruskal-Wallis tests (sign. level 0.01):\n\n")
res <- kruskal.test(mape ~ mode, data = agg.mape)
cat(sprintf("  - MAPE significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
res <- kruskal.test(maie ~ mode, data = agg.maie)
cat(sprintf("  - MAIE significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
# res <- kruskal.test(ib ~ mode, data = agg.ib)
# cat(sprintf("  - IB significant difference due to condition? %s (p-value: %0.5g)\n", res$p.value<0.01, res$p.value))
res <- kruskal.test(lineardrift ~ mode, data = linear.drift)
cat(sprintf("  - lin. D significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
res <- kruskal.test(diffs ~ mode, data = agg.drifts)
cat(sprintf("  - D significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))
res <- kruskal.test(abs(diffs) ~ mode, data = agg.drifts)
cat(sprintf("  - |D| significant difference due to condition? %s (p-value: %0.5g, df = %g, chisquared = %g)\n", res$p.value<0.01, res$p.value, res$parameter, res$statistic))

# cat("- \"Flat\" singer: \n\n")
# "hl"
# with(agg.drifts.all, cat(sprintf("  - drift in normal condition: %0.2f semitones\n", diffs[singer=="hl" & mode=="normal"])))
# with(linear.drift.all, cat(sprintf("  - linear drift in normal condition: %0.2f cents per note\n", 100 * lineardrift[singer=="hl" & mode=="normal"])))
# with(agg.drifts.all, cat(sprintf("  - drift in masked condition: %0.2f semitones\n", diffs[singer=="hl" & mode=="masked"])))
# with(linear.drift.all, cat(sprintf("  - linear drift in masked condition: %0.2f cents per note\n", 100 * lineardrift[singer=="hl" & mode=="masked"])))
# with(agg.drifts.all, cat(sprintf("  - drift in imagined condition: %0.2f semitones\n", diffs[singer=="hl" & mode=="imagined"])))
# with(linear.drift.all, cat(sprintf("  - linear drift in imagined condition: %0.2f cents per note\n\n", 100 * lineardrift[singer=="hl" & mode=="imagined"])))

##############################################################################
# MODEL SECTION OUTPUT
##############################################################################
cat("\n____ MODEL SECTION ____\n\n")

# notevar <- with(mdr[!mdr$ispoor,], aggregate(I(p-s-tstart), list(notenumber=notenumber), function(x)var(x,na.rm=T))[,2])
# notevar2 <- with(data[!data$ispoor,], aggregate(pitcherror, list(notenumber=notenumber), function(x)var(x,na.rm=T))[,2])
# driftvar <- notevar-notevar2
# driftvar.fit <- lm(driftvar[7:75]~I(7:75))
# sqrt(driftvar.fit$coefficients[2]*50)

cat("- existence of 'spread'\n")
var.by.notenumber <- aggregate(data$relativetonic[!data$ispoor], list(data$notenumber[!data$ispoor]), function(x) {var(x, na.rm=T)})
var2.by.notenumber <- aggregate(data$pitcherror[!data$ispoor], list(data$notenumber[!data$ispoor]), function(x) {var(x, na.rm=T)})
counts.by.notenumber <- aggregate(data$relativetonic[!data$ispoor], list(data$notenumber[!data$ispoor]), function(x) {sum(!is.na(x))})
dist.by.notenumber <- aggregate(abs(data$relativetonic[!data$ispoor]), list(data$notenumber[!data$ispoor]), function(x) {mean(x, na.rm=T)})
fit1 <- lm(var.by.notenumber[,2]~var2.by.notenumber[,2], weights=counts.by.notenumber[,2])
fit2 <- lm(fit1$residuals~var.by.notenumber[(burnin+1):75,1], weights=counts.by.notenumber[(burnin+1):75,2])
summary(fit2)
res <- summary(fit2)
pval <- pf(res$fstatistic[1],res$fstatistic[2],res$fstatistic[3],lower.tail=FALSE)
cat(sprintf("  - notenumber: %0.2f%%, p-value: %0.5g, df = %g, F = %g\n\n", res$r.squared * 100, pval, res$fstatistic[1], res$fstatistic[3]))
plot(names(fit1$residuals), fit1$residuals)
abline(fit2)

intvar <- var(data$intervalerror[!data$ispoor], na.rm=T)
expected.D.var <- 50 * intvar
expected.D.sd <- expected.D.var^0.5
pitch.err.sd <- sd(data$pitcherror[!data$ispoor], na.rm=T)
cat(sprintf("- inadequacy of mu = 0 model:\n\n"))
cat(sprintf("  - interval (error) variance: %0.3f\n", intvar))
cat(sprintf("  - expected variance of D (interval variance summed over 50 notes): %0.3f\n", expected.D.var))
cat(sprintf("  - ... and the corresponding standard deviation: %0.3f\n\n", expected.D.sd))
cat(sprintf("  - In comparison the standard deviation of pitch erros is: %0.3f\n\n", pitch.err.sd))

model.mape <- apply(md[!md$ispoor,(1:nmu)+8], 2, function(x) {mean(abs(x),na.rm=T)})
model.mape.singerwise <- aggregate(md[!md$ispoor,(1:nmu)+8], list(md$singer[!md$ispoor]), function(x) {mean(abs(x),na.rm=T)})
model.mape.recordingwise <- aggregate(md[!md$ispoor,(1:nmu)+8], list(md$singer[!md$ispoor],md$mode[!md$ispoor]), function(x) {mean(abs(x),na.rm=T)})

recording.mus <- data.frame(singer = model.mape.recordingwise[,1], mu = mus[apply(model.mape.recordingwise[,-(1:2)], 1, which.min)])
# singer.mus <- data.frame(singer = model.mape.singerwise[,1], mu = mus[apply(model.mape.singerwise[,-1], 1, which.min)])
singer.mus <- aggregate(recording.mus[,2], list(recording.mus[,1]),mean)

cat(sprintf("- Better parameter for mu, firstly via 'model MAPE':\n\n"))
cat(sprintf("  - the smallest error of %0.2f is achieved with model '%s'\n", min(model.mape), names(which.min(model.mape))))
cat(sprintf("  - error at mu = 0: %0.2f\n", model.mape[1]))
cat(sprintf("  - error at mu = 1: %0.2f\n\n", model.mape[length(model.mape)]))

#-------------------------------
mdr <- merge(md, mr[,c(1,2,5)], by=c("singer", "mode"))
# plot(I(p-s-tstart)~factor(notenumber), data=mdr)

sig2 <- var(data$pitcherror[!data$ispoor], na.rm=T)
sig2int <- var(data$intervalerror[!data$ispoor], na.rm=T)
n <- 50
# vd <- var(drifts[!drifts$singer %in% poor.singers,4], na.rm=T)
# vd <- var(drifts$diffs[!drifts$singer %in% poor.singers], na.rm=T)
vd <- var(agg.drifts$diffs)
est.mu <- 1-sqrt(vd/(n*sig2))
  
cat("- Estimate according to drift variance formula:\n\n")
cat(sprintf("  - overall pitch error variance: %0.3f\n", sig2))
# cat(sprintf("  - overall interval error variance: %0.3f", sig2int))
cat(sprintf("  - observed variance over 50-note gap: %0.3f\n", vd))
cat(sprintf("  - estimated mu based on these: %0.3f\n\n", est.mu))

cat("- Stats for mu by singer:\n")
cat(sprintf("  - range: %0.3f to %0.3f\n", range(singer.mus[,2])[1], range(singer.mus[,2])[2]))
cat(sprintf("  - mean: %0.3f, median: %0.3f, sd: %0.3f\n\n", mean(singer.mus[,2]), median(singer.mus[,2]), sd(singer.mus[,2])))


# n <- 50
# mu <- 0.916
# sig2 <- var(data$pitcherror[ !data$ispoor ], na.rm=T)
# vv <- var(drifts$diffs[!drifts$singer %in% poor.singers], na.rm=T)
# vv <- var(drifts$diffs[!drifts$singer %in% poor.singers & drifts$noteid ==iNote], na.rm=T)

# n <- 50
# sig2 <- var(data$pitcherror[ !data$ispoor ], na.rm=T)
# vv <- var(agg.drifts$diffs)
# est.mu <- 1-sqrt(vv/(n*sig2))
# 
# vv <- var(agg.drifts$diffs)
# for (mu in c(0,0.8,0.9,est.mu,  0.95, 1))
# {
#     dvar <- (1-mu)^2 * n * sig2
#     cat(sprintf("Example: mu = %g, sigma^2 = %g, n = %g\n", mu, sig2, n))
#     cat(sprintf(" ... var = %g, sd = %g\n", dvar, sqrt(dvar)))
# }
# CHECK!!!!!

# singermodemu <- c()
# for (singer in singers[!singers %in% poor.singers])
# {
#     for (mode in c("normal","imagined","masked"))
#     {
#         sig2 <- var(data$pitcherror[drifts$singer==singer & drifts$mode==mode], na.rm=T)
#         # sig2int <- var(data$intervalerror[data$singermode==singermode], na.rm=T)
#         n <- 50
#         vd <- var(drifts$diffs[drifts$singer==singer & drifts$mode==mode], na.rm=T)
#         est.mu <- 1-sqrt(vd/(n*sig2))
#         # cat(sprintf("%s: %0.3f\n", singer, est.mu))
#         singermode <- sprintf("%s%s", singer, mode)
#         singermodemu <- rbind(singermodemu, data.frame(singermode=singermode, singer=singer, mode=mode, mu=est.mu))
#     }
# }
# hist(singermodemu$mu, col="gray")

# make drift figures
for (singermode in sort(unique(data$singermode)))
{
	pdf(sprintf("../figures/drift_figure_%s.pdf", singermode), width = 3.5, height = 3.5)
	makeDriftGraph(data, singermode)
	dev.off()
}

for (singermode in sort(unique(data$singermode)))
{
	pdf(sprintf("../figures/drift_figure_model0850_%s.pdf", singermode), width = 5.5, height = 4.5)
	with(data[data$singermode==singermode,], plot(notenumber, newmedian-semitone, pch="+", ylab="pitch", xlab="note"))
	with(md[md$singermode==singermode,], points(notenumber, p-s-err0850, pch=20, type='o'))
	dev.off()
}

countsByInterval <- with(data[!data$ispoor,], aggregate(data.frame(count=data$intervalerror), list(interval=data$expectedinterval), function(x){sum(!is.na(x))})) 

pdf("../figures/interval_error_by_interval_boxplot.pdf", width = 5.5, height = 4)
par(mar=c(2,4,1,.5))
plot(c(1,18), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
abline(h=0)
abline(h=-1, lty=2)
abline(h=1, lty=2)
plot(intervalerror ~ factor(expectedinterval, levels=-5:12), data=data, subset = !ispoor,
	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, xaxt='n', yaxt='n', add=T,
	ylab="interval error (with respect to ET)", xlab="interval in semitones",
	cex.axis=0.8, lty=1)
posi <- c(1:6, 8, 11, 13, 14, 18)
axis(2, at=c(-1,0,1), cex.axis=0.8)
axis(1, at=posi, labels=posi-6, cex.axis=0.8)
# text(counts$interval + 6, rep(2,11), counts$count, cex=0.7)
mtext(c("n=", countsByInterval$count), at=c(-.1,countsByInterval$interval + 6), cex=0.7)
dev.off()

countsByNoteID <- with(data[!data$ispoor,], aggregate(data.frame(count=data$intervalerror), list(noteid=data$noteid), function(x){sum(!is.na(x))}))

pdf("../figures/interval_error_by_notenumber_boxplot.pdf", width = 7, height = 4)
par(mar=c(2,.5,1,1))	
plot(c(1,25), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
abline(h=0)
abline(h=-1, lty=2)
abline(h=1, lty=2)
plot(intervalerror ~ factor(noteid, levels=1:75), data=data, subset = !ispoor,
	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, add=T,
	ylab="interval error (with respect to ET)", xlab="interval in semitones",
	cex.axis=0.8, lty=1, yaxt='n', boxwex=0.6)
# posi <- c(1:6, 8, 11, 13, 14, 18)
# axis(1, at=posi, labels=posi-6, cex.axis=0.8)
# text(countsByNoteID$noteid, rep(2,11), countsByNoteID$count, cex=0.7)
mtext(c(countsByNoteID$count), at=c(countsByNoteID$noteid), cex=0.7)
text(1:25, rep(-1.5,11), c("(-5)", gtdata$interval_ST[-25]), cex=0.7)
text(2.7, -1.35 , "intervals in semitones", cex=0.7)
text(1-.5, 1.3, "phrase 1", pos=4)
abline(v=6.5, lty=1, col='gray')
text(7-.5, 1.3, "phrase 2", pos=4)
abline(v=12.5, lty=1, col='gray')
text(13-.5, 1.3, "phrase 3", pos=4)
abline(v=19.5, lty=1, col='gray')
text(20-.5, 1.3, "phrase 4", pos=4)
dev.off()

pdf("../figures/pitch_error_by_interval_boxplot.pdf", width = 5.5, height = 4)
par(mar=c(2,4,1,.5))
plot(c(1,18), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
abline(h=0)
abline(h=-1, lty=2)
abline(h=1, lty=2)
plot(pitcherror ~ factor(expectedinterval, levels=-5:12), data=data, subset = !ispoor,
	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, xaxt='n', yaxt='n', add=T,
	ylab="pitch error (with respect to ET)", xlab="interval in semitones",
	cex.axis=0.8, lty=1)
posi <- c(1:6, 8, 11, 13, 14, 18)
axis(2, at=c(-1,0,1), cex.axis=0.8)
axis(1, at=posi, labels=posi-6, cex.axis=0.8)
# text(counts$interval + 6, rep(2,11), counts$count, cex=0.7)
# mtext(c("n=", countsByInterval$count), at=c(-.1,countsByInterval$interval + 6), cex=0.7)
dev.off()

countsBySemitone <- with(data[!data$ispoor,], aggregate(data.frame(count=data$pitcherror), list(interval=data$semitone), function(x){sum(!is.na(x))})) 
pdf("../figures/pitch_error_by_semitone_boxplot.pdf", width = 4.5, height = 4)
par(mar=c(2,4,1,.5))
plot(c(.5,13.5), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
abline(h=0)
abline(h=-1, lty=2)
abline(h=1, lty=2)
plot(pitcherror ~ factor(semitone, levels=-5:7), data=data, subset = !ispoor,
	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, xaxt='n', yaxt='n', add=T,
	ylab="pitch error (with respect to ET)", xlab="interval in semitones",
	cex.axis=0.8, lty=1)
posi <- c(1, 3, 5, 6, 8, 10, 11, 13)
axis(2, at=c(-1,0,1), cex.axis=0.8)
axis(1, at=posi, labels=posi-6, cex.axis=0.8)
# text(counts$interval + 6, rep(2,11), counts$count, cex=0.7)
# mtext(c("n=", countsByInterval$count), at=c(-.1,countsByInterval$interval + 6), cex=0.7)
dev.off()

pdf("../figures/pitch_error_by_notenumber_boxplot.pdf", width = 7, height = 4)
par(mar=c(2,.5,1,1))	
plot(c(1,25), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
abline(h=0)
abline(h=-1, lty=2)
abline(h=1, lty=2)
plot(pitcherror ~ factor(noteid, levels=1:75), data=data, subset = !ispoor,
	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, add=T,
	ylab="pitch error (with respect to ET)", xlab="interval in semitones",
	cex.axis=0.8, lty=1, yaxt='n', boxwex=0.6)
# posi <- c(1:6, 8, 11, 13, 14, 18)
# axis(1, at=posi, labels=posi-6, cex.axis=0.8)
# text(countsByNoteID$noteid, rep(2,11), countsByNoteID$count, cex=0.7)
# mtext(c(countsByNoteID$count), at=c(countsByNoteID$noteid), cex=0.7)
text(1:25, rep(-1.5,11), gtdata$pitch_ion_cents_et/100, cex=0.7)
text(2.7, -1.35 , "pitch relative to tonic", cex=0.7)
text(1-.5, 1.3, "phrase 1", pos=4)
abline(v=6.5, lty=1, col='gray')
text(7-.5, 1.3, "phrase 2", pos=4)
abline(v=12.5, lty=1, col='gray')
text(13-.5, 1.3, "phrase 3", pos=4)
abline(v=19.5, lty=1, col='gray')
text(20-.5, 1.3, "phrase 4", pos=4)
dev.off()

# linear fit example plot
pdf("../figures/pitch_error_example.pdf", width = 6, height = 5)
par(mar=c(4,4,.1,.1))
ind <- which(data$singer == "jp" & data$mode == "imagined" & data$run == 1)
plot(data$newmedian[ind], pch = 5,
	xlab="note number", ylab="pitch",
	xlim = c(0,25), las = 1)
points((data$newmedian-data$semitone)[ind], pch=18)
fit <- lm(I(newmedian-semitone)~notenumber, data = data[ind,])
for (ii in 1:25)
{
	lines(c(ii,ii), c(fit$fitted.values[ii], (data$newmedian-data$semitone)[ind[ii]]))
}
abline(a=fit$coefficients[1], b=fit$coefficients[2])
text(0, fit$coefficients[1]+0.3, 0, cex = 0.7)
for (offs in 1:8)
{
	abline(a=fit$coefficients[1]+offs, b=fit$coefficients[2], lty = 2)
	abline(a=fit$coefficients[1]-offs, b=fit$coefficients[2], lty = 2)
	text(0, fit$coefficients[1]+offs+0.3, +offs, cex = 0.7)
	text(0, fit$coefficients[1]-offs+0.3, -offs, cex = 0.7)
}
dev.off()

# drift histogram
pdf("../figures/drift_histogram.pdf", width = 6, height = 4)
par(mar=c(4,4,.3,.5))
hist(agg.drifts$diffs, seq(-1.5,1.5,by=0.2), col="dark gray", border="black", 
	xlab = "drift", ylab = "number of recordings", main = "")
dev.off()

# mean absolute pitch error histogram
pdf("../figures/mape_histogram.pdf", width = 3.3, height = 4)
par(mar=c(4,4,.3,.5))
hist(agg.mape$mape, seq(0,0.4,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24))
abline(h=5, col = "gray")
abline(h=10, col = "gray")
abline(h=15, col = "gray")
abline(h=20, col = "gray")
hist(agg.mape$mape, seq(0,0.4,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), add = T)
dev.off()

# mean absolute interval error histogram
pdf("../figures/maie_histogram.pdf", width = 3.3, height = 4)
par(mar=c(4,1,.3,.5))
hist(agg.maie$maie, seq(0,.5,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n")
abline(h=5, col = "gray")
abline(h=10, col = "gray")
abline(h=15, col = "gray")
abline(h=20, col = "gray")
hist(agg.maie$maie, seq(0,.5,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n", add = T)
dev.off()
mean(agg.maie$maie)
median(agg.maie$maie)
sd(agg.maie$maie)

# MAR
pdf("../figures/mar_histogram.pdf", width = 3.3, height = 4)
par(mar=c(4,1,.3,.5))
hist(agg.rise$rise, seq(0,.5,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n")
abline(h=5, col = "gray")
abline(h=10, col = "gray")
abline(h=15, col = "gray")
abline(h=20, col = "gray")
hist(agg.rise$rise, seq(0,.5,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n", add = T)
dev.off()
mean(agg.rise$rise)
median(agg.rise$rise)
sd(agg.rise$rise)

# drift
pdf("../figures/d_histogram.pdf", width = 3.3, height = 4)
par(mar=c(4,1,.3,.5))
hist(agg.drifts$diffs, seq(-.6,.6,by=0.1), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n")
abline(h=5, col = "gray")
abline(h=10, col = "gray")
abline(h=15, col = "gray")
abline(h=20, col = "gray")
hist(agg.drifts$diffs, seq(-.6,.6,by=0.1), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n", add = T)
dev.off()

# absolute drift
pdf("../figures/abs_d_histogram.pdf", width = 3.3, height = 4)
par(mar=c(4,1,.3,.5))
hist(abs(agg.drifts$diffs), seq(0,.6,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n")
abline(h=5, col = "gray")
abline(h=10, col = "gray")
abline(h=15, col = "gray")
abline(h=20, col = "gray")
hist(abs(agg.drifts$diffs), seq(0,.6,by=0.05), col="dark gray", border="black", 
	xlab = "semitones", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n", add = T)
dev.off()

# linear drift
pdf("../figures/lineard_histogram.pdf", width = 3.3, height = 4)
par(mar=c(4,1,.3,.5))
hist(100*linear.drift$lineardrift, seq(-1.2,1.2,by=0.2), col="dark gray", border="black", 
	xlab = "cents", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n")
abline(h=5, col = "gray")
abline(h=10, col = "gray")
abline(h=15, col = "gray")
abline(h=20, col = "gray")
hist(100*linear.drift$lineardrift, seq(-1.2,1.2,by=0.2), col="dark gray", border="black", 
	xlab = "cents", ylab = "number of recordings", main = "",
	ylim=c(0,24), yaxt="n", add = T)
dev.off()

# mu
pdf("../figures/mu_histogram.pdf", width = 4.3, height = 3)
par(mar=c(4,4,.3,.5))
hist(singer.mus[,2], seq(0.6, 1.0, 0.4/5), col="dark gray", border="black", 
	xlab = expression(mu), ylab = "number of singers", main = "",
	add = F)
dev.off()



# drift significance plot
pdf("../figures/drift_significance_plot.pdf", width = 6, height = 3.5)
plot(agg.drifts$diffs, agg.drifts.p$p, log="y", col = rgb(0,0,0,0.4), 
	pch=19, yaxt = "n", ylab = "p-value                  ", 
	xlab = "drift magnitude in semitones",
	xlim=c(-1,1))
axis(2, at = 10^(-(1:4)), label = c("0.1","0.01","0.001","0.0001"),las=1)
abline(h=0.01, lty = 2)
abline(v=0.0)
text(0.75,0.0000005,"significant\nupward\ndrift")
text(-0.75,0.0000005,"significant\ndownward\ndrift")
dev.off()

# drift significance plot
pdf("../figures/lineardrift_significance_plot.pdf", width = 6, height = 4.8)
plot(linear.drift$lineardrift, linear.drift$p.value, log="y", col = rgb(0,0,0,0.4), 
	pch=19, yaxt = "n", ylab = "p-value                  ", 
	xlab = "linear drift per note in semitones")
axis(2, at = 10^(-(1:4)), label = c("0.1","0.01","0.001","0.0001"),las=1)
abline(h=0.01, lty = 2)
abline(v=0.0)
text(0.0035,0.0000005,"significant\nupward\ndrift")
text(-0.0035,0.0000005,"significant\ndownward\ndrift")
dev.off()


# how many are significantly adrift?
# sum(agg.drifts.p$p<0.01)
# mean(agg.drifts.p$p<0.01)

# drift vs mape
pdf("../figures/drift_vs_mape.pdf", width = 5.5, height = 4)
par(mar=c(4,4,.3,.5))
outl <- agg.mape$mape>0.35
fit1 <- lm(agg.mape$mape~agg.drifts$diffs)
fit2 <- lm(agg.mape$mape~agg.drifts$diffs, subset=!outl)
plot(agg.drifts$diffs, agg.mape$mape, type="n", 
	xlab="drift magnitude in semitones",
	ylab="mean absolute pitch error")
points(agg.drifts$diffs[outl], agg.mape$mape[outl], pch=19)
points(agg.drifts$diffs[!outl], agg.mape$mape[!outl], pch=21)
abline(fit1, lty=2)
abline(fit2)
dev.off()

# drift vs maie
pdf("../figures/drift_vs_maie.pdf", width = 5.5, height = 4)
par(mar=c(4,4,.3,.5))
outl <- agg.mape$mape>0.35
fit1 <- lm(agg.maie$maie~agg.drifts$diffs)
fit2 <- lm(agg.maie$maie~agg.drifts$diffs, subset=!outl)
plot(agg.drifts$diffs, agg.maie$maie, type="n", 
	xlab="drift magnitude in semitones",
	ylab="mean absolute pitch error")
points(agg.drifts$diffs[outl], agg.maie$maie[outl], pch=19)
points(agg.drifts$diffs[!outl], agg.maie$maie[!outl], pch=21)
abline(fit1, lty=2)
abline(fit2)
dev.off()

# drift vs IB
pdf("../figures/drift_vs_ib.pdf", width = 5.5, height = 4)
par(mar=c(4,4,.3,.5))
fit1 <- lm(agg.ib$ib~agg.drifts$diffs)
plot(agg.drifts$diffs, agg.ib$ib, type="n", 
	xlab="drift magnitude in semitones",
	ylab="mean interval error")
points(agg.drifts$diffs, agg.ib$ib, pch=21)
abline(fit1, lty=2)
dev.off()


# cor(merged.singeragg$maie, merged.singeragg$mape)

pdf("../figures/mape_maie_drift.pdf", width = 5, height = 5)
# par(mar=c(1,1,1,1))
plot(data.frame(MAPE=agg.mape$mape, MAIE=agg.maie$maie, Drift=agg.drifts$diffs))
dev.off()
# and their correlations
cor.test(agg.mape$mape, agg.maie$maie)
cor.test(agg.mape$mape, agg.drifts$diffs)
cor.test(agg.maie$maie, agg.drifts$diffs)

# corrgram(data.frame(merged.singeragg$maie, merged.singeragg$mape, merged.singeragg$diffs, as.numeric(merged.singeragg$sing.ability), as.numeric(merged.singeragg$musical.background), as.numeric(merged.singeragg$choir.experience)))
# corrgram(data.frame(merged.agg$maie, merged.agg$mape, merged.agg$diffs, as.numeric(merged.agg$sing.ability), as.numeric(merged.agg$musical.background), as.numeric(merged.agg$choir.experience)))

# pdf("../figures/corrgram.pdf", width = 5, height = 5)
# corrgram(data.frame(MAIE=merged.agg$maie, ib=merged.agg$ib, MAPE=merged.agg$mape, D=merged.agg$diffs, singing.ability=as.numeric(merged.agg$sing.ability), musical.bg=as.numeric(merged.agg$musical.background), choir.exp=as.numeric(merged.agg$choir.experience)), lower.panel=panel.conf, upper.panel=panel.pie)
# dev.off()

# slope figure
pdf("../figures/rise_by_notenumber_boxplot.pdf", width = 5.5, height = 4)
par(mar=c(2,3.8,.1,.1))	
plot(c(1,25), c(-1.6,1.6), type='n', xaxt='n', xlab = '', ylab='', yaxt='n')
abline(h=0)
abline(h=-1, lty=2)
abline(h=1, lty=2)
plot(I(slope * duration * 44100/1024) ~ factor(noteid, levels=1:75), data=data, subset = !ispoor,
	pch=20, col=rgb(0.9,0.9,0.9), ylim = c(-2,2), outline=F, add=T,
	ylab="rise (semitones)", xlab="interval in semitones",
	cex.axis=0.8, lty=1, boxwex=0.6)
# posi <- c(1:6, 8, 11, 13, 14, 18)
# axis(1, at=posi, labels=posi-6, cex.axis=0.8)
# text(countsByNoteID$noteid, rep(2,11), countsByNoteID$count, cex=0.7)
mtext(c(countsByNoteID$count), at=c(countsByNoteID$noteid), cex=0.7)
# text(1:25, rep(-1.5,11), c("(-5)", gtdata$interval_ST[-25]), cex=0.7)
# text(2.7, -1.35 , "intervals in semitones", cex=0.7)
text(1-.5, 1.3, "phrase 1", pos=4)
abline(v=6.5, lty=1, col='gray')
text(7-.5, 1.3, "phrase 2", pos=4)
abline(v=12.5, lty=1, col='gray')
text(13-.5, 1.3, "phrase 3", pos=4)
abline(v=19.5, lty=1, col='gray')
text(20-.5, 1.3, "phrase 4", pos=4)
dev.off()

# mu versus model MAPE
pdf("../figures/memory_vs_model_mape.pdf", width = 4.3, height = 3)
par(mar=c(4,4,.3,.5))
plot(mus,
     apply(md[!md$ispoor,(1:nmu)+8], 2, function(x) {mean(abs(x),na.rm=T)}), 
     type='l', 
     pch=20,
     xlab=expression(mu),
     ylab="mean absolute pitch error",
     ylim=c(0.18,0.30))
abline(h=mean(abs(data$pitcherror[!data$ispoor]), na.rm=T), lty=2)
dev.off()