# R version: 3.3.2
library(lattice) # version: 0.20-34
library(ez) # version: 4.4-0

df = read.table("df.csv")

############
## Define panels with horizontal lines
#
panel0to1 = function (...) {
  for (i in seq(0, 1, 0.1)) {
    panel.abline(h=i, col = 'grey')
  }
  panel.bwplot(...)
}

panel0to50 = function (...) {
  for (i in seq(0, 50, 10)) {
    panel.abline(h=i, col='grey')
  }
  panel.bwplot(...)
}

panel0to500 = function (...) {
  for (i in seq(0, 500, 50)) {
    panel.abline(h=i, col = 'grey')
  }
  panel.bwplot(...)
}

######################
## Define data subsets
#

# V1 is o, o:, u, or u:
subset.ou = df$v1 %in% c("o", "o:", "u", "u:")

# the only "minimal pair" which differs as VC vs. V:C: but hast the same C1 and V2
subset.taka = df$word %in% c("taka", "taakka")


##################
## Proportional vowel duration
#

ylim = c(0,1)

# Grouped by category (V1 and C2 quantity), condition and age
bwplot   (data = df, propVowelDurationWithRelease ~ category | condition * age, panel = panel0to1, ylim = ylim)
aggregate(data = df, propVowelDurationWithRelease ~ category * condition * age, length)



#####
# Figure 1
# Restricted to o, o:, u, and u: words
#
# Export this plot as SVG, size 400x450
bwplot   (data = df[subset.ou,], propVowelDurationWithRelease ~ category | condition * age, panel = panel0to1, ylim = ylim, ylab = "Proportional vowel duration (PVD)")
aggregate(data = df[subset.ou,], propVowelDurationWithRelease ~ category * condition * age, length)
df.anova.pvd = aggregate (data = df[subset.ou,], propVowelDurationWithRelease ~ speaker * category * condition * age, FUN = mean)
ezANOVA (df.anova.pvd, .(propVowelDurationWithRelease), .(speaker), .(category, condition), between = .(age))


#####
# Figure 2
# Restricted to taka and taakka
#
# Export this plot as SVG, size 400x450 
bwplot   (data = df[subset.taka,], propVowelDurationWithRelease ~ category | condition * age, panel = panel0to1, ylim = c(0.3 , 0.7), ylab = "Proportional vowel duration (PVD)")
aggregate(data = df[subset.taka,], propVowelDurationWithRelease ~ category * condition * age, length)
df.anova.pvd.taka = aggregate (data = df[subset.taka,], propVowelDurationWithRelease ~ speaker * category * condition * age, FUN = mean)
ezANOVA (df.anova.pvd.taka, .(propVowelDurationWithRelease), .(speaker), .(category, condition), between = .(age))


###########################
## Geminate/singleton ratio
#

# build a data frame with these columns:
# speaker, age, condition, q_v1, q_c2, N, v1dur, c2dur_total
df.duration              = aggregate(data = df[subset.ou,], v1dur       ~ speaker * age * condition * q_v1 * q_c2, FUN = mean)
df.duration$c2dur_total  = aggregate(data = df[subset.ou,], c2dur_total ~ speaker * age * condition * q_v1 * q_c2, FUN = mean)$c2dur_total

# build a data frame with these columns:
# speaker, age, condition, q_c2, shortV1dur, longV1dur, V1ratio
df.duration.v1           = df.duration[df.duration$q_v1 == 'short', c(1:3, 5:6)]
df.duration.v1$longV1dur = df.duration[df.duration$q_v1 == 'long',  c(1:3, 5:6)]$v1dur
colnames(df.duration.v1)[5] = "shortV1dur"
df.duration.v1$V1ratio = df.duration.v1$longV1dur / df.duration.v1$shortV1dur

# and another one for C2:
# spekaer, age, condition, q_v1, shortC2dur, longC2dur, C2ratio
df.duration.c2                 = df.duration[df.duration$q_c2 == 'short', c(1:4, 7)]
df.duration.c2$longC2dur_total = df.duration[df.duration$q_c2 == 'long',  c(1:4, 7)]$c2dur_total
colnames(df.duration.c2)[5] = "shortC2dur_total"
df.duration.c2$C2ratio = df.duration.c2$longC2dur_total / df.duration.c2$shortC2dur_total

# test V1 ratio
dotplot   (data = df.duration.v1, V1ratio ~ q_c2 | condition * age)
aggregate (data = df.duration.v1, V1ratio ~ q_c2 * condition * age, FUN = length)

ezANOVA (df.duration.v1, .(V1ratio), .(speaker), .(q_c2, condition), between = .(age))

# test C2 ratio
dotplot(data = df.duration.c2, C2ratio ~ q_v1 | condition * age)
ezANOVA (df.duration.c2, .(C2ratio), .(speaker), .(q_v1, condition), between = .(age))


#################################
## Absolute duration of V1 and C2
#

legend.category = df$category
levels(legend.category) = list(
  "VːC"  = "VːC",
  "VC"   = "VC",
  "VːCː" = "VːCː",
  "VCː"  = "VCː"
)

colors = c('black', 'blue', 'red', 'green')
xlim = c(0, 500)
ylim = c(0, 500)
xlab = "Duration of C2 [ms]"
ylab = "Duration of V1 [ms]"
key = list(space = "top",
           points = list(col = colors, pch = 16),
           text = list(levels(legend.category)),
           between = 1,
           columns = 2
)

#####
# Figure 3
# All words included
#
# Export this plot as SVG, size 400x480
xyplot (v1dur ~ c2dur_total | condition * age,
        df,
        groups = legend.category,
        key = key, col = colors, aspect = 1,
        xlim = xlim,  ylim = ylim,
        xlab = xlab,  ylab = ylab)


xyplot ( data = df, v1dur ~ c2dur_total | speaker * condition, groups = interaction(q_v1, q_c2), col=c('black', 'blue', 'red', 'green'))