# R version: 4.2.2

library(tidyverse) # version 1.3.1
library(ggpubr) # version 0.4.0
library(lmerTest) # version 3.1-3
library(emmeans) # version 1.8.2

# 3.1 Postvocalic stops in word-medial position ---------------------------

medial = read.csv("medial.csv")


# sort groups
medial$Group = factor(medial$Group, c("Standard",
                                      "younger Bavarian",
                                      "older Bavarian"))

# Figure 1 - closure duration 
ggplot(medial) +
  aes(x = Group, y = pClo, fill = Phonotactics) +
  geom_boxplot() +
  facet_grid(Vowel~Stop) +
  scale_fill_manual(values=c("grey", "white")) +
  labs(fill = "Phonotactics") +
  ylab(bquote(pClosure)) +
  xlab("Speaker group") +
  theme(axis.text.x = element_text(size = 12, angle=15, hjust=1),
        axis.text.y = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14),
        legend.position = "top")

# LMM - closure duration
pClo_med.lmer = lmer(pClo ~ Group * Vowel * Stop + 
                        (1 + Vowel * Stop | session) +
                        (Group | word), 
                      data = medial)
pClo_med.step = step(pClo_med.lmer)
anova(get_model(pClo_med.step))

# pairwise comparisons
emmeans(pClo_med.lmer,
        pairwise ~ Group * Stop | Vowel,
        lmer.df = "satterthwaite")

# Figure 2 - VOT 
ggplot(medial) +
  aes(x = Group, y = pVOT, fill = Phonotactics) +
  geom_boxplot() +
  facet_grid(Vowel~Stop) +
  scale_fill_manual(values=c("grey", "white")) +
  labs(fill = "Phonotactics") +
  ylab(bquote(pVOT)) +
  xlab("Speaker group") +
  theme(axis.text.x = element_text(size = 12, angle=15, hjust=1),
        axis.text.y = element_text(size = 12),
        axis.title = element_text(size = 14),
        strip.text = element_text(size = 14),
        legend.text = element_text(size = 12),
        legend.title = element_text(size = 14),
        legend.position = "top")

# LMM VOT
pVOT_med.lmer = lmer(pVOT ~ Group * Vowel * Stop + 
                        (1 + Vowel * Stop | session)+
                        (Group | word), 
                      data = medial)
pVOT_med.step = step(pVOT_med.lmer)
anova(get_model(pVOT_med.step))

# pairwise comparisons
emmeans(pVOT_med.lmer,
        pairwise ~ Group * Stop | Vowel,
        lmer.df = "satterthwaite")


# 3.2 Prevocalic stops in word-initial position ---------------------------

initial = read.csv("initial.csv")

# sort groups
initial$Group = factor(initial$Group, c("Standard",
                                        "younger Bavarian",
                                        "older Bavarian"))

# Figure 3 

# Closure duration 
Fig_3.1 = ggplot(initial) +
  aes(x = Group, y = pClo) +
  geom_boxplot()  + facet_grid(~ Stop) +
  ylab(bquote(pClosure)) +
  xlab("Speaker group") +
  theme(axis.text.x = element_text(size = 10, angle=15, hjust=1))

# VOT
Fig_3.2 = ggplot(initial) +
  aes(x = Group, y = pVOT) +
  geom_boxplot()  + facet_grid(~ Stop) +
  ylab(bquote(pVOT)) +
  xlab("Speaker group") +
  theme(axis.text.x = element_text(size = 10, angle=15, hjust=1))

# combine plots and prepare final Figure
Fig_3.3= ggarrange(Fig_3.1 + rremove("xlab"),
                   Fig_3.2 + rremove("xlab"),
                   ncol = 1,nrow = 2,
                   common.legend = TRUE)
annotate_figure(Fig_3.3,
                bottom = text_grob("Speaker group",
                                   size = 12, face= "plain")
                )


# LMM closure duration
pClo_ini.lmer = lmer(pClo ~ Group * Stop + 
                      (1 + Stop | session) +
                      (1 + Group | word), 
                    data = initial)
step(pClo_ini.lmer)
anova(get_model(step(pClo_ini.lmer)))

# pairwise comparison
emmeans(pClo_ini.lmer,
        pairwise ~ Group | Stop,
        lmer.df = "satterthwaite")

# LMM VOT 
pVOT_ini.lmer = lmer(pVOT ~ Group * Stop + 
                      (1 + Stop | session) +
                      (1 + Group | word), 
                    data = initial)
step(pVOT_ini.lmer)
anova(get_model(step(pVOT_ini.lmer)))

# pairwise comparison
emmeans(pVOT_ini.lmer,
        pairwise ~ Group | Stop,
        lmer.df = "satterthwaite")

# 4. Discussion and Conclusion: Correlation  ------------------------------
# correaltion between closure duration and VOT in tokens with underlying fortis stops
temp1 = medial$Stop== "fortis"
fortis = medial[temp1,]

# Figure (not shown in paper)
ggplot(fortis, aes(pClo, pVOT)) +
  geom_point(aes(color = Group, shape = Group)) +
  scale_shape_manual(values = c(1, 2, 3)) +
  geom_smooth(aes(linetype = Group, color = Group), method = "lm", se = F) 

# Pearson correlation coefficients per speaker group
fortis %>% 
  group_by(Group) %>% 
  summarise(corr = cor(pVOT, pClo)) %>%
  ungroup()