---
title: "rMA thesis"
author: "ChengjiaYe"
date: "2023-05-30"
output: html_document
---

```{r}
yng <- read.delim("results_young.txt", stringsAsFactors=TRUE)
yng <- yng[, -c(12,13,19,20)]
yng
```

```{r}
levels(yng$vowel)
is.integer(yng$Duration)

result_duration <- data.frame(vowel = character(),
                        mean = numeric(),
                        sd = numeric(),
                        stringsAsFactors = FALSE)

for (vowel in levels(yng$vowel)){
  dt.duration <- yng[yng$vowel == vowel,"Duration"]
  vowel_mean <- round(mean(dt.duration))
  vowel_sd <- round(sd(dt.duration))
  
  new_row <- data.frame(vowel = vowel,
                        mean = vowel_mean,
                        sd = vowel_sd,
                        stringsAsFactors = FALSE)
  
  # Append the new row to the result data frame
  result_duration <- rbind(result_duration, new_row)
}
print(result_duration)
```
```{r}
for (vowel in levels(yng$vowel)){
  dt.F1midpoint <- yng[yng$vowel == vowel,"F1_midpoint"]
  vowel_mean <- round(mean(dt.duration))
  vowel_sd <- round(sd(dt.duration))
  
  new_row <- data.frame(vowel = vowel,
                        mean = vowel_mean,
                        sd = vowel_sd,
                        stringsAsFactors = FALSE)
  
  # Append the new row to the result data frame
  result_duration <- rbind(result_duration, new_row)
}
print(result_duration)
```

```{r}
result_formant <- data.frame(vowel = character(),
                        mean.female_f1_50 = numeric(), mean.female_f2_50 = numeric(),
                        mean.female_f1_75 = numeric(), mean.female_f2_75 = numeric(),
                        mean.male_f1_50 = numeric(), mean.male_f2_50 = numeric(),
                        mean.male_f1_75 = numeric(), mean.male_f2_75 = numeric(),
                        sd.female_f1_50 = numeric(), sd.female_f2_50 = numeric(),
                        sd.female_f1_75 = numeric(), sd.female_f2_75 = numeric(),
                        sd.male_f1_50 = numeric(), sd.male_f2_50 = numeric(),
                        sd.male_f1_75 = numeric(), sd.male_f2_75 = numeric(),
                        stringsAsFactors = FALSE)

for (vowel in levels(yng$vowel)){
  dt.male_f1_50 <- yng[yng$vowel == vowel & yng$gender == "male","F1_midpoint"]
  dt.male_f2_50 <- yng[yng$vowel == vowel & yng$gender == "male","F2_midpoint"]
  dt.male_f1_75 <- yng[yng$vowel == vowel & yng$gender == "male","F1_0.75point"]
  dt.male_f2_75 <- yng[yng$vowel == vowel & yng$gender == "male","F2_0.75point"]
  dt.female_f1_50 <- yng[yng$vowel == vowel & yng$gender == "female","F1_midpoint"]
  dt.female_f2_50 <- yng[yng$vowel == vowel & yng$gender == "female","F2_midpoint"]
  dt.female_f1_75 <- yng[yng$vowel == vowel & yng$gender == "female","F1_0.75point"]
  dt.female_f2_75 <- yng[yng$vowel == vowel & yng$gender == "female","F2_0.75point"]
  male_f1_50_mean <- round(mean(dt.male_f1_50))
  male_f1.50_sd <- round(sd(dt.male_f1_50))
  male_f2_50_mean <- round(mean(dt.male_f2_50))
  male_f2.50_sd <- round(sd(dt.male_f2_50))
  male_f1_75_mean <- round(mean(dt.male_f1_75))
  male_f1.75_sd <- round(sd(dt.male_f1_75))
  male_f2_75_mean <- round(mean(dt.male_f2_75))
  male_f2.75_sd <- round(sd(dt.male_f2_75))
  female_f1_50_mean <- round(mean(dt.female_f1_50))
  female_f1.50_sd <- round(sd(dt.female_f1_50))
  female_f2_50_mean <- round(mean(dt.female_f2_50))
  female_f2.50_sd <- round(sd(dt.female_f2_50))
  female_f1_75_mean <- round(mean(dt.female_f1_75))
  female_f1.75_sd <- round(sd(dt.female_f1_75))
  female_f2_75_mean <- round(mean(dt.female_f2_75))
  female_f2.75_sd <- round(sd(dt.female_f2_75))
  
  new_row <- data.frame(vowel = vowel,
                        mean.female_f1_50 = female_f1_50_mean, mean.female_f2_50 = female_f2_50_mean,
                        mean.female_f1_75 = female_f1_75_mean, mean.female_f2_75 = female_f2_75_mean,
                        mean.male_f1_50 = male_f1_50_mean, mean.male_f2_50 = male_f2_50_mean,
                        mean.male_f1_75 = male_f1_75_mean, mean.male_f2_75 = male_f2_75_mean,
                        sd.female_f1_50 = female_f1.50_sd, sd.female_f2_50 = female_f2.50_sd,
                        sd.female_f1_75 = female_f1.75_sd, sd.female_f2_75 = female_f2.75_sd,
                        sd.male_f1_50 = male_f1.50_sd, sd.male_f2_50 = male_f2.50_sd,
                        sd.male_f1_75 = male_f1.75_sd, sd.male_f2_75 = male_f2.75_sd,
                        stringsAsFactors = FALSE)
  
  # Append the new row to the result data frame
  result_formant <- rbind(result_formant, new_row)
}
print(result_formant)
```

```{r}
library(lmerTest)
contrast.syllable <- cbind(c(-0.5,+0.5))
colnames(contrast.syllable) <- c("-closed+open")
contrasts(yng$syllable) <- contrast.syllable
contrasts(yng$syllable)

contrast.gender <- cbind(c(+0.5,-0.5))
colnames(contrast.gender) <- c("+f-m")
contrasts(yng$gender) <- contrast.gender
contrasts(yng$gender)

contrast.register <- cbind(c(-0.5,+0.5))
colnames(contrast.register) <- c("-h+l")
contrasts(yng$register) <- contrast.register
contrasts(yng$register)

contrast.highness <- cbind(c(1/3,1/3,-2/3), c(-0.5,0.5,0)) #high,low,mid
colnames(contrast.highness) <- c("+HL-M","-H+L")
contrasts(yng$highness) <- contrast.highness
contrasts(yng$highness)

contrast.backness <- cbind(c(1/3,1/3,-2/3), c(-0.5,0.5,0)) #back,front,mid
colnames(contrast.backness) <- c("+BF-M","-B+F")
contrasts(yng$backness) <- contrast.backness
contrasts(yng$backness)

duration.mdl <- lmerTest::lmer(Duration ~ syllable * gender * register + (syllable * register|participant), data = yng)
summary(duration.mdl)

duration.mdl <- lmerTest::lmer(Duration ~ syllable * gender * register + (syllable + register|participant), data = yng)
summary(duration.mdl)
```

```{r}
library(emmeans)
summary(pairs(emmeans(duration.mdl, ~ syllable | gender | register), adjust = "bonferroni"))
summary(pairs(emmeans(duration.mdl, ~ register | syllable), adjust = "bonferroni"))
summary(pairs(emmeans(duration.mdl, ~ syllable | register), adjust = "bonferroni"))
```

```{r}
library(viridis)
library(ggplot2)
library(hrbrthemes)
library(extrafont)
loadfonts(device = "win")
```

```{r}
ggplot(yng, aes(x=factor(vowel, levels = c("i","I","y","Y","e","round_e", "E","schwa", "A", "u", "ou", "o","unround_o","low_o","a"), labels = c("i","ɪʔ","y","ʏʔ","e","ø", "ɛ","əʔ", "Aʔ", "u", "o", "oʔ","ɤ","ɔ","ɑ")), y=Duration, fill=vowel)) +
    geom_violin() +
    scale_fill_viridis(discrete = TRUE, alpha=0.6, option="A") +
    theme_ipsum() +
    theme(
      legend.position = "none",
      plot.title = element_text(family="serif", face="bold", size=12), 
      axis.title.x  = element_text(family = "serif"),
      axis.title.y  = element_text(family = "serif"),
      axis.text.x   = element_text(family = "serif"),
      axis.text.y   = element_text(family = "serif")
    ) +
    ggtitle("Violin chart of the distribution of vowel duration") +
    xlab("")
```

```{r}
ggplot(yng, aes(x=factor(vowel, levels = c("i","I","y","Y","e","round_e", "E","schwa", "A", "u", "ou", "o","unround_o","low_o","a"), labels = c("i","ɪʔ","y","ʏʔ","e","ø", "ɛ","əʔ", "Aʔ", "u", "o", "oʔ","ɤ","ɔ","ɑ")), y=Duration, fill=register)) +
    ylab("duration (ms)") +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.5, option="G") +
    theme_ipsum() +
    theme(
      legend.position = "bottom",
      #plot.title = element_text(family="Times", face="bold", size=12), 
      axis.title.x  = element_text(family = "serif"),
      axis.title.y  = element_text(family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 12),
      axis.text.y   = element_text(family = "serif", size = 12),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", face="bold", size = 12)
    ) +
    ggtitle("Boxplot of the distribution of vowel duration") +
    xlab("")
```


```{r}
ggplot(yng, aes(x=factor(vowel, levels = c("i","I","y","Y","e","round_e", "E","schwa", "A", "u", "ou", "o","unround_o","low_o","a"), labels = c("i","ɪʔ","y","ʏʔ","e","ø", "ɛ","əʔ", "Aʔ", "u", "o", "oʔ","ɤ","ɔ","ɑ")), y=Duration, fill=vowel)) +
    xlab("") + ylab("duration (ms)") + 
    ggtitle("Boxplot of the distribution of vowel duration") +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.7, option="D") +
    theme_ipsum() +
    theme(
      legend.position = "none",
      #plot.title = element_text(family="Times", face="bold", size=12), 
      axis.title.x  = element_text(family = "serif"),
      axis.title.y  = element_text(family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 12),
      axis.text.y   = element_text(family = "serif", size = 12)
    ) 
```

```{r}
library(dplyr)
means <- yng %>%
    group_by(vowel) %>%
    summarise(mean_F1 = mean(F1_midpoint),
              mean_F2 = mean(F2_midpoint))
  
ggplot(yng, aes(x = F2_midpoint, y = F1_midpoint
, color = vowel, label = factor(vowel, levels = c("i","I","y","Y","e","round_e", "E","schwa", "A", "u", "ou", "o","unround_o","low_o","a"), labels = c("i","ɪʔ","y","ʏʔ","e","ø", "ɛ","əʔ", "Aʔ", "u", "o", "oʔ","ɤ","ɔ","ɑ")))) + 
    geom_point(alpha = 0.1) + 
    stat_ellipse(level = 0.68) +
    geom_text(data = means, aes(x = mean_F2, y = mean_F1, family="serif", fontface="bold"), size=12) + 
    scale_x_reverse(limits = c(3000, 500)) + scale_y_reverse(limits = c(1200, 200)) + 
    scale_color_discrete() +
    guides(color = FALSE) + 
    theme_classic() +
    theme(
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))


means75 <- yng %>%
    group_by(vowel) %>%
    summarise(mean_F175 = mean(F1_0.75point),
              mean_F275 = mean(F2_0.75point))
ggplot(yng, aes(x = F2_0.75point, y = F1_0.75point
, color = vowel, label = factor(vowel, levels = c("i","I","y","Y","e","round_e", "E","schwa", "A", "u", "ou", "o","unround_o","low_o","a"), labels = c("i","ɪʔ","y","ʏʔ","e","ø", "ɛ","əʔ", "Aʔ", "u", "o", "oʔ","ɤ","ɔ","ɑ")))) + 
    geom_point(alpha = 0.1) + 
    stat_ellipse(level = 0.68) + 
    geom_text(data = means75, aes(x = mean_F275, y = mean_F175, family="serif", fontface="bold"), size=12) + 
    scale_x_reverse(limits = c(3000, 500)) + scale_y_reverse(limits = c(1200, 200)) + 
    scale_color_discrete() +
    guides(color = FALSE) + 
    theme_classic() +
    theme(
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))
```


```{r}
library(dplyr)
library(ggpubr)
means50 <- yng %>%
    group_by(vowel) %>%
    summarise(mean_F150 = mean(F1_midpoint),
              mean_F250 = mean(F2_midpoint))
means75 <- yng %>%
    group_by(vowel) %>%
    summarise(mean_F175 = mean(F1_0.75point),
              mean_F275 = mean(F2_0.75point))


open_mid <- ggplot(yng[yng$syllable == "open",], aes(x = F2_midpoint, y = F1_midpoint
, color = vowel, label = factor(vowel, levels = c("i","y","e","round_e", "E", "u", "ou", "unround_o","low_o","a"), labels = c("i","y","e","ø", "ɛ", "u", "o", "ɤ","ɔ","ɑ")))) + 
    #xlab("F2 at the midpoint (Hz)") +
    ylab("Open-ended vowels \n F1 (Hz)") +
    ggtitle("Vowel space at the midpoint") +
    geom_point(size = 0.95, alpha = 0.2) + 
    stat_ellipse(level = 0.68) + 
    geom_text(data = means50, aes(x = mean_F250, y = mean_F150, family="serif", 
                                #fontface="bold"
                                ), size = 10) + 
    scale_x_reverse(limits = c(3000, 500)) + scale_y_reverse(limits = c(1200, 200)) + 
    scale_color_discrete() +
    guides(color = FALSE) + 
    theme_classic() +
    theme(
    plot.title    = element_text(family = "serif", size = 14, hjust = 0.5),
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))

open_upper <- ggplot(yng[yng$syllable == "open",], aes(x = F2_0.75point, y = F1_0.75point
, color = vowel, label = factor(vowel, levels = c("i","y","e","round_e", "E", "u", "ou", "unround_o","low_o","a"), labels = c("i","y","e","ø", "ɛ", "u", "o", "ɤ","ɔ","ɑ")))) + 
    #xlab("F2 at the upper quartile (Hz)") +
    #ylab("Open-ended vowels \n F1 (Hz)") +
    ggtitle("Vowel space at the upper quartile") +
    geom_point(size = 0.95, alpha = 0.2) + 
    stat_ellipse(level = 0.68) + 
    geom_text(data = means75, aes(x = mean_F275, y = mean_F175,  family="serif", 
                                  #fontface="bold"
                                  ), size = 10) + 
    scale_x_reverse(limits = c(3000, 500)) + scale_y_reverse(limits = c(1200, 200)) + 
    scale_color_discrete() +
    guides(color = FALSE) + 
    theme_classic() +
    theme(
    plot.title    = element_text(family = "serif", size = 14, hjust = 0.5),
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))

closed_mid <- ggplot(yng[yng$syllable == "closed",], aes(x = F2_midpoint, y = F1_midpoint
, color = vowel, label = factor(vowel, levels = c("I","Y","schwa","A", "o"), labels = c("ɪ","ʏ","ə", "A", "o")))) + 
    xlab("F2 (Hz)") +
    ylab("Glottal-ended vowels \n F1 (Hz)") +
    #ggtitle("Vowel space at the midpoint of glottal-ended vowels") +
    geom_point(size = 0.95, alpha = 0.2) + 
    stat_ellipse(level = 0.68) + 
    geom_text(data = means50, aes(x = mean_F250, y = mean_F150, family="serif", 
                                #fontface="bold"
                                ), size = 10) + 
    scale_x_reverse(limits = c(3000, 500)) + scale_y_reverse(limits = c(1200, 200)) + 
    scale_color_discrete() +
    guides(color = FALSE) + 
    theme_classic() +
    theme(
    plot.title    = element_text(family = "serif", size = 14, hjust = 0.5),
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))

closed_upper <- ggplot(yng[yng$syllable == "closed",], aes(x = F2_0.75point, y = F1_0.75point
, color = vowel, label = factor(vowel, levels = c("I","Y","schwa","A", "o"), labels = c("ɪ","ʏ","ə", "A", "o")))) + 
    xlab("F2 (Hz)") +
    #ylab("F1 at the upper quartile (Hz)") +
    #ggtitle("Vowel space at the upper quartile of glottal-ended vowels") +
    geom_point(size = 0.95, alpha = 0.2) + 
    stat_ellipse(level = 0.68) + 
    geom_text(data = means75, aes(x = mean_F275, y = mean_F175,  family="serif", 
                                  #fontface="bold"
                                  ), size = 10) + 
    scale_x_reverse(limits = c(3000, 500)) + scale_y_reverse(limits = c(1200, 200)) + 
    scale_color_discrete() +
    guides(color = FALSE) + 
    theme_classic() +
    theme(
    plot.title    = element_text(family = "serif", size = 14, hjust = 0.5),
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))

ggarrange(open_mid + rremove("x.title"), open_upper + rremove("x.title") +rremove("y.title") , closed_mid, closed_upper+rremove("y.title") ,
          ncol = 2, nrow = 2)

```

```{r}
ggplot(yng[yng$syllable == "open",], aes(x = log(F2_midpoint), y = log(F1_midpoint)
, color = vowel, label = factor(vowel, levels = c("i","y","e","round_e", "E", "u", "ou", "unround_o","low_o","a"), labels = c("i","y","e","ø", "ɛ", "u", "o", "ɤ","ɔ","ɑ")))) + 
    xlab("F2 at the midpoint (Hz)") +
    ylab("F1 at the midpoint (Hz)") +
    geom_point(alpha = 0.1) + 
    stat_ellipse(level = 0.68) + 
    geom_text(data = means, aes(x = log(mean_F2), y = log(mean_F1), family="serif", fontface="bold"), size=12) + 
    scale_x_reverse(limits = c(log(3000), log(500))) + scale_y_reverse(limits = c(log(1200), log(200))) + 
    scale_color_discrete() +
    guides(color = FALSE) + 
    theme(
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))+
    theme_classic()
```

```{r}
F1mid.mdl <- lmerTest::lmer(F1_midpoint ~ syllable * gender * register * highness + (syllable + register|participant), data = yng)
summary(F1mid.mdl)
```

```{r}
summary(pairs(emmeans(F1mid.mdl, ~ register | highness | syllable), adjust = "bonferroni"))
```

```{r}
summary(pairs(emmeans(F1mid.mdl, ~ highness | gender | syllable), adjust = "bonferroni"))
```

```{r}
lowplot <- ggplot(yng[yng$highness == "low", ], aes(x=gender, labels = gender, y=F1_midpoint, fill=syllable)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.5, option="C") +
    theme_classic() + scale_y_reverse(limits = c(1200, 200)) +
    theme(
      legend.position = "bottom",
      axis.title.x  = element_text(family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 12),
      axis.text.y   = element_text(family = "serif", size = 12),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", face="bold", size = 12)
    ) +
    xlab("low vowels") +
    ylab("F1 at the midpoint (Hz)")

midplot <- ggplot(yng[yng$highness == "mid", ], aes(x=gender, labels = gender, y=F1_midpoint, fill=syllable)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.5, option="C") +
    theme_classic() + scale_y_reverse(limits = c(1200, 200)) +
    theme(
      legend.position = "bottom",
      #plot.title = element_text(family="Times", face="bold", size=12), 
      axis.title.x  = element_text(family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 12),
      axis.text.y   = element_text(family = "serif", size = 12),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", face="bold", size = 12)
    ) +
    xlab("mid vowels")

highplot <- ggplot(yng[yng$highness == "high", ], aes(x=gender, labels = gender, y=F1_midpoint, fill=syllable)) +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.5, option="C") +
    theme_classic() + scale_y_reverse(limits = c(1200, 200)) +
    theme(
      legend.position = "bottom",
      axis.title.x  = element_text(family = "serif", face="bold", size = 12),
      axis.title.y  = element_text(family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 12),
      axis.text.y   = element_text(family = "serif", size = 12),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", face="bold", size = 12)
    ) +
    xlab("high vowels")

ggarrange(lowplot, midplot + rremove("y.title") , highplot + rremove("y.title") ,
          ncol = 3, nrow = 1)
```



```{r}
F1upper.mdl <- lmerTest::lmer(F1_0.75point ~ syllable * gender * register * highness + (syllable |participant), data = yng)
summary(F1upper.mdl)
```
```{r}
summary(pairs(emmeans(F1upper.mdl, ~ syllable | register | highness), adjust = "bonferroni"))
summary(pairs(emmeans(F1upper.mdl, ~ register | highness | syllable), adjust = "bonferroni"))
```

```{r}
summary(pairs(emmeans(F1upper.mdl, ~ gender | register| highness ), adjust = "bonferroni"))
```

```{r}
F2mid.mdl <- lmerTest::lmer(F2_midpoint ~  gender * syllable * register * backness + (syllable|participant), data = yng)
summary(F2mid.mdl)
```
```{r}
summary(pairs(emmeans(F2mid.mdl, ~ backness | syllable ), adjust = "bonferroni"))

summary(pairs(emmeans(F2mid.mdl, ~ syllable | backness ), adjust = "bonferroni"))

```


```{r}
F2upper.mdl <- lmerTest::lmer(F2_0.75point ~  gender * syllable * register * backness + (syllable|participant), data = yng)
summary(F2upper.mdl)
```
```{r}
summary(pairs(emmeans(F2upper.mdl, ~ gender | backness | syllable ), adjust = "bonferroni"))
```

```{r}
# t-test for F1
for (vowel in levels(yng$vowel)){
  t_test <- t.test(yng[yng$vowel == vowel,"F1_midpoint"], yng[yng$vowel == vowel,"F1_0.75point"],paired = TRUE)
  print(vowel)
  print(t_test) 
}
```

```{r}
# t-test for F2
for (vowel in levels(yng$vowel)){
  t_test <- t.test(yng[yng$vowel == vowel,"F2_midpoint"], yng[yng$vowel == vowel,"F2_0.75point"],paired = TRUE)
  print(vowel)
  print(t_test) 
}
```

```{r}
breathy <- lmerTest::lmer(H1.H2 ~  gender * syllable * register + (syllable + register |participant), data = yng)
summary(breathy)
```

```{r}
creaky <- glmer(creakiness ~ gender * syllable * register + (register |participant) , data = yng ,family = "binomial")
summary(creaky)
summary(pairs(emmeans(creaky, ~ register | syllable ), adjust = "bonferroni"))
```

```{r}
F1a <- yng[,c("vowel","F1_midpoint")]
colnames(F1a) <- c("vowel","F1")
F1b <- yng[,c("vowel","F1_0.75point")]
colnames(F1b) <- c("vowel","F1")
overalF1 <- rbind(F1a,F1b)
F2a <- yng[,c("vowel","F2_midpoint","syllable")]
colnames(F2a) <- c("vowel","F2","syllable")
F2b <- yng[,c("vowel","F2_0.75point","syllable")]
colnames(F2b) <- c("vowel","F2","syllable")
overalF2 <- rbind(F2a,F2b)
overalF2 <- overalF2[, -1]
#merge two data frames by vowel
overal <- cbind(overalF1, overalF2)
```

```{r}
means <- overal %>%
    group_by(vowel) %>%
    summarise(mean_F1 = mean(F1),
              mean_F2 = mean(F2))
ggplot(overal, aes(x = F2, y = F1, color = vowel, label = factor(vowel, levels = c("i","I","y","Y","e","round_e", "E","schwa", "A", "u", "ou", "o","unround_o","low_o","a"), labels = c("i","ɪʔ","y","ʏʔ","e","ø", "ɛ","əʔ", "ɐʔ", "ɯ", "u", "oʔ","ɤ","ɔ","ɑ")))) + 
    xlab("F2 (Hz)") +
    ylab("F1 (Hz)") +
    #ggtitle("Vowel space at the midpoint of glottal-ended vowels") +
    geom_point(size = 0.95, alpha = 0.15) + 
    stat_ellipse(level = 0.50, na.rm = TRUE, lwd = 1.0) + 
    geom_text(data = means, aes(x = mean_F2, y = mean_F1, family="serif"), size=9) +
    scale_x_reverse(limits = c(3000, 500)) + scale_y_reverse(limits = c(1200, 200)) +
    
    scale_color_manual(values=c("#21918c", "#443983", "#31688e", "#90d743","#fde725","#443983","#443983","#90d743", "#440154","#440154", "#fde725","#35b779","#31688e","#35b779", "#21918c"
                                 #"turquoise4","goldenrod","gold2","turquoise1","dodgerblue2", "darkorange","deepskyblue","royalblue","royalblue4"
                                 # "#fde725", "#31688e", "#90d743", "#443983", "#35b779", "#440154", "#21918c"
                                 )) +
    guides(color = FALSE) + 
    theme_classic() +
    theme(
    axis.title.x  = element_text(family = "serif", size = 12),
    axis.title.y  = element_text(family = "serif", size = 12),
    axis.text.x   = element_text(family = "serif", size = 12),
    axis.text.y   = element_text(family = "serif", size = 12),
    legend.title  = element_text(family = "serif", face="bold", size = 12))
```
```{r}
ggplot(yng, aes(x=syllable, y=H1.H2, fill=register)) +
    ylab("H1-H2 difference (dB/Hz)") +
    geom_boxplot() +
    scale_fill_viridis(discrete = TRUE, alpha=0.5, option="G") +
    theme_ipsum() +
    theme(
      legend.position = "bottom",
      #plot.title = element_text(family="Times", face="bold", size=12), 
      axis.title.x  = element_text(family = "serif"),
      axis.title.y  = element_text(family = "serif", face="bold", size = 12),
      axis.text.x   = element_text(family = "serif", size = 12),
      axis.text.y   = element_text(family = "serif", size = 12),
      legend.title  = element_text(family = "serif", face="bold", size = 12),
      legend.text = element_text(family = "serif", face="bold", size = 12)
    ) +
    # ggtitle("Boxplot of the distribution of vowel duration") +
    xlab("")
```





```{r}
median(yng$spectrum_difference, na.rm = T)
quantile(yng$spectrum_difference, na.rm = T)
upper.bound <- 0.05499226 + 1.5* (0.05499226 - 0.01196686)
lower.bound <- 0.01196686 - 1.5* (0.05499226 - 0.01196686)
quantile(yng$spectrum_difference, 0.95, na.rm = T)

ggplot(yng, aes(x=spectrum_difference, na.rm = T)) + 
    geom_histogram(aes(y=..density..),      # Histogram with density instead of count on y-axis
                   binwidth=.01,
                   colour="black", fill="white") +
    geom_density(alpha=.2, fill="#FF6666")  # Overlay with transparent density plot


```

```{r}
yng <- yng[yng$spectrum_difference < 0.15, ]
breathy2 <- lmerTest::lmer(H1.H2 ~  gender * syllable * register + (syllable + register |participant), data = yng)
summary(breathy2)
```




