#' ---
#' title: "Interaction between Phrasal Structure and Vowel Tenseness in German: An Acoustic and Articulatory Study"
#' author: "Malte Belz, Oksana Rasskazova, Jelena Krivokapic, Christine Mooshammer"
#' date: 2021
#' output:
#'    html_document:
#'      toc: true
#'      toc_depth: 4
#'      toc_float: true
#'      number_sections: true
#' ---

options(OutDec= ".")
set.seed(1)

#' Packages needed:
library(ggplot2)
library(lme4)
library(lmerTest)
library(effects)
library(car)
library(xtable)
library(dplyr)
library(tidyr)
library(texreg)
library(MuMIn)
library(ggpubr)
library(knitr)
library(emmeans)

#' # Acoustic analysis
#### Acoustic analysis ####
# read monosyllabic data
mono <- read.table("mono.txt",sep="\t")
mono$tenseness<-factor(mono$tenseness, levels = c("tense","lax"))
mono$position<-factor(mono$position, levels = c("phrase-medial","phrase-final"))
levels(mono$iteration)<-c("1","1","2","2","3","3","4","4","5","5")


# read disyllabic data
di <- read.table("di.txt",sep="\t")
di$tenseness<-factor(di$tenseness, levels = c("tense","lax"))
di$position<-factor(di$position, levels=c("phrase-medial","phrase-final"))
levels(di$iteration)<-c("1","1","2","2","3","3","4","4","5")

# mono and di contain
# columns 1:16 from emuR
# column 17: subject id
# column 18: phrasal position
# column 19: target word
# column 20: iteration: number of repetition
# column 21: pausedur: duration of the following pause
# column 22: labelsdur: duration of the phonetic segment (not to use: to gather one data point per row, the actual duration of
#             e.g., aspiration of /b/, is not included in this column. It can be deducted implicitly by subtracting sylposdur with sylseglab == "bh"
#             from labelsdur with labelsfac == "b")
# column 23: labelsfac: phonological phonetic segment (factor coded), comprising segments with and without aspiration
# column 24: sylpos: syllable position (0: onset, N: nucleus, C: coda, O1: first syllable, O2: second syl., N1/N2: first/second nucleus)
# column 25: tenseness: tense or lax
# column 26: sylposdur: duration of syllable constituent (sylpos). It is important to note that this is the sum of both
#             the stop closure and stop aspiration (corr. to labelsfac)
# column 27: sylseglabel: actual phonetic label of the syllable position (combinations): [Sth, bh, th]


#' ## Data used
#' Unequal data points for respective target words are due to annotation decisions
#'
#' **Monosyllabic**
table(mono[mono$position=="phrase-medial",]$word,mono[mono$position=="phrase-medial",]$sylpos)
table(mono[mono$position=="phrase-final",]$word,mono[mono$position=="phrase-final",]$sylpos)
addmargins(table(mono$position,mono$sylpos))
addmargins(table(mono$word,mono$sylpos))
m<-addmargins(table(mono$word,mono$sylpos))


#' 784 target words were analyzed in total
#'
#' **Disyllabic**
table(di[di$position=="phrase-medial",]$word,di[di$position=="phrase-medial",]$sylpos)
table(di[di$position=="phrase-final",]$word,di[di$position=="phrase-final",]$sylpos)
#' 197 target words were analyzed
addmargins(table(di$position,di$sylpos))
addmargins(table(di$word,di$sylpos))
addmargins(table(di$word,di$position))
d<-addmargins(table(di$word,di$sylpos))

#' **Total data excluded**
# Monosyllabic: 2 position X 8 targets X 5 iterations X 10 subjects X 3 syllable positions
sprintf("%1.2f%%",round((1 - m[9,4]/(2 * 8 * 5 * 10 * 3))*100,3))
# Disyllabic: 2 position X 4 targets X 5 iterations X 10 subjects X 4 syllable positions
sprintf("%1.2f%%",round((1 - d[5,5]/(2 * 4 * 5 * 10 * 4))*100,3))
# Total total:
sprintf("%1.2f%%",round((1 - (m[9,4]+d[5,5])/((2 * 8 * 5 * 10 * 3) + (2 * 4 * 5 * 10 * 4)))*100,3))


#' ## Analysis of monosyllabic target words
### Analysis ####

#' Normalisation to speech rate: sylposdur (ms) divided by 1000 (sec)., multiplied with (sil/s), results in syllables
meta<-read.table("meta.speech.dur.txt",sep="\t", fileEncoding = "UTF-8")
meta$subject

# merge data with speech rate
mono<-merge(mono, meta[,c("subject","sprechgeschw")], by = "subject")
di<-merge(di, meta[,c("subject","sprechgeschw")], by = "subject")

# normalisation
mono<-mono %>% group_by(subject) %>%
  mutate(sylposdurnorm = sprechgeschw*(sylposdur/1000))
di<-di %>% group_by(subject) %>%
  mutate(sylposdurnorm = sprechgeschw*(sylposdur/1000))

head(as.data.frame(mono))
mono$sylpos<-factor(mono$sylpos, levels=c("O","N","C"))

# Monosyllabic data without Stiel/still
nostiel<-mono %>% filter(word %in% c("Bahn","Bann","Beet","Bett","Ruhm","Rum")) %>% droplevels()

# Overall model monosyllabic acoustic ####
m1<-lmer(sylposdur ~ position + sylpos + tenseness + (1|subject) + (1|word), data = mono)
plot(allEffects(m1))
qqPlot(resid(m1))
m2<-lmer(sylposdur ~ position * sylpos + tenseness + (1|subject) + (1|word), data = mono)
anova(m1,m2)
m3<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word), data = mono)
anova(m2,m3)
m4<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position|subject) + (0+position|word), data = mono)
anova(m3,m4,refit=F)
plot(allEffects(m4))
summary(m4)
#models with additional random slopes do not converge

qqPlot(resid(m4))
abline(h = cut<- 105)
mono2<-droplevels(subset(mono, resid(m4)<cut))
m4.1<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position|subject) + (0+position|word), data = mono2)
# after cut model m4.1 does not converge any more

mono2$position.num <- as.numeric(mono2$position)
mono2$tenseness.num <- as.numeric(mono2$tenseness)
mono2$sylpos.num <- as.numeric(mono2$sylpos)


m5<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position+tenseness|subject) + (0+position|word), data = mono2)
m5.r2dummy<-lmer(sylposdur ~ position.num * sylpos.num * tenseness.num +
           (1+position.num+tenseness.num|subject) + (0+position.num|word), data = mono2)
r.squaredGLMM(m5.r2dummy)
qqPlot(resid(m5))
summary(m5)
# best model is m5


# model withouth stiel
no5<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+position+tenseness|subject) + (0+position|word), data = nostiel)
plot(allEffects(no5))
plot(allEffects(m5))
summary(no5)


#' ## Overall lengthening for each syllable constituent ####

# Without tenseness group
mono %>% group_by(sylpos, position) %>%
  summarise(mean.ms = mean(sylposdur),sd.ms=sd(sylposdur))

mono %>% group_by(sylpos, position) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))

# With tenseness group
mono %>% group_by(sylpos, position,tenseness) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))

# normalized
mono %>% group_by(sylpos, position) %>%
  summarise(mean.norm = mean(sylposdurnorm)) %>%
  spread(position, mean.norm) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))

# Disyllabic without tenseness
di %>%  mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position, sylpos) %>%
  summarise(mean.ms = mean(sylposdur),sd.ms=sd(sylposdur))

di %>%  mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position, sylpos) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))

# Disyllabic with tenseness
di %>%  mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position, sylpos,tenseness) %>%
  summarise(mean.ms = mean(sylposdur)) %>%
  spread(position, mean.ms) %>%
  mutate(percent = `phrase-final` / `phrase-medial`) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))


#' ### Excursus: Normalization
#' Normalization of dependent variable for speech rate does not change the effects considerably
m5.norm<-lmer(sylposdurnorm ~ position * sylpos * tenseness + (1|subject) + (1|word) +
                (0+position+sylpos|subject) + (0+position|word), data = mono)
plot(allEffects(m5.norm))
plot(allEffects(m5))

#' ### Figure 2 (normalized stacked mean duration per tenseness and condition)
leng.perc<-mono %>% group_by(position,word, sylpos) %>%
  summarise(mean.dur = mean(sylposdurnorm)) %>%
  spread(position, mean.dur) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))
leng.perc.di<-di %>%  mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position,word, sylpos) %>%
  summarise(mean.dur = mean(sylposdurnorm)) %>%
  spread(position, mean.dur) %>%
  rename(medial = `phrase-medial` , final = `phrase-final`) %>%
  mutate(lengthening = paste(round((final/medial-1)*100),"%"))

### normalized duration in syllables (not in paper)
(plot.ac.mono<- mono %>% group_by(position,word,tenseness, sylpos) %>%
  dplyr::summarise(mean.dur = mean(sylposdurnorm), sd.dur = sd(sylposdurnorm)) %>%
  filter(position=="phrase-final") %>%
  ungroup() %>%
  mutate(lengthening = leng.perc$lengthening) %>%
  group_by(word) %>%
  mutate(ymin=cumsum(mean.dur) -sd.dur/2, ymax=cumsum(mean.dur) + sd.dur/2) %>%
  ggplot(aes(x=word, y=mean.dur, fill=sylpos))+
  geom_bar(
    stat="identity",
            position=position_stack(reverse=TRUE),alpha=0.9
           )+
  geom_errorbar(aes(ymin=ymin, ymax=ymax),
                stat="identity",
                width=.2                    # Width of the error bars
                )+
  coord_flip()+
  ylab("")+
  geom_text(aes(x=word, y=mean.dur, label = lengthening),
            color="black", hjust= 2,size = 3, position=position_stack(reverse=TRUE))+
  annotate("text", x=1:8,y= 2.7,label=rep(c("tense","lax"),4),size=3)+
  # ylim(0,2.7)+
  # scale_fill_grey(start=0,end=0.6, name="")+
  xlab("Monosyllabic targets")+
  scale_fill_brewer(palette="Set3",name="",labels=c("Onset","Nucleus","Coda"))+
  theme_bw(base_size = 10) +
  theme( legend.position="top"))
(plot.ac.di<-
  di %>% mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position,word,tenseness, sylpos) %>%
  summarise(mean.dur = mean(sylposdurnorm), sd.dur = sd(sylposdurnorm)) %>%
  filter(position=="phrase-final") %>%
  ungroup() %>%
  mutate(lengthening = leng.perc.di$lengthening) %>%
  group_by(word) %>%
  mutate(ymin=cumsum(mean.dur) -sd.dur/2, ymax=cumsum(mean.dur) + sd.dur/2) %>%
  ggplot(aes(x=word, y=mean.dur, fill=sylpos))+
  geom_bar(
    stat="identity",
    position=position_stack(reverse=T),alpha=0.9
  )+
  geom_errorbar(aes(ymin=ymin, ymax=ymax),
                stat="identity",
                width=.2                    # Width of the error bars
  )+
  coord_flip()+
  ylab("Normalized duration in syllables")+
  geom_text(aes(x=word, y=mean.dur, label = lengthening),
            color="black", vjust= 0,hjust=2,size = 3, position=position_stack(reverse=T))+
  annotate("text", x=1:4,y= 2.7,label=rep(c("tense","lax"),2),size=3)+
  # scale_fill_grey(start=0,end=0.6, name="")+
  xlab("Disyllabic targets")+
  scale_fill_brewer(palette="Set2",name="",labels=c("Onset 1","Nucleus 1","Onset 2","Nucleus 2"))+
  theme_bw(base_size = 10) +
  theme( legend.position="top"))

#' Figure 3: Acoustic mean duration of syllable constituents per phrasal position
# Duration in milliseconds
(plot.ac.mono<- mono %>% group_by(position,word,tenseness, sylpos) %>%
  dplyr::summarise(mean.dur = mean(sylposdur), sd.dur = sd(sylposdur)) %>%
  filter(position=="phrase-final") %>%
  ungroup() %>%
  mutate(lengthening = leng.perc$lengthening) %>%
  group_by(word) %>%
  mutate(ymin=cumsum(mean.dur) -sd.dur/2, ymax=cumsum(mean.dur) + sd.dur/2) %>%
  ggplot(aes(x=word, y=mean.dur, fill=sylpos))+
  geom_bar(
    stat="identity",
    position=position_stack(reverse=TRUE),alpha=0.9
  )+
  geom_errorbar(aes(ymin=ymin, ymax=ymax),
                stat="identity",
                width=.2                    # Width of the error bars
  )+
  coord_flip()+
  ylab("")+
  geom_text(aes(x=word, y=mean.dur, label = lengthening),
            color="black", hjust= 2,size = 3, position=position_stack(reverse=TRUE))+
  annotate("text", x=1:8,y= 500,label=rep(c("tense","lax"),4),size=3)+
  # ylim(0,2.7)+
  # scale_fill_grey(start=0,end=0.6, name="")+
  xlab("Monosyllabic targets")+
  scale_fill_brewer(palette="Set3",name="",labels=c("Onset","Nucleus","Coda"))+
  theme_bw(base_size = 10) +
  theme( legend.position="top"))
(plot.ac.di<-
  di %>% mutate(sylpos=factor(sylpos, levels = c("O1","N1","O2","N2"))) %>%
  group_by(position,word,tenseness, sylpos) %>%
  summarise(mean.dur = mean(sylposdur), sd.dur = sd(sylposdur)) %>%
  filter(position=="phrase-final") %>%
  ungroup() %>%
  mutate(lengthening = leng.perc.di$lengthening) %>%
  group_by(word) %>%
  mutate(ymin=cumsum(mean.dur) -sd.dur/2, ymax=cumsum(mean.dur) + sd.dur/2) %>%
  ggplot(aes(x=word, y=mean.dur, fill=sylpos))+
  geom_bar(
    stat="identity",
    position=position_stack(reverse=T),alpha=0.9
  )+
  geom_errorbar(aes(ymin=ymin, ymax=ymax),
                stat="identity",
                width=.2                    # Width of the error bars
  )+
  coord_flip()+
  ylab("Acoustic duration (ms)")+
  geom_text(aes(x=word, y=mean.dur, label = lengthening),
            color="black", vjust= 0,hjust=2,size = 3, position=position_stack(reverse=T))+
  annotate("text", x=1:4,y= 500,label=rep(c("tense","lax"),2),size=3)+
  # scale_fill_grey(start=0,end=0.6, name="")+
  xlab("Disyllabic targets")+
  scale_fill_brewer(palette="Set2",name="",labels=c("Onset 1","Nucleus 1","Onset 2","Nucleus 2"))+
  theme_bw(base_size = 10) +
  theme( legend.position="top"))
ggpubr::ggarrange(plot.ac.mono, plot.ac.di, common.legend = F,labels = c("a)","b)"),
          font.label = list(size=11,face="plain"),align = "h", widths = c(1.1,1), ncol = 1)
ggsave(file="v2.stackedacousticfinal.pdf", width = 200, height = 160, units = "mm")


# Plots of acoustic vowel length####
mono$wordpair <- mono$word
levels(mono$wordpair) <- c(rep("[a:] vs. [a]",2),
                           rep("[e:] vs. [ɛ]",2),
                           rep("[u:] vs. [ʊ]",2),
                           rep("[i:] vs. [ɪ]",2))
di$wordpair <- di$word
levels(di$wordpair) <- c(rep("[y:] vs. [ʏ]",2),rep("[i:] vs. [ɪ]",2))

(plot.ac.mono.data<-
  ggplot(mono %>% filter(sylpos=="N"), aes(x=tenseness, y=sylposdur,fill=position))+
  geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
             alpha=1,size=.1,col="black")+
  geom_boxplot(alpha=.3, notch = T, outlier.shape = NA, lwd = .3)+
  facet_wrap(~wordpair,nrow=1)+
  scale_fill_brewer(palette = "Set1")+
  xlab("Tenseness")+
  ylab("Acoustic duration (ms)")+
  theme_bw(base_size=10)+
  theme(legend.position = "top",legend.title = element_blank())
)
(plot.ac.di.data<-ggplot(di %>% filter(sylpos=="N1"), aes(x=tenseness, y=sylposdur, fill=position))+
  geom_boxplot(alpha=.3, notch=T, outlier.shape = NA, lwd = .3)+
  geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
             alpha=1,size=.1,col="black")+
  facet_wrap(~wordpair,nrow=1)+
  scale_fill_brewer(palette = "Set1")+
  # scale_fill_grey()+
  xlab("Tenseness")+
  ylab("")+
  theme_bw(base_size=10)+
  theme(legend.position = "top",legend.title = element_blank())
)
cairo_pdf(filename = "ac-vowel-tenseness-position.pdf",width=8,height=5)
ggpubr::ggarrange(plot.ac.mono.data, plot.ac.di.data, common.legend = T,labels = c("a)","b)"),
                  font.label = list(size=11,face="plain"),align = "h", widths = c(1.6,1))
dev.off()



#' ### Subset model acoustic monosyllabic for lax only
mono.lax<-droplevels(subset(mono, tenseness=="lax"))
mono.lax.nucl<-droplevels(subset(mono, tenseness=="lax" & sylpos=="N"))
mono.tense<-droplevels(subset(mono, tenseness=="tense"))
mono.tense.nucl<-droplevels(subset(mono, tenseness=="tense" & sylpos=="N"))

l1<-lmer(sylposdur ~ position + sylpos  + (1|subject) + (1|word), data = mono.lax)
l2<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word), data = mono.lax)
anova(l1,l2)
l3<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax)
anova(l2,l3,refit=F) #better
qqPlot(resid(l3))
abline(h = cut <- 50)
mono.lax.sub<-droplevels(subset(mono.lax, resid(l3) < cut))
l4<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
qqPlot(resid(l4))
abline(h = cut <- -50)
mono.lax.sub<-droplevels(subset(mono.lax.sub, resid(l4) > cut))
l5<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
qqPlot(resid(l5))
plot(allEffects(l5))
summary(l5)
lattice::dotplot(ranef(l5))


#' ### Proportional lengthening of vowels in nuclei, monosyllabic, with and without Stiel/still
tn1<-lmer(sylposdur ~ position  + (1|subject) + (0+position|word), data = mono.tense.nucl)
qqPlot(resid(tn1))
summary(tn1)
tn1.coef<-summary(tn1)$coefficients
# Proportion lengthening of tense vowel with Stiel
(tn1.coef[1,1]+tn1.coef[2,1])/tn1.coef[1,1] # tense with Stiel

# Model of proportional lengthening in monosyllabic nuclei
mono.proplen <- mono %>%
  filter(sylpos=="N") %>%
  group_by(subject,tenseness,wordpair,position) %>%
  summarise(mean = mean(sylposdur)) %>%
  tidyr::spread(position, mean) %>%
  mutate(proplen = ((`phrase-final`/`phrase-medial`)-1)*100)

mp1 <- lmer(proplen ~ tenseness + (1|subject) + (1|wordpair), data=mono.proplen)
qqPlot(resid(mp1))
summary(mp1)
plot(allEffects(mp1))

# lax only
ln1<-lmer(sylposdur ~ position  + (1|subject) + (1|word), data = mono.lax.nucl)
qqPlot(resid(ln1))
abline(h = cut <- c(-40,40))
mono.lax.nucl.sub<-droplevels(subset(mono.lax.nucl, resid(ln1) > cut[1] & resid(ln1) < cut[2]))
ln2<-lmer(sylposdur ~ position  + (1|subject) + (1|word), data = mono.lax.nucl.sub)
qqPlot(resid(ln2))
summary(ln2)
ln2.coef<-summary(ln2)$coefficients
# Proportion lengthening of lax vowels with Stiel
(ln2.coef[1,1]+ln2.coef[2,1])/ln2.coef[1,1]


#without Stiel/stiel tense
mono.tense.nucl.nostiel<-droplevels(subset(mono, tenseness=="tense" & sylpos=="N" & !(word%in%c("Stiel","still")) ))
mono.lax.nucl.nostiel<-droplevels(subset(mono, tenseness=="lax" & sylpos=="N" & !(word%in%c("Stiel","still")) ))

#tense
tn.nostiel<-lmer(sylposdur ~ position  + (1|subject) + (0+position|word), data = mono.tense.nucl.nostiel)
qqPlot(resid(tn.nostiel))
summary(tn.nostiel)
tn1.coef.nostiel<-summary(tn.nostiel)$coefficients
# Proportion lengthening of tense vowel without Stiel
(tn1.coef.nostiel[1,1]+tn1.coef.nostiel[2,1])/tn1.coef.nostiel[1,1]

#lax
ln.nostiel<-lmer(sylposdur ~ position  + (1|subject) + (0+position|word), data = mono.lax.nucl.nostiel)
qqPlot(resid(ln.nostiel))
summary(ln.nostiel)
ln1.coef.nostiel<-summary(ln.nostiel)$coefficients
# Proportion lengthening of lax vowel without Stiel
(ln1.coef.nostiel[1,1]+ln1.coef.nostiel[2,1])/ln1.coef.nostiel[1,1]



l2<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word), data = mono.lax)
anova(l1,l2)
l3<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax)
anova(l2,l3,refit=F) #better
qqPlot(resid(l3))
abline(h = cut <- 50)
mono.lax.sub<-droplevels(subset(mono.lax, resid(l3) < cut))
l4<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
qqPlot(resid(l4))
abline(h = cut <- -50)
mono.lax.sub<-droplevels(subset(mono.lax.sub, resid(l4) > cut))
l5<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = mono.lax.sub)
qqPlot(resid(l5))
plot(allEffects(l5))
summary(l5)
lattice::dotplot(ranef(l5))


# Model of proportional lengthening in disyllabic nuclei
di.proplen <- di %>%
  filter(sylpos=="N1") %>%
  group_by(subject,tenseness,wordpair,position) %>%
  summarise(mean = mean(sylposdur)) %>%
  tidyr::spread(position, mean) %>%
  mutate(proplen = ((`phrase-final`/`phrase-medial`)-1)*100)

dp1 <- lmer(proplen ~ tenseness + (1|subject) + (1|wordpair), data=di.proplen)
qqPlot(resid(dp1))
summary(dp1)
plot(allEffects(dp1))

# Table with proplen
x<-capture.output(
  texreg::texreg(list(texreg::extract(mp1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(dp1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                      # texreg::extract(mf3.1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
  ),
  digits = 1,booktabs = T,custom.model.names = c("Monosyllabic \\%","Disyllabic \\%"),
  float.pos = "!htbp", label = "tab:proplen", use.packages = F,single.row = T,caption.above = T,
  caption="Model predictions for the proportional lengthening of vowels  in  mono- and disyllabic words
   (standards error in parentheses) for the fixed effects tenseness (tense vs. lax) and random intercepts for subjects and word pairs.",
  dcolumn = T, fontsize = "scriptsize",center=T)
)
x[5]<-paste0(x[5], "\\centering")
x[16]<-paste0(x[16],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(mp1)[1],2),"/",round(r.squaredGLMM(mp1)[2],2),
              "&", round(r.squaredGLMM(dp1)[1],2), "/", round(r.squaredGLMM(dp1)[2],2),
              "\\\\ \n")
cat(x, sep = "\n")


#' ### Post-hoc tests (multiple comparisons, Tukey-corrected) for mono, lax
l5em <- emmeans(l5, pairwise~position*sylpos)
l5em.tab<-l5em$contrasts %>%
  summary(infer = T)
l5em.tab$p.value<-lapply(l5em.tab$p.value, FUN = function (x) {
  # levels of p-values
  if (x>=0.05)  p <- "n.\\,s." # not significant
  if (x<0.05)  p <- "<0,05"  # significant
  if (x<=0.01)  p <- "<0,01"  # significant
  if (x<=0.001) p <- "<0,001" # significant
  return(p)
})
l5em.tab<-as.data.frame(l5em.tab)
colnames(l5em.tab)<-c("Comparison","$\\beta$","s","df","Lower CI","Upper CI","t","p")
l5em.tab.sig<-subset(l5em.tab, p!="n.\\,s.")
l5em.tab.sig


#' ### Monosyllabic Effect plots
#' Preparation
ef1b<-as.data.frame(effect("position*sylpos*tenseness",m5))
ef1b$sylpos<-factor(ef1b$sylpos, levels=c("O","N","C"))
ef1b$position<-factor(ef1b$position, levels=c("phrase-medial","phrase-final"))
ef1b$tenseness<-factor(ef1b$tenseness, levels=c("tense","lax"))
ef1b$pos.tense<-paste0(ef1b$tenseness," ",ef1b$position)
ef1b$pos.tense<-as.factor(ef1b$pos.tense)

(ef1b.plot<-
  ef1b %>% ggplot(aes(x=sylpos, y=fit,col=position,group=position,lty=position)) +
  geom_point() +
  geom_line()+
  # geom_line(aes(group=1),col="red")+
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=.8) +
  facet_grid(~tenseness)+
  xlab("Syllable position")+
  scale_color_brewer(palette = "Set1")+
  ylab("Fit for duration (ms)")+
  ylim(30,160)+
  theme_bw(base_size=10)+
  theme(legend.position = "top",legend.title = element_blank()))

(ef1b.plot.pos.tens<-
    ef1b %>% ggplot(aes(x=sylpos, y=fit,col=pos.tense,group=pos.tense,lty=tenseness)) +
    geom_point() +
    geom_line()+
    # geom_line(aes(group=1),col="red")+
    geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=.8) +
    # facet_grid(~tenseness)+
    xlab("Syllable position")+
    scale_color_brewer(palette = "Set1")+
    ylab("Fit for duration (ms)")+
    ylim(30,160)+
    theme_bw(base_size=10)+
    theme(legend.position = "top",legend.title = element_blank()))

#' ## Disyllabic target words
#### Disyllabic acoustic target words ####

d1<-lmer(sylposdur ~ position + sylpos + tenseness + (1|subject) + (1|word), data = di)
qqPlot(resid(d1))
d2<-lmer(sylposdur ~ position * sylpos + tenseness + (1|subject) + (1|word), data = di)
anova(d1,d2)
d3<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word), data = di)
anova(d2,d3)
d4<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
           (0+tenseness|subject) + (0+position+tenseness+sylpos|word), data = di)
anova(d3,d4,refit=F)
#best model is d4, models with additional random slopes do not converge
qqPlot(resid(d4))
abline(h = cut<- c(-50,60))
di2<-droplevels(subset(di, resid(d4)<cut[2] & resid(d4)> cut[1]))

di2$sylpos <- factor(di2$sylpos, levels = c("O1","N1","O2","N2"))
d5<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
             (0+tenseness|subject) + (0+position+tenseness+sylpos|word), data = di2)
d5.r2dummy<-lmer(sylposdur ~ position * sylpos * tenseness +
               (1+tenseness|subject) + (1+position+tenseness+sylpos|word), data = di2)
anova(d5, d5.r2dummy, refit=F)
# after cut model does not converge any more
qqPlot(resid(d5))
summary(d5)

#' ### Excursus: Normalization
#' Normalization of dependent variable for speech rate does not change the effects considerably
d5.norm<-lmer(sylposdur ~ position * sylpos * tenseness + (1|subject) + (1|word) +
                (0+tenseness|subject) + (0+position+tenseness+sylpos|word), data = di2)
plot(allEffects(d5.norm))
summary(d5.norm)
###


#' ### Post-hoc tests (multiple comparisons, Tukey-corrected):
d5em <- emmeans(d5, pairwise~position*sylpos*tenseness)
d5em.tab<-d5em$contrasts %>%
  summary(infer = T)
d5em.tab$p.value<-lapply(d5em.tab$p.value, FUN = function (x) {
  # levels of p-values
  if (x>=0.05)  p <- "n.\\,s." # not significant
  if (x<0.05)  p <- "<0.05"  # significant
  if (x<=0.01)  p <- "<0.01"  # significant
  if (x<=0.001) p <- "<0.001" # significant
  return(p)
})
d5em.tab<-as.data.frame(d5em.tab)
colnames(d5em.tab)<-c("Comparison","$\\beta$","s","df","Lower CI","Upper CI","t","p")
d5em.tab.sig<-subset(d5em.tab, p!="n.\\,s.")

rownames(d5em.tab.sig)<-seq(1, nrow(d5em.tab.sig), 1)

#' ### Subset model acoustic disyllabic for lax only
di.lax<-droplevels(subset(di, tenseness=="lax"))
dl1<-lmer(sylposdur ~ position + sylpos  + (1|subject) + (1|word), data = di.lax)
dl2<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word), data = di.lax)
anova(dl1,dl2)
dl3<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = di.lax)
anova(dl2,dl3,refit=F) #better
qqPlot(resid(l3))
abline(h = cut <- 50)
di.lax.sub<-droplevels(subset(di.lax, resid(dl3) < cut))
dl4<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position+sylpos|word), data = di.lax.sub)
qqPlot(resid(dl4))
abline(h = cut <- -50)
di.lax.sub<-droplevels(subset(di.lax.sub, resid(dl4) > cut))
dl5<-lmer(sylposdur ~ position * sylpos  + (1|subject) + (1|word) + (0+position|word), data = di.lax.sub)
qqPlot(resid(dl5))
plot(allEffects(dl5))
summary(dl5)
lattice::dotplot(ranef(dl5))


#' ### Post-hoc tests (multiple comparisons, Tukey-corrected):
dl5em <- emmeans(dl5, pairwise~position*sylpos)
dl5em.tab<-dl5em$contrasts %>%
  summary(infer = T)
dl5em.tab$p.value<-lapply(dl5em.tab$p.value, FUN = function (x) {
  # levels of p-values
  if (x>=0.05)  p <- "n.\\,s." # not significant
  if (x<0.05)  p <- "<0,05"  # significant
  if (x<=0.01)  p <- "<0,01"  # significant
  if (x<=0.001) p <- "<0,001" # significant
  return(p)
})
dl5em.tab<-as.data.frame(dl5em.tab)
colnames(dl5em.tab)<-c("Comparison","$\\beta$","s","df","Lower CI","Upper CI","t","p")
dl5em.tab.sig<-subset(dl5em.tab, p!="n.\\,s.")
dl5em.tab.sig


#' ### Disyllabic effect plot
#' Preparation
ef2b<-as.data.frame(effect("position*sylpos*tenseness",d5))
ef2b$sylpos<-factor(ef2b$sylpos, levels=c("O1","N1","O2","N2"))
ef2b$position<-factor(ef2b$position, levels=c("phrase-medial","phrase-final"))
ef2b$tenseness<-factor(ef2b$tenseness, levels=c("tense","lax"))

(ef2b.plot<-
  ef2b %>% ggplot(aes(x=sylpos, y=fit,col=position,group=position,lty=position)) +
  geom_point() +
  geom_line()+
  # geom_line(aes(group=1),col="red")+
  geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
  facet_grid(~tenseness)+
  scale_color_brewer(palette = "Set1")+
  xlab("Syllable position")+
  ylab("")+
  ylim(30,160)+
  theme_bw(base_size=10)+
  theme(legend.position = "top",legend.title = element_blank()))

#' ## Figure 3: Mono- and disyllabic effect plot
ggpubr::ggarrange(ef1b.plot, ef2b.plot, common.legend = T,
                  labels = c("a)","b)"),font.label = list(size=11,face="plain"),align = "h",nrow=1)
ggsave(file="mono-di-acoustic-effects.pdf", width = 170, height = 100, units = "mm")

#' ## Table 3 and 4
x<-capture.output(
   texreg::texreg(list(texreg::extract(m5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(d5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(mf3.1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   ),
   digits = 1,booktabs = T,custom.model.names = c("Duration in monosyllabic words"),
   float.pos = "!htbp", label = "tab:acousticresults.mono", use.packages = F,single.row = T,caption.above = T,
   caption="Model predictions for the acoustic duration of syllable position in  monosyllabic words
   (standards error in parentheses) for the fixed effects  phrasal position (phrase-medial vs. phrase-final),
   syllable constituent (\\emph{sylcon}, O = onset, N = nucleus, C = coda), and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "scriptsize",center=F)
 )
 x[5]<-paste0(x[5], "\\centering")
 x[27]<-paste0(x[27],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(m5.r2dummy)[1],2),"/",round(r.squaredGLMM(m5.r2dummy)[2],2),
               # "&", round(r.squaredGLMM(mf3.1)[1],2), "/", round(r.squaredGLMM(mf3.1)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")


 x<-capture.output(
   texreg::texreg(list(texreg::extract(d5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(d5, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
                       # texreg::extract(mf3.1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   ),
   digits = 1,booktabs = T,custom.model.names = c("Duration in disyllabic words"),
   float.pos = "!htbp", label = "tab:acousticresults.di", use.packages = F,single.row = T,caption.above = T,
   caption="Model predictions for the acoustic duration of syllable position in disyllabic words
   (standards error in parentheses) for the fixed effects  phrasal position (phrase-medial vs. phrase-final),
   syllable position (\\emph{sylcon}, O1/O2 = first/second onset, N1/N2 = first/second nucleus), and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "scriptsize",center=F)
 )
 x[5]<-paste0(x[5], "\\centering")
 x[31]<-paste0(x[31],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(d5.r2dummy)[1],2),"/", round(r.squaredGLMM(d5.r2dummy)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")




#' ## Table 9 and 10: Post-hoc tests (multiple comparisons, Tukey-corrected):
 m5em <- emmeans(m5, pairwise~position*sylpos*tenseness)
 m5em.tab<-m5em$contrasts %>%
   summary(infer = T)
 m5em.tab$p.value<-lapply(m5em.tab$p.value, FUN = function (x) {
   # levels of p-values
   if (x>=0.05)  p <- "n.\\,s." # not significant
   if (x<0.05)  p <- "<0.05"  # significant
   if (x<=0.01)  p <- "<0.01"  # significant
   if (x<=0.001) p <- "<0.001" # significant
   return(p)
 })
 m5em.tab<-as.data.frame(m5em.tab)
 colnames(m5em.tab)<-c("Comparison","$\\beta$","s","df","Lower CI","Upper CI","t","p")
 m5em.tab.sig<-subset(m5em.tab, p!="n.\\,s.")
 m5em.tab.sig<-m5em.tab.sig[-grep(pattern = ",O,",x = m5em.tab.sig$Comparison),]

 rownames(m5em.tab.sig)<-seq(1,nrow(m5em.tab.sig),1)

print(xtable::xtable(m5em.tab.sig,label="tab:acoustic.mono.post",
                      caption = "Significant contrasts for monosyllabic target words with duration as dependent variable,
                      the predictors position (phrase-medial/-final), syllable position (O = Onset, N = Nucleus, C = Coda),
                      and tenseness (tense/lax), as well as  estimated value $\\beta$, standard deviation (s), degrees of freedom (df), 95\\,\\% confidence
                      intervals, and t-values, Tukey-corrected for  multiple comparisons.",
                      align = "llrrrrrrr",digits = c(0,0,1,1,1,2,2,2,2)),
       include.colnames = T,include.rownames=T,booktabs = T, size = "\\scriptsize",
       caption.placement = "top", sanitize.text.function = function(x) x)

 print(xtable::xtable(d5em.tab.sig,label="tab:acoustic.di.post",
                      caption = "Significant contrasts for disyllabic target words with duration as dependent variable,
                      the predictors position (phrase-medial/-final), syllable position (O1/O2 = first/second onset, N1/N2 = first/second nucleus, C = coda),
                      and tenseness (tense/lax), as well as  estimated value $\\beta$, standard deviation (s), degrees of freedom (df), 95\\,\\% confidence
                      intervals, and t-values, Tukey-corrected for  multiple comparisons.",
                      align = "llrrrrrrr",digits = c(0,0,1,1,1,2,2,2,2)),
       include.colnames = T,include.rownames=T,booktabs = T, size = "\\scriptsize",
       caption.placement = "top", sanitize.text.function = function(x) x)




#' # Articulatory analysis
# Articulatory analysis ####
# read mono- and disyllabic data
 mono.art<-read.table("mono.art.txt", sep="\t", fileEncoding = "utf-8")
 mono.art$tenseness<-factor(mono.art$tenseness,levels=c("tense","lax"))
 mono.art$position<-factor(mono.art$position,levels=c("phrase-medial","phrase-final"))

 di.art<-read.table("di.art.txt", sep="\t", fileEncoding = "utf-8")
 di.art$tenseness<-factor(di.art$tenseness,levels=c("tense","lax"))
 di.art$position<-factor(di.art$position,levels=c("phrase-medial","phrase-final"))


 # data contains the following columns:
 # 1: Source name
 # 2: Trajectory name
 # 3: Comment
 # 4: Gesture onset at time point
 # 5: Peak velocity of closing at time point
 # 6: nucleus (plateau) onset at time point
 # 7: Maximum constriction at time point
 # 8: nucleus (plateau) offset at timepoint
 # 9: Peak velocity of opening at time point
 # 10: Gesture offset at time point
 # 11: Gesture onset of anterior-posterior dimension at position
 # 12: Peak velocity of closing of anterior-posterior dimension at position
 # 13: nucleus (plateau) onset of anterior-posterior dimension at position
 # 14: Maximum constriction of anterior-posterior dimension at position
 # 15: nucleus (plateau) offset of anterior-posterior dimension at position
 # 16: Peak velocity of closing of anterior-posterior dimension at position
 # 17: Gesture offset of anterior-posterior dimension at position
 # 18: Gesture onset of inferior-superior dimension at position
 # 19: Peak velocity of closing of inferior-superior dimension at position
 # 20: nucleus (plateau) onset of inferior-superior dimension at position
 # 21: Maximum constriction of inferior-superior dimension at position
 # 22: nucleus (plateau) offset of inferior-superior dimension at position
 # 23: Peak velocity of closing of inferior-superior dimension at position
 # 24: Gesture offset of inferior-superior dimension at position
 # 25-31: velocities in cm/sec
 # 32: subject id
 # 33: sex
 # 34: phrasal position
 # 35: target word
 # 36: repetition
 # 37: plateau duration (platdur)
 # 38: opening duration (opdur)
 # 39: closing duration (closdur)
 # 40: opening displacement (opdisp)
 # 41: closing displacement (closdisp)


#' ## Data removal and data used
#' Remove all items with closing or opening displacement below 1mm
#' **Monosyllabic**
#' 389 target words were originally included:
 addmargins(table(mono.art$word,mono.art$position))
#' **Disyllabic**
#' 198 target words were originally included:
 addmargins(table(di.art$word,di.art$position))

 ggplot(mono.art, aes(x=word, y=closdisp))+
   geom_point()+
   geom_hline(yintercept=1)
 ggplot(mono.art, aes(x=word, y=opdisp))+
   geom_point()+
   geom_hline(yintercept=1)
 ggplot(di.art, aes(x=word, y=opdisp))+
   geom_point()+
   geom_hline(yintercept=1)

 mono.art.norm<- mono.art[mono.art$opdisp>1 & mono.art$closdisp>1,]
 nrow(mono.art)-nrow(mono.art.norm)
 #' 748 stimuli after data exclusion:
 addmargins(table(mono.art.norm$word,mono.art.norm$position))

 # nrow(di.art[di.art$opdisp<1 | di.art$closdisp<1 ,])
 di.art.norm<- di.art[di.art$opdisp>1 & di.art$closdisp>1,]
 nrow(di.art)-nrow(di.art.norm)

 # invalid data of closing duration outlier, remove:
 di.art.norm<-di.art.norm[di.art.norm$closdur>15,]
 nrow(di.art)-nrow(di.art.norm)

#' 385 stimuli after data exclusion:
 addmargins(table(di.art.norm$word,di.art.norm$position))

#' **Total data excluded in the articulatory analysis **
 # nrow(mono.art)/nrow(mono.art.norm)-1
#' **Monosyllabic**
 sprintf("%1.2f%%",round((1- nrow(mono.art.norm)/800)*100),2)
#' **Disyllabic**
 #nrow(di.art)/nrow(di.art.norm)-1
 sprintf("%1.2f%%",round((1- nrow(di.art.norm)/400)*100),2)
#' ** Total total **
 sprintf("%1.2f%%",round((1- ((nrow(di.art.norm)+nrow(mono.art.norm))/(400+800)) )*100,2))


#' Number of target words used in the articulatory analysis
 addmargins(table(rbind(mono.art.norm,di.art.norm)$word, rbind(mono.art.norm,di.art.norm)$position))

 mono.art$wordpair<-NULL
 mono.art$wordpair[mono.art$word=="Bahn" | mono.art$word=="Bann"] <-"Bahn/Bann"
 mono.art$wordpair[mono.art$word=="Beet" | mono.art$word=="Bett"] <-"Beet/Bett"
 mono.art$wordpair[mono.art$word=="Ruhm" | mono.art$word=="Rum"] <-"Ruhm/Rum"
 mono.art$wordpair[mono.art$word=="Stiel" | mono.art$word=="still"] <-"Stiel/still"


 di.art$wordpair<-NULL
 di.art$wordpair[di.art$word=="Huete" | di.art$word=="Huette"] <-"Huete/Huette"
 di.art$wordpair[di.art$word=="Miete" | di.art$word=="Mitte"] <-"Miete/Mitte"

#' ## Monosyllabic target words
 # mono.art.norm$wordpair<-NULL
 # mono.art.norm$wordpair[mono.art.norm$word=="Bahn" | mono.art.norm$word=="Bann"] <-"Bahn/Bann"
 # mono.art.norm$wordpair[mono.art.norm$word=="Beet" | mono.art.norm$word=="Bett"] <-"Beet/Bett"
 # mono.art.norm$wordpair[mono.art.norm$word=="Ruhm" | mono.art.norm$word=="Rum"] <-"Ruhm/Rum"
 # mono.art.norm$wordpair[mono.art.norm$word=="Stiel" | mono.art.norm$word=="still"] <-"Stiel/still"


#' ### Closing duration
 # Closing duration mono ####
 clodur1<-lmer(log(closdur) ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art)
 qqPlot(resid(clodur1))
 abline(h = cut <- .6)
 mono.art.sub1<-droplevels(subset(mono.art, resid(clodur1) < cut))
 clodur2<-lmer(log(closdur) ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art.sub1)
 qqPlot(resid(clodur2))
 summary(clodur2)
 plot(allEffects(clodur2))

 clodur3.int<-lmer(log(closdur) ~ position  * tenseness + (1|subject) + (1|word), data = mono.art.sub1)
 anova(clodur2, clodur3.int) # not better
 # failed

 clodur2.ms<-lmer(closdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art.sub1)
 summary(clodur2.ms)
 plot(allEffects(clodur2.ms))

 #subset lax
 lax.cd.mono<-mono.art.sub1 %>% filter(tenseness == "lax") %>% droplevels()
 cd.lax.mono.1<-lmer(log(closdur) ~ position   + (1|subject)  + (1|word), data = lax.cd.mono)
 qqPlot(resid(cd.lax.mono.1))
 plot(allEffects(cd.lax.mono.1))
 summary(cd.lax.mono.1)

#' ### Plateau duration
 platdur1<-lmer(log(platdur) ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art)
 qqPlot(resid(platdur1))
 abline(h = cut <- c(-1.1,1.3))
 mono.art.plat1<-droplevels(subset(mono.art, resid(platdur1) > cut[1] & resid(platdur1) < cut[2]))
 platdur2<-lmer(log(platdur) ~ position  + tenseness + (1|subject) + (1|word), data = mono.art.plat1)
 qqPlot(resid(platdur2))
 summary(platdur2)
 plot(allEffects(platdur2))
 platdur3<-lmer(log(platdur) ~ position  * tenseness + (1|subject) + (1|word), data = mono.art.plat1)
 anova(platdur2,platdur3) #no


 platdur2.ms<-lmer(platdur ~ position  + tenseness + (1|subject) + (1|word), data = mono.art.plat1)
 summary(platdur2.ms)

#' ### Closing displacement
 closdis1<-lmer(closdisp ~ position  + tenseness +(1+position+tenseness|subject) + (1|word), data = mono.art)
 qqPlot(resid(closdis1))
 abline(h = cut <- 10)
 mono.art.disp1<-droplevels(subset(mono.art, resid(closdis1) < cut))
 closdis2<-lmer(closdisp ~ position  + tenseness +(1+position|subject) + (1|word), data = mono.art.disp1)
 qqPlot(resid(closdis2))
 summary(closdis2)
 plot(allEffects(closdis2))
 # final is more displaced than medial


#' ### Closing velocity
 closvel1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = mono.art)
 qqPlot(resid(closvel1))
 abline(h = cut <- 18)
 mono.art.vel1<-droplevels(subset(mono.art, resid(closvel1) < cut))
 closvel2<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1+position|subject) + (1|word), data = mono.art.vel1)
 qqPlot(resid(closvel2))
 summary(closvel2)
 plot(allEffects(closvel2))
 # dotplot(ranef(closvel2))
 # final is sign. slower than medial


#' In none of the models with position and tenseness as fixed and subjects and words
#' as random effects does tenseness affect the outcome variable

#' ## Articulatory disyllabic analysis ####
#' ### Closing duration
 di.clodur1<-lmer(closdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = di.art)
 qqPlot(resid(di.clodur1))
 abline(h = cut <- 20)
 di.art.sub1<-droplevels(subset(di.art, resid(di.clodur1) < cut))
 di.clodur2<-lmer(closdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = di.art.sub1)
 di.clodur3<-lmer(closdur ~ position  + tenseness +
                    (1+position+tenseness|subject) + (1|word), data = di.art.sub1)
 anova(di.clodur2,di.clodur3,refit=F)
 qqPlot(resid(di.clodur2))
 summary(di.clodur2)
 plot(allEffects(di.clodur2))
 di.clodur.int<-lmer(closdur ~ position  * tenseness + (1+position+tenseness|subject) + (1|word), data = di.art.sub1)
 anova(di.clodur2,di.clodur.int) # not better
 summary(di.clodur.int)
 # longer closing duration in final position


#' ### Plateau duration di
 di.platdur1<-lmer(platdur ~ position  + tenseness + (1+position+tenseness|subject) + (1|word), data = di.art)
 qqPlot(resid(di.platdur1))
 abline(h = cut <- 30)
 di.art.plat1<-droplevels(subset(di.art, resid(di.platdur1)  < cut))
 di.platdur2<-lmer(platdur ~ position  + tenseness + (1+position|subject) + (1|word), data = di.art.plat1)
 qqPlot(resid(di.platdur2))
 summary(di.platdur2)
 di.platdur.int<-lmer(platdur ~ position  * tenseness + (1+position|subject) + (1|word), data = di.art.plat1)
 anova(di.platdur2,di.platdur.int) # not better
 summary(di.platdur.int)
 # final has longer plateau phase than medial


#' ### Closing displacement di
 di.closdis1<-lmer(closdisp ~ position  + tenseness + (1+position|subject) + (1|word), data = di.art)
 qqPlot(resid(di.closdis1))
 summary(di.closdis1)
 plot(allEffects(di.closdis1))

 di.closdis.int<-lmer(closdisp ~ position  * tenseness +
                        (1+position+tenseness|subject) + (1|word), data = di.art)
  anova(di.closdis1,di.closdis.int) # interaction for displacement better, but only slightly
# warrants closer inspection for differences in vowel height
# final is less (!) displaced than medial, other than for mono-words




#' ### Closing velocity di
 di.closvel1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +  (1+position+tenseness|subject) + (1|word), data = di.art)
 qqPlot(resid(di.closvel1))
 summary(di.closvel1)
 plot(allEffects(di.closvel1))
 di.closvel.int<-lmer(PVEL..V.cm.sec. ~ position  * tenseness + (1|subject) + (1|word), data = di.art)
 summary(di.closvel.int)
 anova(di.closvel1,di.closvel.int) # not better
 plot(allEffects(di.closvel.int))

 di.closvel.tenseness.med<-lmer(PVEL..V.cm.sec. ~ tenseness + (1|subject) + (1|word), data = di.art %>% filter(position=="phrase-medial"))
 summary(di.closvel.tenseness.med)
 di.closvel.tenseness.fin<-lmer(PVEL..V.cm.sec. ~ tenseness + (1|subject) + (1|word), data = di.art %>% filter(position=="phrase-final"))
 summary(di.closvel.tenseness.fin)

 di.closvel.position.tense<-lmer(PVEL..V.cm.sec. ~ position + (1|subject) + (1|word), data = di.art %>% filter(tenseness=="tense"))
 summary(di.closvel.position.tense)
 di.closvel.position.lax<-lmer(PVEL..V.cm.sec. ~ position + (1|subject) + (1|word), data = di.art %>% filter(tenseness=="lax"))
 summary(di.closvel.position.lax)

# final is sign. slower than medial

#' ## Follow-up-models for mono-/disyllabic articulatory displacement and velocity

closdis.bahn<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                  data = droplevels(subset(mono.art, wordpair == "Bahn/Bann")))
 qqPlot(resid(closdis.bahn))
plot(allEffects(closdis.bahn))
summary(closdis.bahn)
# for Bahn/Bann no effect of position, but lax vowels are less displaced
closvel.bahn<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1+position|subject),
                   data = droplevels(subset(mono.art, wordpair == "Bahn/Bann")))
qqPlot(resid(closvel.bahn))
plot(allEffects(closvel.bahn))
summary(closvel.bahn)
# velocity is slower in final position, no effect for laxness


vh.beet<-droplevels(subset(mono.art, wordpair == "Beet/Bett"))
closdis.beet<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                   data = vh.beet)
qqPlot(resid(closdis.beet))
abline(h = cut <- -4)
vh.beet1<-droplevels(subset(vh.beet, resid(closdis.beet) > cut))
closdis.beet1<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                    data = vh.beet1)
qqPlot(resid(closdis.beet1))
plot(allEffects(closdis.beet1))
summary(closdis.beet1)
# for Beet/Bett no effect of position, but lax vowels are more displaced
closvel.beet1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                    data = vh.beet1)
qqPlot(resid(closvel.beet1))
plot(allEffects(closvel.beet1))
summary(closvel.beet1)
closvel.beet2<-lmer(PVEL..V.cm.sec. ~ position  * tenseness +(1|subject),
                    data = vh.beet1)
anova(closvel.beet1,closvel.beet2) # no sign. interaction
# final position is slower; lax vowels lead to faster velocity (consistent with greater displacement)


vh.ruhm<-droplevels(subset(mono.art, wordpair == "Ruhm/Rum"))
closdis.ruhm<-lmer(closdisp ~ position  + tenseness +(1|subject),
                   data = vh.ruhm)
qqPlot(resid(closdis.ruhm))
plot(allEffects(closdis.ruhm))
summary(closdis.ruhm)
# for Ruhm/Rum no effect of position, but lax vowels are more displaced
closvel.ruhm<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                   data = vh.ruhm)
qqPlot(resid(closvel.ruhm))
abline(h = cut <- 6)
vh.ruhm.vel<-droplevels(subset(vh.ruhm, resid(closvel.ruhm) < cut))
closvel.ruhm1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                   data = vh.ruhm.vel)
plot(allEffects(closvel.ruhm1))
summary(closvel.ruhm1)
closvel.ruhm1.int<-lmer(PVEL..V.cm.sec. ~ position  * tenseness +(1|subject),
                    data = vh.ruhm.vel)
anova(closvel.ruhm1,closvel.ruhm1.int) # interaction not better


vh.stiel<-droplevels(subset(mono.art, wordpair == "Stiel/still"))
closdis.stiel<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                   data = vh.stiel)
qqPlot(resid(closdis.stiel))
abline(h = cut <- 5)
vh.stiel1<-droplevels(subset(vh.stiel, resid(closdis.stiel) < cut))
closdis.stiel1<-lmer(closdisp ~ position  + tenseness +(1+position|subject),
                    data = vh.stiel1)
plot(allEffects(closdis.stiel1))
summary(closdis.stiel1)
# for Stiel/still more displacement in final position and less displacement after lax vowels
closvel.stiel<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1+position|subject),
                    data = vh.stiel)
qqPlot(resid(closvel.stiel))
abline(h = cut <- 10)
vh.stiel.vel1<-droplevels(subset(vh.stiel, resid(closvel.stiel) < cut))
closvel.stiel1<-lmer(PVEL..V.cm.sec. ~ position  + tenseness +(1|subject),
                     data = vh.stiel.vel1)
qqPlot(resid(closvel.stiel1))
plot(allEffects(closvel.stiel1))
summary(closvel.stiel1)

#' ## Table 6
# part 1
x<-capture.output(
  texreg::texreg(list(texreg::extract(closdis.bahn, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(closdis.beet1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(closdis.ruhm, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                      texreg::extract(closdis.stiel1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
  ),
  digits = 1,booktabs = T,custom.model.names = c("Bahn/Bann","Beet/Bett","Ruhm/Rum","Stiel/still"),
  float.pos = "!htbp", label = "tab:dispvel.mono", use.packages = F,single.row = T,caption.above = T,
  caption="LME model predictions for the articulatory measurements of closing gesture displacement (ClosDisp) and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for   phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax) in all monosyllabic target words.",
  dcolumn = T, fontsize = "footnotesize",center=F)
)
x[3]<-paste0(x[3], "\\centering")
x[15]<-paste0(x[15],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(closdis.bahn)[1],2),"/",round(r.squaredGLMM(closdis.bahn)[2],2),
              "&", round(r.squaredGLMM(closdis.beet1)[1],2),"/", round(r.squaredGLMM(closdis.beet1)[2],2),
              "&", round(r.squaredGLMM(closdis.ruhm)[1],2), "/", round(r.squaredGLMM(closdis.ruhm)[2],2),
              "&", round(r.squaredGLMM(closdis.stiel1)[1],2), "/", round(r.squaredGLMM(closdis.stiel1)[2],2),
              "\\\\ \n")
cat(x, sep = "\n")

# part 2
  x<-capture.output(
    texreg::texreg(list(texreg::extract(closvel.bahn, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                        texreg::extract(closvel.beet1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                        texreg::extract(closvel.ruhm1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                        texreg::extract(closvel.stiel1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
    ),
    digits = 1,booktabs = T,custom.model.names = c("Bahn/Bann","Beet/Bett","Ruhm/Rum","Stiel/still"),
    float.pos = "!htbp", label = "tab:dispvel.mono", use.packages = F,single.row = T,caption.above = T,
    caption="LME model predictions for the articulatory measurements of closing gesture displacement (ClosDisp) and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for   phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax) in all monosyllabic target words.",
    dcolumn = T, fontsize = "footnotesize",center=F)
  )
  x[3]<-paste0(x[3], "\\centering")
  x[15]<-paste0(x[15],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(closvel.bahn)[1],2),"/",round(r.squaredGLMM(closvel.bahn)[2],2),
                "&", round(r.squaredGLMM(closvel.beet1)[1],2),"/", round(r.squaredGLMM(closvel.beet1)[2],2),
                "&", round(r.squaredGLMM(closvel.ruhm)[1],2), "/", round(r.squaredGLMM(closvel.ruhm)[2],2),
                "&", round(r.squaredGLMM(closvel.stiel1)[1],2), "/", round(r.squaredGLMM(closvel.stiel1)[2],2),
                "\\\\ \n")
  cat(x, sep = "\n")

vh.di.miete<-droplevels(subset(di.art, wordpair=="Miete/Mitte"))
disp.miete<-lmer(closdisp ~ position  + tenseness + (1|subject), data = vh.di.miete)
qqPlot(resid(disp.miete))
abline(h = cut <- 4)
vh.di.miete.sub<-subset(vh.di.miete, resid(disp.miete) < cut)
disp.miete1<-lmer(closdisp ~ position  + tenseness + (1|subject), data = vh.di.miete.sub)
qqPlot(resid(disp.miete1))
summary(disp.miete1)
# marginal effect for position; lax has greater displacement
vel.miete<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1|subject), data = vh.di.miete)
qqPlot(resid(vel.miete))
plot(allEffects(vel.miete))
summary(vel.miete)
# slower velocity in final position; faster velocity after lax vowels


vh.di.huete<-droplevels(subset(di.art, wordpair=="Huete/Huette"))
disp.huete<-lmer(closdisp ~ position  + tenseness + (1|subject), data = vh.di.huete)
qqPlot(resid(disp.huete))
summary(disp.huete)
# no effect for position; lax has greater displacement
vel.huete<-lmer(PVEL..V.cm.sec. ~ position  + tenseness + (1|subject), data = vh.di.huete)
qqPlot(resid(vel.huete))
plot(allEffects(vel.huete))
summary(vel.huete)
# slower velocity in final position; faster velocity after lax vowels


#' ## Fig 6: Articulatory effect plots #####
ef.m.closdur.pos<-as.data.frame(effect("position",clodur2.ms)) %>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
levels(ef.m.closdur.pos$position)<-c("medial","final")
ef.m.closdur.ten<-as.data.frame(effect("tenseness",clodur2))%>% mutate(tenseness= factor(tenseness, levels = c("tense","lax")))

 ef.m.platdur.pos<-as.data.frame(effect("position",platdur2.ms)) %>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
 levels(ef.m.platdur.pos$position)<-c("medial","final")
 ef.m.platdur.ten<-as.data.frame(effect("tenseness",platdur2))%>% mutate(tenseness= factor(tenseness, levels = c("tense","lax")))

 ef.m.closdis.pos<-as.data.frame(effect("position",closdis2)) %>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
 levels(ef.m.closdis.pos$position)<-c("medial","final")
 ef.m.closdis.ten<-as.data.frame(effect("tenseness",closdis2))%>% mutate(tenseness= factor(tenseness, levels = c("tense","lax")))

 ef.m.closvel.pos<-as.data.frame(effect("position",closvel2)) %>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
 levels(ef.m.closvel.pos$position)<-c("medial","final")
 ef.m.closvel.ten<-as.data.frame(effect("tenseness",closvel2))%>% mutate(tenseness= factor(tenseness, levels = c("tense","lax")))

 # disyllabic effect plot preparation
 ef.d.closdur.pos<-as.data.frame(effect("position",di.clodur2))%>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
 levels(ef.d.closdur.pos$position)<-c("medial","final")
 ef.d.platdur.pos<-as.data.frame(effect("position",di.platdur2))%>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
 levels(ef.d.platdur.pos$position)<-c("medial","final")
 ef.d.closdis.pos<-as.data.frame(effect("position",di.closdis1))%>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
 levels(ef.d.closdis.pos$position)<-c("medial","final")
 ef.d.closvel.pos<-as.data.frame(effect("position",di.closvel1))%>% mutate(position= factor(position, levels = c("phrase-medial","phrase-final")))
 levels(ef.d.closvel.pos$position)<-c("medial","final")

 plevel <- function (x) {
   if (x>0.05 )  p <- paste0("p = ",round(x,1)) # not significant
   if (x<=0.05)  p <- "p < 0.05"  # significant
   if (x<=0.01)  p <- "p < 0.01"  # significant
   if (x<=0.001) p <- "p < 0.001" # significant
   return(p)
 }

 (ef.m.closdur.pos.plot<-
     ef.m.closdur.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=4.1,label=plevel(summary(clodur2)$coefficients[1,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Mono: ",plevel(summary(clodur2)$coefficients[2,5])))+
     ylab("Closing duration (ms)")+
     ylim(40,80)+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))
 (ef.m.closdur.ten.plot<-
     ef.m.closdur.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("Closing duration (log)")+
     ylab("")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))

 (ef.m.platdur.pos.plot<-
     ef.m.platdur.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=4.1,label=plevel(summary(platdur2)$coefficients[1,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Mono: ",plevel(summary(platdur2)$coefficients[2,5])))+
     ylab("Plateau duration (ms)")+
     ylim(20,120)+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))
 (ef.m.platdur.ten.plot<-
     ef.m.platdur.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("")+
     ylab("Plateau duration (log)")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))

 (ef.m.closdis.pos.plot<-
     ef.m.closdis.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=9,label=plevel(summary(closdis2)$coefficients[2,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Mono: ",plevel(summary(closdis2)$coefficients[2,5])))+
     ylab("Closing displacement (mm)")+
     ylim(5,11)+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))
 (ef.m.closdis.ten.plot<-
     ef.m.closdis.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("")+
     ylab("Closing displacement (mm)")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))

 (ef.m.closvel.pos.plot<-
     ef.m.closvel.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=22,label=plevel(summary(closvel2)$coefficients[2,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Mono: ",plevel(summary(closvel2)$coefficients[2,5])))+
     ylab("Closing velocity (cm/s)")+
     ylim(12,28)+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))
 (ef.m.closvel.ten.plot<-
     ef.m.closvel.ten %>%
     ggplot(aes(x=tenseness, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     scale_color_brewer(palette = "Set1")+
     xlab("")+
     ylab("Closing velocity (cm/s)")+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))


 # disyllabic effect plots
 (ef.d.closdur.pos.plot<-
     ef.d.closdur.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=55,label=plevel(summary(di.clodur2)$coefficients[2,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Di: ",plevel(summary(di.clodur2)$coefficients[2,5])))+
     ylab("Closing duration (ms)")+
     ylim(c(40,80))+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))
 (ef.d.platdur.pos.plot<-
     ef.d.platdur.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=43,label=plevel(summary(di.platdur2)$coefficients[2,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Di: ",plevel(summary(di.platdur2)$coefficients[2,5])))+
     ylab("Plateau duration (ms)")+
     ylim(c(20,120))+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))
 (ef.d.closdis.pos.plot<-
     ef.d.closdis.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=7.1,label=plevel(summary(di.closdis1)$coefficients[2,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Di: ",plevel(summary(di.closdis1)$coefficients[2,5])))+
     ylab("Closing displacement (mm)")+
     ylim(5,11)+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))
 (ef.d.closvel.pos.plot<-
     ef.d.closvel.pos %>%
     ggplot(aes(x=position, y=fit)) +
     geom_point() +
     geom_line(aes(group=1),col="red")+
     geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.4,alpha=.8,size=0.8) +
     # annotate("text",x=1.5,y=22,label=plevel(summary(di.closvel1)$coefficients[2,5]))+
     scale_color_brewer(palette = "Set1")+
     xlab(paste0("Di: ",plevel(summary(di.closvel1)$coefficients[2,5])))+
     ylab("Closing velocity (cm/s)")+
     ylim(12,28)+
     theme_bw(base_size=10)+
     theme(legend.position = "top",legend.title = element_blank()))


 ggpubr::ggarrange(ef.m.closdur.pos.plot, ef.m.platdur.pos.plot, ef.m.closdis.pos.plot, ef.m.closvel.pos.plot,
                   # ef.m.closdur.ten.plot, ef.m.platdur.ten.plot, ef.m.closdis.ten.plot, ef.m.closvel.ten.plot,
                   ef.d.closdur.pos.plot, ef.d.platdur.pos.plot, ef.d.closdis.pos.plot, ef.d.closvel.pos.plot,
                   common.legend = T,
                   labels = c("a)","b)","c)","d)","e)","f)","g)","h)"),
                   font.label = list(size=11,face="plain"),align = "h",nrow=2,ncol=4)
 ggsave(file="mono-art-effects.pdf", width = 200, height = 130, units = "mm")

#' ## Table 5
 x<-capture.output(
   texreg::texreg(list(texreg::extract(clodur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(platdur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(closdis2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(closvel2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   )
   ,
   override.pvalues = list(summary(clodur2)$coefficients[1:3,5],
                           summary(platdur2)$coefficients[1:3,5],
                           summary(closdis2)$coefficients[1:3,5],
                           summary(closvel2)$coefficients[1:3,5])
   ,
   digits = 1,booktabs = T,custom.model.names = c("ClosDur (log)","PlatDur (log)","ClosDisp","ClosVel"),
   float.pos = "!htbp", label = "tab:articresults.mono", use.packages = F,single.row = T,caption.above = T,
   caption="LME model predictions for the articulatory measurements of closing gesture duration  (ClosDur),    plateau duration (PlatDur), closing gesture displacement (ClosDisp), and closing gesture velocity (ClosVel)    of the last consonant in all monosyllabic word pairs (standard error in parentheses)    for the fixed effects  phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "footnotesize",center=F)
 )
 x[3]<-paste0(x[3], "\\centering")
 x[16]<-paste0(x[16],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(clodur2)[1],2),"/",round(r.squaredGLMM(clodur2)[2],2),
               "&", round(r.squaredGLMM(platdur2)[1],2),"/", round(r.squaredGLMM(platdur2)[2],2),
               "&", round(r.squaredGLMM(closdis2)[1],2), "/", round(r.squaredGLMM(closdis2)[2],2),
               "&", round(r.squaredGLMM(closvel2)[1],2), "/", round(r.squaredGLMM(closvel2)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")

#' ## Table 7
 x<-capture.output(
   texreg::texreg(list(texreg::extract(di.clodur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(di.platdur2, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(di.closdis1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F),
                       texreg::extract(di.closvel1, include.variance=F,include.nobs=T,include.bic=F,include.loglik=F)
   ),
   override.pvalues = list(summary(di.clodur2)$coefficients[1:3,5],
                             summary(di.platdur2)$coefficients[1:3,5],
                             summary(di.closdis1)$coefficients[1:3,5],
                             summary(di.closvel1)$coefficients[1:3,5]),
   digits = 1,booktabs = T,custom.model.names = c("ClosDur","PlatDur","ClosDisp","ClosVel"),
   float.pos = "!htbp", label = "tab:articresults.di", use.packages = F,single.row = T,caption.above = T,
   caption="LME model predictions for the articulatory measurements of closing gesture duration  (ClosDur),    plateau duration (PlatDur), closing gesture displacement (ClosDisp), and closing gesture velocity (ClosVel)    of the last consonant in all disyllabic word pairs (standard error in parentheses)    for the fixed effects  phrasal position (phrase-medial vs. phrase-final)    and tenseness (tense vs. lax).",
   dcolumn = T, fontsize = "footnotesize",center=F)
 )
 x[3]<-paste0(x[3], "\\centering")
 x[16]<-paste0(x[16],"\\midrule R$_m^2$/R$_c^2$ & ", round(r.squaredGLMM(di.clodur2)[1],2),"/",round(r.squaredGLMM(di.clodur2)[2],2),
               "&", round(r.squaredGLMM(di.platdur2)[1],2),"/", round(r.squaredGLMM(di.platdur2)[2],2),
               "&", round(r.squaredGLMM(di.closdis1)[1],2), "/", round(r.squaredGLMM(di.closdis1)[2],2),
               "&", round(r.squaredGLMM(di.closvel1)[1],2), "/", round(r.squaredGLMM(di.closvel1)[2],2),
               "\\\\ \n")
 cat(x, sep = "\n")




 # plot articulatory data in boxplots from data subsets used for the models above
 # mono.art.sub1 # clodur mono
 # mono.art.plat1 # platdur mono
 # mono.art # clodis mono
 # mono.art.vel1 # clovel mono
 # di.art.sub1 # clodur di
 # di.art.plat1 # platdur di
 # di.art # clodis di
 # di.art # clovel di

 #' ## Figure 4
 #### position on x-axis, tenseness as fill
 m.clos<-ggplot(mono.art.sub1, aes(x=position, y=closdur,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 m.plat<-ggplot(mono.art.plat1, aes(x=position, y=platdur,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Plateau duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")


 m.disp<-ggplot(mono.art, aes(x=position, y=closdisp,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing displacement (mm)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 m.vel<-ggplot(mono.art.vel1, aes(x=position, y=PVEL..V.cm.sec.,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing velocity (cm/s)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 #### tenseness on x-axis, position as fill
 m.clos<-ggplot(mono.art.sub1, aes(x=tenseness, y=closdur,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 m.plat<-ggplot(mono.art.plat1, aes(x=tenseness, y=platdur,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Plateau duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")


 m.disp<-ggplot(mono.art, aes(x=tenseness, y=closdisp,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = 0),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing displacement (mm)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 m.vel<-ggplot(mono.art.vel1, aes(x=tenseness, y=PVEL..V.cm.sec.,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .1),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing velocity (cm/s)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 ggpubr::ggarrange(m.clos,m.plat,m.disp,m.vel,
                   common.legend = T,
                   labels = c("a)","b)","c)","d)"),
                   font.label = list(size=11,face="plain"),align = "h",nrow=2,ncol=2)
 ggsave(file="mono-art-raw-tenseness-on-x.pdf", width = 250, height = 250, units = "mm")


 #' ## Figure 5
 # di.art.sub1 # clodur di
 # di.art.plat1 # platdur di
 # di.art # clodis di
 # di.art # clovel di

 d.clos<-ggplot(di.art.sub1, aes(x=position, y=closdur,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 d.plat<-ggplot(di.art.plat1, aes(x=position, y=platdur,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Plateau duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")


 d.disp<-ggplot(di.art, aes(x=position, y=closdisp,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing displacement (mm)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 d.vel<-ggplot(di.art, aes(x=position, y=PVEL..V.cm.sec.,fill=tenseness))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing velocity (cm/s)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 #### tenseness on x-axis, position as fill

 d.clos<-ggplot(di.art.sub1, aes(x=tenseness, y=closdur,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 d.plat<-ggplot(di.art.plat1, aes(x=tenseness, y=platdur,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Plateau duration (ms)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")


 d.disp<-ggplot(di.art, aes(x=tenseness, y=closdisp,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing displacement (mm)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")

 d.vel<-ggplot(di.art, aes(x=tenseness, y=PVEL..V.cm.sec.,fill=position))+
   geom_point(position=position_jitterdodge(dodge.width = .75,jitter.width = .3,jitter.height = .2),
              alpha=1,size=.1,col="black")+
   geom_boxplot(alpha=.3,outlier.shape = NA, lwd = .3, notch = T)+
   facet_grid(~wordpair)+
   scale_fill_brewer(name="",palette="Set1")+
   # scale_x_discrete(labels=c("medial","final"))+
   ylab("Closing velocity (cm/s)")+
   xlab("")+
   theme_bw()+
   theme(legend.position="top")



 ggpubr::ggarrange(d.clos,d.plat,d.disp,d.vel,
                   common.legend = T,
                   labels = c("a)","b)","c)","d)"),
                   font.label = list(size=11,face="plain"),align = "h",nrow=2,ncol=2)
 ggsave(file="di-art-raw-tenseness-on-x.pdf", width = 200, height = 200, units = "mm")


#' # Figure of proportional lengthening
 sd.mono<-mono %>%
   group_by(tenseness, position, word, sylpos) %>%
   summarise(sylpos.sd = sd(sylposdur)) %>%
   spread(position, sylpos.sd) %>%
   rename( "medial.sd" = `phrase-medial`, "final.sd" = `phrase-final`)
 a<-mono %>%
   group_by(tenseness, position, word, sylpos) %>%
   summarise(sylpos.mean = mean(sylposdur)) %>%
   spread(position, sylpos.mean) %>%
   rename( "medial" = `phrase-medial`, "final" = `phrase-final`) %>%
   mutate(prop.len = final/medial,prop.len = prop.len-1,
          wordpairs = recode_factor(word, "Bahn" = "Bahn/Bann", "Bann" = "Bahn/Bann",
                                    "Beet"="Beet/Bett","Bett"="Beet/Bett",
                                    "Ruhm"="Ruhm/Rum", "Rum" = "Ruhm/Rum",
                                    "Stiel" = "Stiel/still", "still" = "Stiel/still"),
          sylpos = recode_factor(sylpos, "C" = "Coda", "N" = "Nucleus", "O" = "Onset"),
          sylpos = factor(sylpos, levels=c("Onset","Nucleus","Coda"))) %>%
   ggplot(aes(x=word,y=prop.len,  shape=tenseness,col=sylpos))+
   geom_point()+
   geom_line(aes(group=interaction(wordpairs, sylpos)))+
   scale_y_continuous(labels = scales::percent,breaks = seq(0,2,.2), name = "Lengthening proportion")+
   xlab("Target word")+
   theme_bw()+
   theme(legend.position = "top",legend.title = element_blank(),legend.box = "vertical")
 b<-di %>%
   group_by(tenseness, position, word, sylpos) %>%
   summarise(sylpos.mean = mean(sylposdur)) %>%
   spread(position, sylpos.mean) %>%
   rename( "medial" = `phrase-medial`, "final" = `phrase-final`) %>%
   mutate(prop.len = final/medial,prop.len = prop.len-1,
          wordpairs = recode_factor(word, "Miete" = "Miete/Mitte", "Mitte" = "Miete/Mitte",
                                    "Huete"="Hüte/Hütte","Huette"="Hüte/Hütte"),
          sylpos = recode_factor(sylpos, "O1" = "1st onset", "N1" = "1st nucleus", "O2" = "2nd onset", "N2" = "2nd nucleus"),
          sylpos = factor(sylpos, levels=c("1st onset","1st nucleus","2nd onset","2nd nucleus"))) %>%
   ggplot(aes(x=word,y=prop.len,  shape=tenseness,col=sylpos))+
   geom_point()+
   geom_line(aes(group=interaction(wordpairs, sylpos)))+
   scale_y_continuous(labels = scales::percent,breaks = seq(0,2.6,.2), name = "")+
   xlab("Target word")+
   theme_bw()+
   theme(legend.position = "top",legend.title = element_blank(),legend.box = "vertical",legend.text.align = 1)+
   guides(shape = guide_legend(order = 2),col = guide_legend(order = 1))
 cowplot::plot_grid(a,b,  labels = c("a) Monosyllabic","b) Disyllabic"),label_fontface = "plain", label_size = 11,
                    rel_widths = c(1,1), nrow = 1)