---
title: "/nt, nd/"
author: "Jonathan Harrington"
output: 
  bookdown::html_document2:
    number_sections: TRUE
    toc: true
    theme: flatly
    highlight: pygments
editor_options: 
  markdown: 
    wrap: 72
---

<style>
div.col {background-color: #ffd5c0; border-radius: 5px; padding: 20px;}
body {font-size: 16pt;}
h1 {font-size: 24pt;}
h2 {font-size: 22pt;}
p.caption {font-size: 12pt; text-align: justify;}
code.sourceCode {font-size: 16pt;}
</style>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


Load libraries and data:
```{r, message=FALSE}
# for time_in_sequence() function
library(tidyverse)
library(gridExtra)
library(emuR)
library(lmerTest)
library(emmeans)
library(glmmTMB)
library(bbmle)
#Rda Files
pfad = "/Users/cunha/Documents/projects/MRI_E/language/files"
#emu_DB
pfademu = "/Users/cunha/Documents/projects/MRI_E/DB/version_oct2024"


# directory for saving figures
pfadfigs = "/Users/cunha/Documents/projects/MRI_E/language/figs"
pfadfigs2 = pfadfigs


load(file = file.path(pfad,       
"erc_mrinasals_concatsigs_prenasalraising_v4_rsnonan_prevowelVokalCodapostcoda_A_tdn__IEAVe__i_pbfszS_01.Rda"))
```


# Preliminaries

## Figures

- Fig 2. fig2schematic.pdf
- Fig 3. fig3_nt.pdf
- Fig. 4. fig4ntb.pdf
- Fig. 5. fig5nt.pdf
- Fig. 6. fig.asynch.paper.pdf
- Fig. 7. t_align2.pdf
- Fig. 8. tt_vel_box.pdf
- Fig. 9. fig1paperntrev.pdf
- Fig. 10. fig2nt.pdf
- Fig. 11. fig.ttipaper.pdf
- Fig. 12. figttdisplvel.pdf
- Fig. 13. figttduration.pdf
- Fig. 14. durplot.pdf

These figures are created in nt_prob5_rev.Rmd

- Fig. 15. figalignprob.pdf
- Fig. 16. peryin1.pdf
- Fig. 17. SPL1.pdf
- Fig. 18. tmin.pdf
- Fig. 19. peryinfig.pdf
- Fig. 20. figscheme_new.pdf
- Fig. 21. model_paper.pdf

## keep only monosyllables
```{r}
df_erc_prenasalraising = df_erc_prenasalraising %>%
  filter(secvow == "_")
```

## Create a segment identifier
This is just speaker and trial number pasted together
```{r}
df_erc_prenasalraising = df_erc_prenasalraising  %>%
  mutate(rowind = paste(speaker, trialnumber, sep="."))
```
## Change speaker id to a factor

```{r}
df_erc_prenasalraising = df_erc_prenasalraising  %>%
  mutate(speaker = factor(paste0("S", speaker)))
```

## remove S41
```{r}
df_erc_prenasalraising = df_erc_prenasalraising  %>%
  filter(speaker!="S41") %>%
  mutate(speaker = factor(speaker))
```

## Identify nasality and change vowel labels


```{r, message=FALSE}
dt.df = df_erc_prenasalraising %>%
  mutate(Nasality = 
           case_when(substring(coda, 1, 1) %in% c("m", "n", "N") ~ "n", 
                              TRUE ~ "o")) %>%
  # change 'vowel labels
  mutate(vowel = case_when(vowel == "{" ~ "æ",
                           vowel == "E" ~ "ɛ",
                           vowel == "V" ~ "ʌ",
                            vowel == "I" ~ "ɪ",
                           vowel == "ei" ~ "eɪ")) %>%
  
# change the factor order for vowels
  mutate(vowel = factor(vowel, 
                        levels = 
                          c("æ", "eɪ", "ʌ", "ɛ", "ɪ")))

table(dt.df$Nasality)
table(dt.df$vowel)
```

## give orthographic labels
```{r}
# read in word table
words = read.table(file.path(pfad, "codeprosortho_list_corx.txt"))
# strip off final two characters in dt.df$labs
n = nchar(as.character(dt.df$labs))
w = substring(dt.df$labs, 1, n-2)
# match this to the first column of words
m = match(w, words[,1])
# check there are no NAs
any(is.na(m))
# these are the word labels
wfin = words[m,2]
# bind them in to dt.df
dt.df$orth = factor(wfin)

# divide USE speakers in 3 and 4 regions 
table(dt.df$speaker)
# 3 US regions  - US_northeast, US_south, US_west
dt.df= dt.df %>%
  mutate(region3 = recode (speaker, "S1" = "BR", "S10" = "BR","S11" = "BR","S12" = "BR","S13" = "BR",
        "S14" = "US_northeast","S15" = "BR","S16" = "BR","S17" = "BR","S18" = "BR",
        "S2" = "US_south","S20"="BR", "S21" = "BR","S22" = "BR","S23" = "BR","S24" = "US_northeast",
        "S25" = "US_west", "S26" = "BR", "S27" = "US_west", "S28" = "US_south","S29" = "BR",
        "S3" = "US_south","S30" = "US_northeast","S31" = "US_northeast","S32" = "BR","S33" = "BR",
        "S34" = "US_northeast","S35" = "US_northeast",
        "S36" = "BR","S37" = "BR","S38" = "BR","S39" = "BR","S4" = "US_west","S5" = "US_west",
        "S6" = "US_northeast","S7" = "US_west", "S8" = "BR","S9" = "BR",
        "S40" = "BR", "S42" = "US_west","S43" = "BR", "S44" = "BR","S45" = "BR"))
# split US_northeast (7/16  speakers)
# in US_Atlantic and Midland to obtain roughtly the same amount of speakers per group
dt.df= dt.df %>%
    mutate(region4 = recode (speaker, "S1" = "BR", "S10" = "BR","S11" = "BR","S12" = "BR","S13" = "BR",
        "S14" = "US_Atlantic","S15" = "BR","S16" = "BR","S17" = "BR","S18" = "BR",
        "S2" = "US_south","S20"="BR", "S21" = "BR","S22" = "BR","S23" = "BR","S24" = "US_midland",
        "S25" = "US_west", "S26" = "BR", "S27" = "US_west", "S28" = "US_south","S29" = "BR",
        "S3" = "US_south","S30" = "US_Atlantic","S31" = "US_midland","S32" = "BR","S33" = "BR",
        "S34" = "US_Atlantic","S35" = "US_midland",
        "S36" = "BR","S37" = "BR","S38" = "BR","S39" = "BR","S4" = "US_west","S5" = "US_west",
        "S6" = "US_midland","S7" = "US_west", "S8" = "BR","S9" = "BR",
        "S40" = "BR", "S42" = "US_west","S43" = "BR", "S44" = "BR","S45" = "BR"))
   table(dt.df$region3)
   table(dt.df$region4)
```

## Read in kinematic markers for nasalisation
```{r}
load(file.path(pfad, 
"erc_mrinasals_allsubjects_artaku_allsegandsig_table_v2_rsnonan_addlang_addtt_sig.Rda"))
dim(df_erc_prenasalraising)

# code for segment identification, as before
# code for segment identification, as before

df_erc_prenasalraising = df_erc_prenasalraising %>%
  mutate(rowind = paste(speaker, trialnumber, sep="."))


df_erc = df_erc_prenasalraising %>%
  select(velumopening_gesture_on, 
         velumopening_gesture_off,
         velumopening_maxvel_on,
         velumopening_maxvel_off,
         velumopening_maxcon_on,
         velumopening_nucleus_on,
         velumopening_nucleus_off,
         velumopening_nucleus_dur,
         velum2USV_velumopening_maxvel_onset,
         velum2USV_velumopening_maxvel_offset,
         velum2US_velumopening_maxcon_onset,
         velum2US_velumopening_nucleus_midpoint,
         velum2US_velumopening_gesture_onset,
         Vokal_on,
         rowind
         ) 

# join to dt.df
dt.df = left_join(dt.df, df_erc, 
             by = "rowind")
rm(df_erc, df_erc_prenasalraising)

# obtain time series offset
timeseriesoffset = dt.df %>% 
group_by (rowind) %>%
filter(segment == "Vokal") %>%
slice_head(n = 1) %>%
select(time_in_sequence, rowind) %>%
  rename(timeseriesoffset = time_in_sequence) %>%
ungroup()

# join to dt.df
dt.df = left_join(dt.df, timeseriesoffset, by = "rowind")

# obtain lineuptime and subtract
dt.df = dt.df %>%
  mutate(
  lineuptime = Vokal_on - timeseriesoffset, 
  velumopening_gesture_on = 
           velumopening_gesture_on - lineuptime,
         velumopening_gesture_off = 
           velumopening_gesture_off - lineuptime,
         velumopening_maxvel_on = 
           velumopening_maxvel_on - lineuptime,
         velumopening_maxvel_off = 
           velumopening_maxvel_off - lineuptime,
         velumopening_maxcon_on = 
           velumopening_maxcon_on - lineuptime,
         velumopening_nucleus_on = 
           velumopening_nucleus_on - lineuptime,
         velumopening_nucleus_off = 
           velumopening_nucleus_off - lineuptime)

```

## Read in kinematic markers for tongue tip

```{r}
load(file.path(pfad, 
"erc_mrinasals_allsubjects_artaku_allsegandsig_table_v2_rsnonan_addlang_addtt_sig.Rda"))
dim(df_erc_prenasalraising)
```

```{r}
# code for segment identification, as before
df_erc_prenasalraising = df_erc_prenasalraising %>%
  mutate(rowind = paste(speaker, trialnumber, sep="."))

# code for segment identification, as before
df_erc_prenasalraising = df_erc_prenasalraising %>%
  mutate(rowind = paste(speaker, trialnumber, sep="."))
df_erc = df_erc_prenasalraising %>%
  select(alveolarconstriction_maxcon_on, 
alveolarconstriction_nucleus_on,
alveolarconstriction_nucleus_off,
alveolarconstriction_maxvel_on,
alveolarconstriction_maxvel_off,
alveolarconstriction_nucleus_dur,
alveolarconstriction_maxvel_dur,
# displacements
alvUS_alveolarconstriction_maxvel_onset,
alvUS_alveolarconstriction_maxvel_offset,
alvUS_alveolarconstriction_maxcon_onset,
alvUS_alveolarconstriction_nucleus_onset,
alvUS_alveolarconstriction_nucleus_offset,
# position at gesture onset
alvUS_alveolarconstriction_gesture_onset,
# position at vowel minimum
alvUS_Vokal_min,
# velocities
alvUSV_alveolarconstriction_maxvel_onset,
alvUSV_alveolarconstriction_maxvel_offset,
Vokal_on,
rowind) 

# left join to dt.df
dt.df = left_join(dt.df, df_erc, 
             by = "rowind")
rm(df_erc, df_erc_prenasalraising)

# subtract lineuptime

dt.df = dt.df %>%
  mutate(alveolarconstriction_maxcon_on =
           alveolarconstriction_maxcon_on - lineuptime,
         alveolarconstriction_nucleus_off =
           alveolarconstriction_nucleus_off - lineuptime,
         alveolarconstriction_nucleus_on = 
           alveolarconstriction_nucleus_on - lineuptime,
         alveolarconstriction_maxvel_on =
           alveolarconstriction_maxvel_on - lineuptime,
         alveolarconstriction_maxvel_off =
           alveolarconstriction_maxvel_off - lineuptime) 


```

## Acoustic times

```{r}
times.df = dt.df %>%
  group_by(rowind) %>%
  filter(segment == "Vokal") %>%
  slice_tail(n=1) %>%
  mutate(endv_time = time_in_sequence,
         midv_time = (endv_time+timeseriesoffset)/2) %>%
  select(rowind, midv_time, endv_time) %>%
  ungroup()
# now join this to dt.df
dt.df = left_join(dt.df, times.df, 
             by = "rowind")
```


(Amended 21.06.2024). This is to include a marker for the beginning of acoustic vowel periodicity that is included in the data-frame `periodicV.df.txt`, derived as follows:


`autobuild_linkFromTimes(phys_db, superlevelName = "Wort", sublevelName = "Vokal")
summary(phys_db)`


`autobuild_linkFromTimes(phys_db, superlevelName = "Wort", sublevelName = "vowelfor")
add_linkDefinition(phys_db, "ONE_TO_MANY", "Wort", "vowelfor")`


```{r}
##Emu database
phys_db = load_emuDB(file.path(pfademu, 
                               "ercnasals_audio_emuDB"))
## serve(phys_db, useViewer=F)
## The following functions were carried out
## in order to link tiers:

##autobuild_linkFromTimes(phys_db, 
##                        superlevelName = "Wort", 
##                       sublevelName = "Vokal")

##autobuild_linkFromTimes(phys_db, 
##                        superlevelName = "Wort", 
##                        sublevelName = "vowelfor")
##add_linkDefinition(phys_db, "ONE_TO_MANY", 
##                   "Wort", "vowelfor")

## All segments of Vokal tier
vok.s = query(phys_db, "Vokal =~.*")
## as above without empty "" labels
vok.s = vok.s[vok.s$labels!="",]
## All segments from vowelfor tier
v.s = query(phys_db, "vowelfor =~.*")
## as above without empty "" labels
v.s = v.s[v.s$labels!="",]
## Test matching bundles
all(vok.s$bundle == v.s$bundle)
## get speaker-id
speaker = substring(v.s$bundle, 20, 22)
speaker = substring(speaker, 2, nchar(speaker))
speaker = as.numeric(speaker)
## get Trial-id
trial = as.numeric(substring(v.s$bundle, 
                             nchar(v.s$bundle)-3, 
                           nchar(v.s$bundle)))
##paste speaker.id and trial.id together
## to turn it into a character(!). This is needed
## so that a character vector is stored and then
## a character vector is read back in
id = paste("T", 
          paste(as.character(speaker), 
               as.character(trial), sep="."), sep=".")

## build the data frame
#df = data.frame(Vstart = vok.s$start/1000, 
#                Vend = vok.s$end/1000, 
#                V2start = v.s$start/1000, 
#                V2end = v.s$end/1000, 
#                rowind = id)
## write out the data-frame
#write.table(df, file = file.path(pfad, "periodicV.df.txt"), quote=F)
```

`df` has 4 columns that were obtained with the Emu script from the database `ercnasals_audio_emuDB` containing acoustic segmentations right at the end of this document. The four columns in `periodicV.df.txt` are:

- `Vstart`: start of the vowel as defined in the database tier `Vokal`
- `Vend`: end fo the vowel as defined in the database tier `Vokal`
- `Vstart2`: acoustic onset of periodicity of the vowel extracted from the database tier `vowelfor`
- `Vend2`: end fo the vowel as defined in the database tier `vowelfor` = `Vstart`

```{r}
df = read.table(file.path(pfad, "periodicV.df.txt"))
# make rowind match with that of dt.df
df$rowind = as.character(substring(df$rowind, 3, 
                                   nchar(df$rowind)))
# join to dt.df
dt.df = left_join(dt.df, df, by="rowind")
```
Verify that `Vend - Vstart` is the same as `endv_time - timeseriesoffset`. They are except for a 2 ms difference on average (i.e. a rounding error). 

```{r}
# dt.df %>%
# mutate(v1 = Vend - Vstart, 
# v2 = endv_time - timeseriesoffset, 
#d = v1 - v2) %>%
# summarise(mean(d), sd(d))
```
Now define the onset of periodicity in relation to `timeseriesoffset`. To do this, subtract `Vstart` from `V2start` then add that to `timeseriesoffset`

```{r}
dt.df = dt.df %>% 
  # this is the duration from the start of the vowel to the the beginning of periodicity
  mutate(d = V2start - Vstart, 
         # add this duration to `timeseriesoffset` which is also the start of the vowel but on a different time scale
         Vperiodiconset = timeseriesoffset + d)
```

Verify for any segment that `Vperiodiconset` corresponds to the start of the vowel in the `vowelfor` tier. This is done for segment `1.51` which is a production of 'pant'. This says that the duration from the start of the vowel to the onset of periodicity is 0.1065625 s.

```{r}
#speaker.i = 1
#trial.i = 51
#i = paste(speaker.i, trial.i, sep=".")
#dt.df %>%
#  filter(rowind %in% i) %>%
#  mutate(d = Vperiodiconset - timeseriesoffset) %>% 
#  pull(d) %>% unique()
```
Get the corresponding bundle from the Emu database. This matches the above time exactly. 
```{r}
#temp = (speaker == speaker.i) & (trial == trial.i)
#(v.s$start[temp] - vok.s$start[temp])/1000
```

Add an additional column, `Vowelonset2` which has the same start times as `timeseriesoffset` for /p, b/ initial words, otherwise `Vperiodiconset`

```{r}
dt.df = dt.df %>%
  mutate(Vowelonset2 = 
    ifelse(onset %in% c("p", "b"), 
           timeseriesoffset, Vperiodiconset))
```


Are there any non-values in the times of velum gesture? Yes, 10 of them:
```{r}
  gaps = dt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  filter(is.na(velumopening_maxcon_on) |
           is.na(velumopening_maxvel_off) | 
           is.na(velumopening_maxvel_on)) %>%
  pull(rowind)
gaps

# remove these segments from further consideration
dt.df = dt.df %>%
  filter(!rowind %in% gaps) 
```

Add CV variable

```{r}
dt.df = dt.df %>%
  mutate(CV = interaction(onset, vowel))
```

## data-frame of tongue tip data

```{r}
tt.rowind = dt.df %>%
  filter(Nasality == "n") %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  ungroup() %>%
  # choose ones that have values for the kinematic events
  filter(!is.na(alveolarconstriction_nucleus_on)) %>%
  filter(!is.na(alveolarconstriction_nucleus_off)) %>%
  pull(rowind)

tt.df = dt.df %>%
  filter(rowind %in% tt.rowind)

# number of segments: 1337
tt.df %>% pull(rowind) %>% n_distinct()

tt.df$orth = factor(tt.df$orth)
tt.df$onset = factor(tt.df$onset)
tt.df$vowel= factor(tt.df$vowel)
tt.df$coda = factor(tt.df$coda)
```

```{r}
tt.df %>%
   filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  # remove 2 NA
  filter(!(is.na(velum2USV_velumopening_maxvel_offset))) %>%
  rename(Dialect = language) %>% nrow()
  
#write.table(dt.df, file.path(pfad, 
#                             "dtphysnt.df.txt"), 
#            quote=F)
#write.table(tt.df, file.path(pfad, 
#                             "ttphysnt.df.txt"), 
#            quote=F)
```
A count of the data

```{r}
dt.df %>%
    filter(coda %in% c("d", "t", "nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1)   %>%
  ungroup() %>%
  group_by(coda, language) %>%
summarise(count = n()) %>%
  ungroup() %>%
  pivot_wider(names_from=coda, values_from=count) %>%
  relocate(language, nt, nd, t, d)
```


# Themes
```{r}
theme_michigan1 <-
function(n = 24){
  theme(
axis.text = element_text(size=n), 
axis.title.x = element_text(size=n), 
axis.title.y = element_text(size=n), 
text = element_text(size=n))
}
theme_phon = 
function(atextx = 18, atexty = 18, atitlex=
                 18, atitley = 18, ptitle=24, strip.text.x = 24,
             lpos = "bottom", othertext=24, ltitle =24){
        # e.g. atexty = element_blank()
        theme(
            axis.text.x = element_text(size=atextx),
            axis.title.x = element_text(size=atitley),
            axis.text.y = element_text(size=atexty),
            axis.title.y = element_text(size=atitley),
            plot.title = element_text(size = ptitle),
            legend.position=lpos,
            legend.title = ltitle,
            strip.text.x =element_text(size=strip.text.x),
            text = element_text(size=othertext))
    }
theme_michigan2 <-
function(n = 24){
    theme(
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        legend.title = element_blank(),
        axis.text.y = element_text(size=24),
        plot.title = element_text(size = 32),
        legend.position="top",
        text = element_text(size=n))
}
theme_michigan2epaper <-
function(n = 32, k = 28, position="bottom"){
    theme(
        axis.text.x = element_text(size=k),
        axis.title.x = element_text(size=k),
        axis.text.y = element_text(size=k),
axis.title.y = element_text(size=k),
        strip.text.x = element_text(size = n),
legend.title = element_text(size = n),
        legend.position=position,
        text = element_text(size=n))
}
theme_michigan2x <-
function(n = 24){
    theme(
text = element_text(size=n),
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        legend.title = element_blank(),
        axis.text.y = element_text(size=20),
axis.title.y = element_text(size=n),
        plot.title = element_text(size = 28),
        legend.position = "none")
}
theme_michigan2dpaper <-
function(n = 24, k = 18){
    theme(
        axis.text.x = element_text(size=k),
        axis.title.x = element_text(size=k),
        axis.text.y = element_text(size=k),
axis.title.y = element_text(size=k),
        plot.title = element_text(size = n),
        legend.position="bottom",
        text = element_text(size=n))
}
theme_michigan2apaper <-
function(n = 24, legend.pos="none", k = 18){
    theme(
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size=k),
axis.title.y = element_text(size=k),
        plot.title = element_text(size = n),
        legend.position=legend.pos,
        text = element_text(size=n))
}
theme_michigan2fpaper <-
function(n = 24, k = 18){
    theme(
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size=k),
axis.title.y = element_text(size=k),
        plot.title = element_text(size = n))
}
theme_michigan2gpaper <-
function(n = 24, k = 18){
    theme(
        axis.text.x = element_text(size=k),
        axis.title.x = element_text(size=k),
        axis.text.y = element_text(size=k),
axis.title.y = element_text(size=k),
        plot.title = element_text(size = n))
}

theme_michigan3 <-
function(n = 24){
    theme(
        legend.position = "none",
        axis.text.y = element_text(size=20),
axis.title.y = element_text(size=n),
        plot.title = element_text(size = 28), 
 text = element_text(size=n))
}
theme_michigan5 = 
function(n = 18, m = 32, i = 20){
    theme(    
        axis.text = element_text(size=n),
	axis.title = element_text(size=i),

        plot.title = element_text(size = m),
        legend.position="top",
        text = element_text(size=n))
}
theme_michigan8 <-
function(n = 48){
    theme(
        strip.text.x = element_text(size = 40),
        legend.title = element_text(size=40),
        legend.text = element_text(size=40),
        axis.title.x = element_text(size=36),
        axis.title.y = element_text(size=44),
        axis.text.x = element_text(size=32),
        axis.text.y = element_text(size=32),
        legend.position="top",
        text = element_text(size=n))
}
theme_michigan9 <-
function(n = 36){
    theme(
        axis.title.x = element_text(size=n),
        axis.title.y = element_text(size=n),
        axis.text.x = element_text(size=26),
        axis.text.y = element_text(size=26),
        legend.position="top",
        text = element_text(size=n))
}
theme_michigan10 <-
function(n = 36){
  theme(
axis.text = element_text(size=n), 
axis.title.x = element_text(size=n), 
axis.title.y = element_text(size=n), 
legend.position="top",
text = element_text(size=n))
}
theme_michigan10t = 
function(n = 36){
  theme(
axis.text = element_text(size=n), 
axis.title.x = element_text(size=n), 
axis.title.y = element_text(size=n), 
legend.position="top",
legend.title = element_blank(),
text = element_text(size=n))
}
theme_michigan11t = 
function(n = 32, k = 28, position="bottom"){
    theme(
        axis.text.x = element_text(size=k),
        axis.title.x = element_text(size=n),
        axis.text.y = element_text(size=k),
axis.title.y = element_text(size=n),
        strip.text.x = element_text(size = n),
        legend.position=position,
legend.title = element_blank(),
        text = element_text(size=n))
}

theme_michigan2hpaper =
function(n = 32, k = 28, position="bottom"){
        theme(
            axis.text = element_blank(),
            axis.title = element_text(size=k),
            legend.title = element_text(size = n),
            legend.position=position,
            text = element_text(size=n))
}
theme_michigan2xpaper =
function(n = 32, k = 28, position="bottom"){
        theme(
            axis.text = element_blank(),
            axis.title = element_text(size=k),
            legend.position="none",
            text = element_text(size=n))
}
```

write out `dt.df`
```{r}
#write.table(dt.df, file.path(pfad, "dt.df.txt"), quote=F)
```

# Intra-gestural velum analysis
```{r}
cols.coda = c("orange", "darkgreen")
cols.dial = c("red", "blue")
```

## Method
```{r, fig.width=25/2.54, fig.height=20/2.54}
j = "24.63"

# get some kinematic parameters
maxcon = dt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxcon_on) %>%
  unique()
velon = dt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxvel_on) %>%
  unique()
veloff = dt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxvel_off) %>%
  unique()

# times of plateau of peak displacement
nucleuson = dt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_nucleus_on) %>%
  unique()
nucleusoff = dt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_nucleus_off) %>%
  unique()

# average over plateau
plateau = dt.df %>%
  filter(rowind == j) %>%
filter(time_in_sequence >= velumopening_nucleus_on) %>%
   filter(time_in_sequence <= velumopening_nucleus_off) %>%
   summarise(displ = mean(velum2US)) %>%
  pull(displ)

# the peak displacement value
pkdisp = dt.df %>%
  filter(rowind == j) %>%
  pull(velum2US) %>% max()

# get the peak velocity value
pkvel = dt.df %>%
  filter(rowind == j) %>%
   mutate(velocity = 500.5 * c(stats::filter(velum2US, c(.5, 0, -.5)))) %>%
  pull(velocity) %>% max(na.rm=T)

# get the peak velocity minimum
pkvel.raise = dt.df %>%
  filter(rowind == j) %>%
   mutate(velocity = 500.5 * c(stats::filter(velum2US, c(.5, 0, -.5)))) %>%
  pull(velocity) %>% min(na.rm=T)
# readjust slightly for plotting purposes
pkvel.raise = -7.6
j = "24.63"

d1paper = dt.df %>%
  filter(rowind == j) %>%
  ggplot +
  aes(y = velum2US, x = time_in_sequence) +
  geom_line(linewidth=1.5) +
#geom_segment(aes(y = velum2US, x = time_in_sequence), 
#               x = maxcon, y = plateau, 
#               xend = 0, yend = plateau, lty=2,
#               arrow=arrow(ends="last", type="closed")) +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = (nucleuson+nucleusoff)/2, y = plateau, 
               xend = (nucleuson+nucleusoff)/2, yend = 0, lwd=1.5,
               arrow=arrow(ends="both", type="closed")) +
  ylab("Displacement") +
  geom_vline(xintercept = c(nucleuson, nucleusoff), 
           lty=2, linewidth=1.5) +
  theme_michigan2apaper(k = 24) +
scale_colour_manual(values = cols)

v1paper = dt.df %>%
  filter(rowind == j) %>%
   mutate(velocity = 500.5 * c(stats::filter(velum2US, c(.5, 0, -.5)))) %>%
  ggplot +
  aes(y = velocity, x = time_in_sequence) +
  geom_line(linewidth=1.5)  +

  geom_vline(
    xintercept = c(velon, veloff), 
             lty=2, linewidth=1.5) +
  geom_segment(aes(y = velocity, x = time_in_sequence), 
               x = velon, y = pkvel, 
               xend = velon, yend = 0, lwd=1.5,
               arrow=arrow(ends="first", type="closed")) +
   geom_segment(aes(y = velocity, x = time_in_sequence), 
               x = veloff, y = pkvel.raise, 
               xend = veloff, yend = 0, lwd=1.5, 
               arrow=arrow(ends="first", type="closed")) +
geom_hline(yintercept = 0) +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = velon, y = 0, 
               xend = veloff, yend = 0, lty=2,
               arrow=arrow(ends="both", type="closed")) +
  ylab("Velocity (units/s)") +
  xlab("Time (s)") + 
  theme_michigan2dpaper(k = 24) +
scale_colour_manual(values = cols)

# we need this to align the x-axis of the two
# plots. Without it, v1paper is right shifted
# because of the extra space taken up
# by the y-axis tick-mark-labels
gA <- ggplotGrob(d1paper)
gB <- ggplotGrob(v1paper)

grid::grid.newpage()
fig3_nt = grid.arrange(rbind(gA, gB))
fig3_nt
ggsave(filename = file.path(pfadfigs, "fig3_nt.png"), 
    plot = fig3_nt, 
      width = 25, 
       height = 20, 
   units = "cm")

```

```{r}
## Aligned at acoustic vowel onset (all vowels)
z = dt.df %>%
  mutate(Vowelonset2 = 
    ifelse(onset %in% c("p", "b"), 
           timeseriesoffset, Vperiodiconset))

v.df = dt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  mutate(reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  
                                velumopening_maxcon_on), 5), 
         coda = factor(coda))

aconset.t = v.df %>%
  rename(Dialect = language) %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval2 = endv_time - velumopening_maxcon_on, 
         intval1 = Vowelonset2 - 
           velumopening_maxcon_on) %>%
  group_by(vowel, coda, Dialect) %>%
  summarise(m2 = mean(intval2, na.rm=T) * 1000, 
            m1 = mean(intval1, na.rm=T) * 1000) %>%  ungroup()
```

```{r, fig.width=45/2.54, fig.height=30/2.54}
# Aligned at peak velum lowering
v.df = dt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  mutate(reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  
                                velumopening_maxcon_on), 5), 
         coda = factor(coda))


aconset.t = v.df %>%
  rename(Dialect = language) %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval2 = endv_time - velumopening_maxcon_on, 
         intval1 = Vowelonset2 - 
           velumopening_maxcon_on) %>%
  group_by(vowel, coda, Dialect) %>%
  summarise(m2 = mean(intval2, na.rm=T) * 1000, 
            m1 = mean(intval1, na.rm=T) * 1000) %>%  ungroup()


n =48; k = 32; i = 28
lwd = 1.5
fig4ntb = v.df %>%
group_by(language, coda, vowel, reltime)  %>%
    summarise(velum = mean(velum2US)) %>%
  ungroup() %>%
  rename(Dialect = language) %>%
    ggplot + 
    aes(y = velum, x = reltime, col = coda, 
        group = interaction(Dialect, coda, vowel)) +
  geom_line(linewidth=lwd) +
  scale_x_continuous("Time (ms)", breaks = c(-200,  0, 100)) +
  coord_cartesian(xlim =c(-300, 150)) +
    ylab("Velum displacement") +
    facet_wrap(Dialect ~ vowel, ncol=5) +
  theme_phon(lpos="top",
             atextx=i,
             atexty=k,
             atitlex=k,
             atitley=k,
             strip.text.x = k,
             ptitle=n,
             ltitle=element_blank(),
             othertext=n) + 
    geom_vline(xintercept=0, lty=2, linewidth=lwd) +
   geom_vline(data = aconset.t, 
             aes(xintercept=m2,  col = coda),  
             lty=3,linewidth=lwd) +
  geom_vline(data = aconset.t, 
             aes(xintercept=m1,  col = coda),  
             lty=1,linewidth=.75) +
  scale_colour_manual(values = cols.coda) 
fig4ntb
ggsave(filename = file.path(pfadfigs2, "fig4ntb.pdf"), 
      plot = fig4ntb, 
      device=cairo_pdf,
       width = 45, 
       height = 30, 
      units = "cm")
```






## Displacement

Displacement is averaged over the plateau velum movement
```{r}

displ.df = dt.df %>% 
  filter(!is.na(velum2US_velumopening_gesture_onset))  %>%
  mutate(velum2US = 
           velum2US - velum2US_velumopening_gesture_onset) %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(time_in_sequence >= velumopening_nucleus_on) %>%
   filter(time_in_sequence <= velumopening_nucleus_off) %>%
  group_by(rowind, language, 
           coda, CV, speaker) %>%
   summarise(displ = mean(velum2US)) %>%
  rename(Dialect = language) %>%
  ungroup() 
```


```{r}
displ.df %>%
     lmer(displ ~  coda  * Dialect +
            (Dialect|CV)  +
            (coda | speaker), .) %>%
     step()
```
 
 The model found was `displ ~ coda + (Dialect | CV) + (coda | speaker)`. Further simplification of RF so that it converges

```{r}
displ.df %>%
     lmer(displ ~ coda + (Dialect | CV) + (coda | speaker), .) %>%
     anova()
```
Estimated marginal means
```{r}
displ.df %>%
     lmer(displ ~ coda + (Dialect | CV) + (coda | speaker), .) %>%
  emmeans(., ~coda)
```


## peak velocity of velum lowering and raising

### figure

Calculate velocity from displacement signal for plotting purposes only. 
```{r}
dt.df = dt.df %>%
  group_by(rowind) %>%
  # mutliplied by the sample rate
   mutate(velum2USV = 500.5 * 
            c(stats::filter(velum2US, c(.5, 0, -.5))))  %>%
  ungroup()
```

Aggregate plot of velum velocity

```{r, fig.width=45/2.54, fig.height=30/2.54}
n =48; k = 32; i = 28
lwd = 1.5
fig5nt = dt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  mutate(reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  
                                velumopening_maxcon_on), 5), 
         coda = factor(coda)) %>%
group_by(language, coda, vowel, reltime)  %>%
    summarise(velum = mean(velum2USV, na.rm=T)) %>%
  ungroup() %>%
  rename(Dialect = language) %>%
    ggplot + 
    aes(y = velum, x = reltime, col = coda, 
        group = interaction(Dialect, coda, vowel)) +
  geom_line(linewidth=lwd) +
  scale_x_continuous("Time (ms)", breaks = c(-100, 0, 100)) +
 coord_cartesian(xlim =c(-200, 200)) +
    ylab("Velocity") +
    facet_wrap(Dialect ~ vowel, ncol=5) +
  theme_phon(lpos="top",
             atextx=i,
             atexty=k,
             atitlex=k,
             atitley=k,
             strip.text.x = k,
             ptitle=n,
             ltitle=element_blank(),
             othertext=n)  +
  scale_colour_manual(values = cols.coda) 
fig5nt
ggsave(filename = file.path(pfadfigs2, "fig5nt.pdf"), 
      plot = fig5nt, 
      device=cairo_pdf,
       width = 45, 
       height = 30, 
      units = "cm")
```


### statistics
#### lowering
```{r}
pklower.df = dt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  # remove 2 NA
  filter(!(is.na(velum2USV_velumopening_maxvel_onset))) %>%
  rename(Dialect = language)
pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ 
       coda * Dialect +
       (Dialect|CV) +
     ( coda | speaker), .) %>%
step()

# model found
pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ coda + (1 | CV) + (1 | speaker), .) %>%
anova()
# EMM
pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ coda + (1 | CV) + (1 | speaker), .) %>%
  emmeans(., ~coda)


```

Get the means for coda and vowel. For coda: /nd/ = 6.17; /nt/ = 6.63.



#### raising
```{r}
pkraise.df = dt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  # remove 2 NA
  filter(!(is.na(velum2USV_velumopening_maxvel_offset))) %>%
  rename(Dialect = language)
```



Model found: `velum2USV_velumopening_maxvel_offset ~ (Dialect | CV) + (coda | speaker)`
```{r}
pkraise.df %>%
lmer(velum2USV_velumopening_maxvel_offset ~ 
       coda *  Dialect +
       (Dialect|CV) +
     ( coda | speaker), .) %>%
step()

```




##  articulatory duration
```{r}
artdur.df = dt.df %>%
  group_by(rowind) %>%
  slice_head(n = 1) %>%
  ungroup() %>%
  filter(coda %in% c("nt", "nd")) %>%
  mutate(art_time =
           velumopening_maxvel_off - 
           velumopening_maxvel_on
           ) %>%
  rename(Dialect = language)

artdur.df %>%
  lmer(art_time ~ coda * Dialect +
        (Dialect|CV) +
     (coda|speaker), .) %>%
step()

# model found

artdur.df %>%
  lmer(art_time ~ coda + (Dialect | CV) 
       + (coda | speaker), .) %>%
 anova()

# EMM
artdur.df %>%
  lmer(art_time ~ coda + (Dialect | CV) 
       + (coda | speaker), .) %>%
 emmeans(., ~coda)
```


# Inter-gestural analysis of tongue tip and velum 
## method (fig.asynch.paper.png)
Word is `band`, BRE speaker.


```{r}
j = "1.239"
# velum kinematic parameters
maxcon.vel = tt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxcon_on) %>%
  unique()
velon.vel = tt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxvel_on) %>%
  unique()
veloff.vel = tt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxvel_off) %>%
  unique()
plateau.on.tt = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_nucleus_on) %>%
  unique()
plateau.off.tt = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_nucleus_off) %>%
  unique()
ttmax.t = (plateau.on.tt + plateau.off.tt)/2
# TT raising pk velocity max
velon.tt = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_maxvel_on) %>%
  unique()
# set the time scale to be the same for
# TT and velum
 xlim.trade  = tt.df %>%
  filter(rowind == j) %>%
  pull(time_in_sequence) %>%
  range()
```

```{r, fig.width=45/2.54, fig.height=30/2.54}

n = 54; k = 40
yendline = 1
velum_asynch_paper = tt.df %>%
  filter(rowind == j) %>%
  filter(time_in_sequence >= timeseriesoffset) %>%
  filter(time_in_sequence <= velumopening_maxvel_off) %>%
  ggplot +
  aes(y = velum2US, x = time_in_sequence) +
  geom_line(linewidth = lwd) +
  coord_cartesian(xlim = c(0.1, veloff.vel), ylim = c(0.18, 1.12)) +
  scale_y_continuous("Displacement", breaks = c(.4, .6, .8)) +
 
ggtitle("Velum") +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = maxcon.vel, y = .1, 
               xend = maxcon.vel, yend = yendline, lty=2, lwd=lwd, col="black") +
   geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = velon.tt, y = .1, 
               xend = velon.tt, yend = yendline, lty=2, lwd=lwd, col="slategray") +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = ttmax.t, y = .1, 
               xend = ttmax.t, yend = yendline, lty=2, lwd=lwd, col="slategray") +
  
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = velon.tt, y = .5, 
               xend = maxcon.vel, yend = .5, col="black",
               arrow=arrow(ends="both", type="closed")) +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = ttmax.t, y = .3, 
               xend = maxcon.vel, yend = .3, col="black",
               arrow=arrow(ends="both", type="closed")) +
   annotate("text", x=mean(c(velon.tt, maxcon.vel)), y=0.57, 
           label="italic(i[1])", 
           parse=TRUE, 
           size=n/.pt) +
  annotate("text", x=mean(c(maxcon.vel, ttmax.t)), y=0.37, 
           label="italic(i[2])", 
           parse=TRUE, 
           size=n/.pt) +
  annotate("text", x=c(velon.tt,  maxcon.vel, ttmax.t), 
           y = rep(1.09, 3),
           label=c("italic(t[1])", "italic(t[2])", "italic(t[3])"),
           parse=TRUE, 
           size=n/.pt) +
theme_michigan2fpaper(n= n, k = k)


# draw tongue tip with kinematic markers
tt_asynch_paper = tt.df %>%
  filter(rowind == j) %>%
  filter(time_in_sequence >= timeseriesoffset) %>%
   filter(time_in_sequence <= velumopening_maxvel_off) %>%
  ggplot +
  aes(y = alvUS, x = time_in_sequence) +
  geom_line(col="slategray",  linewidth = lwd) +
  coord_cartesian(xlim = c(0.1, veloff.vel)) +
  geom_vline(xintercept = c(velon.tt, ttmax.t), 
             lty=2, linewidth = lwd, 
             col = c("slategray", "slategray")) +
  xlab("Time (s)") +
   ylab("Displacement") + 
  ggtitle("Tongue tip") +
  theme_michigan2gpaper(n = n, k = k)

fig.asynch.paper = 
  grid.arrange(velum_asynch_paper, 
               tt_asynch_paper, nrow=2)
ggsave(filename = file.path(pfadfigs,
                            "fig.asynch.paper.pdf"), 
 plot = fig.asynch.paper, 
 device = cairo_pdf,
       width = 45, 
      height = 30, 
       units = "cm")
fig.asynch.paper
```

## Aggregated plot
Parameters:

```{r}

x2.df = tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  mutate(
    peakVE = (velumopening_nucleus_on +
      velumopening_nucleus_off)/2,
    peakTT = (alveolarconstriction_nucleus_on +
      alveolarconstriction_nucleus_off)/2,
    zerot = peakVE, 
    reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  zerot), 5), 
         coda = factor(coda))

tt2.maxr.t = x2.df %>%
  rename(Dialect = language) %>%
  group_by(vowel, coda, Dialect) %>%
  summarise(m1 = mean(alveolarconstriction_maxvel_on -  
                       zerot, na.rm=T) * 1000, 
            m2 = mean(peakTT -  
                       zerot, na.rm=T) * 1000, 
            m3 = mean(endv_time -  
                       zerot, na.rm=T) * 1000) %>%  
  ungroup()
```

Plot:
```{r, fig.width=40/2.54, fig.height=30/2.54}
t_align2 = x2.df %>%
group_by(language, coda, vowel, reltime)  %>%
    summarise(tt = mean(alvUS - alvUS_Vokal_min)) %>%
  ungroup() %>%
  rename(Dialect = language) %>%
    ggplot + 
    aes(y = tt, x = reltime, col = coda, 
        group = interaction(Dialect, coda, vowel)) +
  geom_line(linewidth=1.2) +
  scale_x_continuous("Time (ms)", breaks = c(-100, 0, 100)) +
    ylab("Tongue tip height") +
    facet_wrap(Dialect ~ vowel, ncol=3) +
  coord_cartesian(xlim = c(-150, 150)) +
    geom_vline(xintercept=0, lty=2) +
   theme_michigan2epaper(position="top") + 
   theme(legend.title = element_blank()) +
  geom_vline(data = tt2.maxr.t, 
             aes(xintercept=m2,  col = coda),  lty=2, 
             linewidth=1.5) + 
   geom_vline(data = tt2.maxr.t, aes(xintercept=m1,  
                                    col = coda), lty=2,
              linewidth=1.5) + 
  geom_vline(data = tt2.maxr.t, 
             aes(xintercept=m3,  col = coda),  lty=1, 
             linewidth=.5) + 
 #  ggtitle("Tongue tip height", 
 #          subtitle = "Aligned at peak velum opening") +
  scale_colour_manual(values = cols.coda) 

ggsave(filename = file.path(pfadfigs,
                            "t_align2.pdf"), 
     plot = t_align2, 
     device=cairo_pdf,
       width = 40, 
      height = 30, 
       units = "cm")
t_align2
```

Parameters i1, i2
```{r}
z.df = tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  mutate(
    peakVE = (velumopening_nucleus_on +
      velumopening_nucleus_off)/2,
    peakTT = (alveolarconstriction_nucleus_on +
      alveolarconstriction_nucleus_off)/2, 
    i1 = alveolarconstriction_maxvel_on - peakVE, 
    i2 = peakTT - peakVE, 
    i3 = endv_time - peakVE)
```


Boxplots showing `i1` and `i2`

```{r,fig.width = 30/2.54, fig.height=20/2.54}
lwd = 1.4
tt_vel_a = z.df %>%
  mutate(i1 = i1 * 1000) %>%
 ggplot + 
 aes(y = i1, col = coda, x = Dialect) +
     geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 3) +
  theme_michigan2(n=26) +
scale_y_continuous(expression(italic(i[1])~ms), 
                   breaks=c(-100, -50, 0)) +
 theme(axis.text.y = 
         element_text(angle = 90, size=20,
                      vjust = 1, hjust=1))   +
  xlab("") +
coord_cartesian(ylim = c(-175, 50)) + 
scale_colour_manual(values = cols.coda)

tt_vel_b = z.df %>%
  mutate(i2 = i2 * 1000) %>% 
 ggplot + 
 aes(y = i2, col = coda, x = Dialect) +
     geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 3) +
  theme_michigan3(n=26) +
theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
scale_y_continuous(expression(italic(i[2])~ms), 
                   breaks=c(0, 100, 200)) +
  xlab("") +
  coord_cartesian(ylim = c(-50, 200)) + 
scale_colour_manual(values = cols.coda)


tt_vel_box = 
  grid.arrange(tt_vel_a, tt_vel_b, nrow=2)
ggsave(filename = file.path(pfadfigs2,
                            "tt_vel_box.pdf"), 
      plot = tt_vel_box, 
      device=cairo_pdf,
       width = 30, 
       height = 20, 
      units = "cm")
```





## Statistical test on i1 and i2
i1: 

- BRE > USE in BOTH /nt, nd/ 
- nd > nt for BRE only

```{r}
z.df %>%
 lmer(i1   ~ 
         vowel * coda * Dialect + 
         (Dialect | onset) + 
         (vowel + coda| speaker), .) %>%
  step()

# model found
z.df %>%
 lmer(i1 ~ coda + Dialect + (coda | speaker) + 
        (1 | CV) + coda:Dialect, .) %>%
  anova()

z.df %>%
 lmer(i1 ~ coda + Dialect + (coda | speaker) + 
        (1 | CV) + coda:Dialect, .) %>%
emmeans(.,  ~ Dialect|coda) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
```

i2: 

- model found is `i2 ~ vowel + coda + (coda | speaker)`
- thus no dialect differences, only vowel and coda
- for coda, duration from velum peak to TT peak is greater for /nt/ than for /nd/.

```{r}
z.df %>%
 lmer(i2   ~ 
         coda * Dialect + 
         (Dialect | CV) + 
         (coda| speaker), .) %>%
  step()

z.df %>%
 lmer(i2 ~ coda + (coda | speaker) + (1 | CV), .) %>%
  anova()

# EMM
z.df %>%
 lmer(i2 ~ coda + (coda | speaker) + (1 | CV), .) %>%
  emmeans(., ~ coda)

```




# Nasalization relative to the acoustic VN boundary 
## Method

extract the acoustic waveform from the Emu database
```{r}
temp = speaker == 24 & trial == 63
segment.s = vok.s[temp,]
segment.s =  segment.s %>%
  mutate(start = start -200,
         end = end + 300)
segment.wav = get_trackdata(phys_db, segment.s, "MEDIAFILE_SAMPLES")
# time of periodic onset of the vowel in Emudatabase: 0.7880938 * 1000
acvowelon.t = 
  df %>% filter(rowind == "24.63") %>% pull(V2start) * 1000

# Change time scale in segment.wav to make the periodic onset in the
# vowel have a time of zero ms
segment.wav = segment.wav %>%
  mutate(times_orig  =times_orig -acvowelon.t)
```

Alignment of phys. data relative to the acoustic vowel onset
```{r}
lwd = 1.5
j = "24.63"
j.df =  dt.df %>%
  mutate(
    zerot = Vperiodiconset,  
    reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  zerot), 5)) %>%
   filter(rowind == j)


j.df =  dt.df %>%
  mutate(
    zerot = Vperiodiconset,  
    reltime =  1000 * (time_in_sequence -  zerot)) %>%
   filter(rowind == j)
```
Some time markers
```{r}
# acoustic vowel offset time
VNbound24.63 = j.df %>%
  mutate(endv_time = 1000 * (endv_time - zerot)) %>%
  pull(endv_time) %>% unique()

# time of peak velocity of velum lowering
pkvellower24.63 = 
j.df %>%
  mutate(velumopening_maxvel_on = 
           1000 * (velumopening_maxvel_on - zerot)) %>%
  pull(velumopening_maxvel_on) %>% unique()

# time of velum lowering peak
pkvel24.63 = 
j.df %>%
  mutate(velumopening_maxcon_on = 
           1000 * (velumopening_maxcon_on - zerot)) %>%
  pull(velumopening_maxcon_on) %>% unique()

# time of peak velocity of velum raising
pkvelraise24.63 = 
j.df %>%
  mutate(velumopening_maxvel_off = 
           1000 * (velumopening_maxvel_off - zerot)) %>%
  pull(velumopening_maxvel_off) %>% unique()
```

Set the scale on the x-axis
```{r}
timeaxis = c(-50, pkvelraise24.63 + 50)
```

Data frames for drawing the areas
```{r}
data1 = j.df %>%
  filter(reltime >= VNbound24.63) %>%
  filter(reltime <= pkvelraise24.63)

data2 = j.df %>%
  filter(reltime <= VNbound24.63) %>%
  filter(reltime >= pkvellower24.63)
```


```{r,fig.width = 30/2.54, fig.height=20/2.54}
ystart = -.1; yend24.63 = .75; yend24.63text = yend24.63 + .02
textsize = 24
 nasal_trajectory =  j.df %>%
  ggplot + 
  aes(y = velum2US - velum2US_velumopening_gesture_onset, 
      x = reltime) +
   coord_cartesian(xlim = timeaxis, ylim = c(-.05, .77)) + 
  geom_line() +
  xlab("Time (ms)") +
  ylab("Size of velum lowering") + 
  geom_area(
    aes(y = velum2US - velum2US_velumopening_gesture_onset, 
        x = reltime), 
        data = data1, 
    fill = 4,
            alpha = 0.5,
            color = 1,    
            lwd = 0.5,    
            linetype = 1) +
  geom_area(
    aes(y = velum2US - 
          velum2US_velumopening_gesture_onset, x = reltime), 
        data = data2, 
    fill = "orange",
            alpha = 0.5,
            color = 1,    
            lwd = 0.5,    
            linetype = 1) +
theme_michigan1(n = textsize) +
scale_colour_manual(values = cols) +
  # acoustic vowel offset
  geom_segment(aes(y = velum2US - 
          velum2US_velumopening_gesture_onset, x = reltime), 
             x = VNbound24.63, y = ystart, xend = VNbound24.63,
             yend = yend24.63, lty=3, linewidth=lwd) +
   annotate("text", x=VNbound24.63, y=yend24.63text, 
           label="italic(t[2])", 
           parse=TRUE, 
           size=textsize/.pt) +
   # peak velum lowering
   geom_segment(aes(y = velum2US - 
          velum2US_velumopening_gesture_onset, x = reltime), 
             x = pkvel24.63, y = ystart, xend = pkvel24.63,
             yend = yend24.63, lty=2, linewidth=lwd) +
    annotate("text", x=pkvel24.63, y=yend24.63text, 
           label="italic(t[3])", 
           parse=TRUE, 
           size=textsize/.pt) +
   # pk vel velum lowering
   geom_segment(aes(y = velum2US - 
          velum2US_velumopening_gesture_onset, x = reltime), 
             x = pkvellower24.63, y = ystart, xend = pkvellower24.63,
             yend = yend24.63, lty=2, linewidth=lwd) +
   annotate("text", x=pkvellower24.63, y=yend24.63text, 
           label="italic(t[1])", 
           parse=TRUE, 
           size=textsize/.pt) +
     # pk vel velum raising
      geom_segment(aes(y = velum2US - 
          velum2US_velumopening_gesture_onset, x = reltime), 
             x = pkvelraise24.63, y = ystart, xend = pkvelraise24.63,
             yend = yend24.63, lty=2, linewidth=lwd) +
   annotate("text", x=pkvelraise24.63, y=yend24.63text, 
           label="italic(t[4])", 
           parse=TRUE, 
           size=textsize/.pt) +
   # acoustic vowel onset
   geom_segment(aes(y = velum2US - 
          velum2US_velumopening_gesture_onset, x = reltime), 
             x = 0, y = ystart, xend = 0,
             yend = yend24.63, lty=3, linewidth=lwd) +
   annotate("text", x=0, y=yend24.63text, 
           label="italic(t[0])", 
           parse=TRUE, 
           size=textsize/.pt) 

waveform = segment.wav %>%
  ggplot +
  aes(y = T1, x = times_orig) +
  geom_line() +
   coord_cartesian(xlim = timeaxis, ylim=c(-12000, 12000)) +
  theme(
        axis.text = element_blank(),
        axis.title = element_blank(),
        legend.title = element_blank(),
        legend.position = "none") +
  # acoustic vowel onset and offset
  geom_vline(xintercept = c(0, VNbound24.63),  lty=3, linewidth=lwd) 

gA <- ggplotGrob(nasal_trajectory)
gB <- ggplotGrob(waveform)

grid::grid.newpage()
fig1paperntrev = grid.arrange(rbind(gA, gB))

fig1paperntrev
ggsave(filename = file.path(pfadfigs2, "fig1paperntrev.pdf"), 
       plot = fig1paperntrev, 
       device=cairo_pdf,
      width = 30, 
       height = 20, 
       units = "cm")
```

Sum of nasalisation in coda

```{r}
propc.df = dt.df %>%
  mutate(velum2US = 
           velum2US - velum2US_velumopening_gesture_onset) %>%
  filter(coda %in% c("nd", "nt")) %>%
    filter(time_in_sequence >= endv_time) %>%
  filter(time_in_sequence <= velumopening_maxvel_off) %>%
  #filter(time_in_sequence <= velumopening_nucleus_off) %>%
  group_by(rowind, vowel, coda, language, CV, speaker) %>%
  summarise(cprop = sum(velum2US)) %>%
  ungroup()
```

Sum of nasalisation between points of maximum velocity.
```{r}
nasdur.df = dt.df %>%
  mutate(velum2US = 
           velum2US - velum2US_velumopening_gesture_onset) %>%
  filter(coda %in% c("nd", "nt")) %>%
filter(time_in_sequence >= velumopening_maxvel_on) %>%
 # filter(time_in_sequence >= timeseriesoffset) %>%
 filter(time_in_sequence <= velumopening_maxvel_off) %>%
  #filter(time_in_sequence <= velumopening_nucleus_off) %>%
  group_by(rowind) %>%
  summarise(nasalprop = sum(velum2US)) %>%
  ungroup()

```

Ratio
```{r}
propcnas.df = left_join(propc.df, nasdur.df, 
             by = "rowind") %>%
  mutate(ratio = cprop/nasalprop)
```

## Results
### Plot
```{r,fig.width = 30/2.54, fig.height=20/2.54}
lwd = 1.4
fig2ant = dt.df %>%
  rename(Dialect = language) %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(vn2.num = 
           velumopening_maxcon_on - Vperiodiconset, 
         vn2.den = endv_time - Vperiodiconset, 
         vn2 = vn2.num/vn2.den) %>%
  ggplot + 
  aes(y = log(vn2), col = coda, x = Dialect) +
  geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 5) +
  theme_michigan2(n=26) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
 scale_y_continuous("Proportional offset", breaks=c(-0.5,  0, 0.5)) +
  xlab("") +
  # ggtitle("Proportional alignment of peak velum lowering") +
 coord_cartesian(ylim = c(-1,1.0)) + 
 geom_hline(yintercept=0, lty=2) +
scale_colour_manual(values = cols.coda)

fig2bnt = propcnas.df %>%
  filter(ratio <1) %>%
  rename(Dialect = language) %>% 
  ggplot +
  aes(y = ratio, x = Dialect, col=coda) +
  geom_boxplot(lwd=lwd) +
   geom_hline(yintercept=0.5, lty=2) + 
  facet_wrap(~ vowel, ncol = 5) +
  ylab("Ratio") +
  xlab("") +
theme_michigan3(n = 26) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
  scale_colour_manual(values = cols.coda) 

fig2nt = 
  grid.arrange(fig2ant, fig2bnt, nrow=2)
ggsave(filename = file.path(pfadfigs2, "fig2nt.pdf"), 
      plot = fig2nt, 
      device=cairo_pdf,
       width = 30, 
       height = 20, 
      units = "cm")
```

### Statistics row 1 (tprop)

```{r}

dt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  mutate(vn2.num = 
        (velumopening_nucleus_on + 
           velumopening_nucleus_off)/2 - Vperiodiconset,
       vn2.den = endv_time - Vperiodiconset,
         vn2 = vn2.num/vn2.den) %>%
  lmer(log(vn2) ~ Dialect * coda + 
         (coda | speaker) + (Dialect | CV), .) %>% step()

# model found
dt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  mutate(vn2.num = 
        (velumopening_nucleus_on + 
           velumopening_nucleus_off)/2 - Vperiodiconset,
       vn2.den = endv_time - Vperiodiconset,
         vn2 = vn2.num/vn2.den) %>%
  lmer(log(vn2) ~ Dialect + (coda | speaker) + 
         (Dialect | CV),  .) %>% anova()

# EMMs
dt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  mutate(vn2.num = 
        (velumopening_nucleus_on + 
           velumopening_nucleus_off)/2 - Vperiodiconset,
       vn2.den = endv_time - Vperiodiconset,
         vn2 = vn2.num/vn2.den) %>%
  lmer(log(vn2) ~ Dialect + (coda | speaker) + 
         (Dialect | CV),  .) %>% 
  emmeans(., ~ Dialect)
```

Dialect sig., not coda

### Statistics row 2

```{r}
m1 = 
glmmTMB(ratio ~ Dialect * coda + 
       (coda|speaker) +
      (Dialect|CV), 
      family=beta_family(),
      data = propcnas.df %>%
        filter(ratio < 1) %>%
      rename(Dialect = language)) 
summary(m1)
emmeans(m1, ~ Dialect | coda, 
        type = "response") %>% 
  pairs(simple = "each", combine = TRUE)
```

# Intra-gestural tongue tip analysis
The object here is to determine whether the tongue tip gesture is more reduced in /nt/ than in /nd/. 

## Method

```{r}
j = "1.239"
# TT displacement maximum
ttmax = tt.df %>%
  filter(rowind == j) %>%
  pull(alvUS) %>% max()

# Time of TT nucleus onset
ttnon.t = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_nucleus_on) %>% max()

# Time of TT nucleus offset
ttnoff.t = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_nucleus_off) %>% max()


# Time of TT displacement maximum
ttmax.t = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_maxcon_on) %>% max()

# Time of TT raising pk velocity max
velon.tt = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_maxvel_on) %>%
  unique()

# TT raising pk velocity max
ttvel = tt.df %>%
  filter(rowind == j) %>%
   mutate(velocity = 
            500.5 *
            c(stats::filter(alvUS, c(.5, 0, -.5)))) %>%
  pull(velocity) %>% max(na.rm=T)

# time of TT lowering pk velocity max
veloff.tt = tt.df %>%
  filter(rowind == j) %>%
  pull(alveolarconstriction_maxvel_off) %>%
  unique()

# Time of acoustic vowel offset
offset.vowel = tt.df %>%
  filter(rowind == j) %>%
  pull(endv_time) %>%
  unique()

# Time of VELUM pk. velocity raising
offset.velum = tt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxvel_off) %>%
  unique()  

xlim.tt = tt.df %>%
  filter(rowind == j) %>%
  pull(time_in_sequence) %>%
  range()
```


```{r,fig.width = 25/2.54, fig.height=20/2.54}
tt_dis_paper =  tt.df %>%
  filter(rowind == j) %>%
  ggplot +
  aes(y = alvUS, x = time_in_sequence) +
  ylab("Displacement") + 
  geom_line(linewidth = lwd) +
  geom_segment(aes(y = alvUS, x = time_in_sequence), 
               x = (ttnon.t+ttnoff.t)/2, y = 0, lwd=1, 
               xend = (ttnon.t+ttnoff.t)/2, yend = ttmax, 
               arrow=arrow(ends="last", type="closed")) +
  geom_segment(aes(y = alvUS, x = time_in_sequence), 
               x = ttnon.t, y = 0, lwd=1, lty=2, 
               xend = ttnon.t, yend = .9) +
  geom_segment(aes(y = alvUS, x = time_in_sequence), 
               x = ttnoff.t, y = 0, lwd=1, lty=2, 
               xend = ttnoff.t, yend = .9) +
   annotate("text", x = .395, y = .79, label = "1", size=24/.pt)  +
  geom_vline(xintercept = velon.tt, col="slategray", lwd=1) +
geom_segment(aes(y = alvUS, x = time_in_sequence), 
               x = velon.tt, y = .3, 
               xend = (ttnon.t+ttnoff.t)/2, yend = .3, 
             col="black", lwd=1, 
               arrow=arrow(ends="both", type="closed")) +
   annotate("text", x = .335, y = .35, label = "4", size=24/.pt) +
theme_michigan2apaper(k = 24) +
#  geom_vline(xintercept = veloff.tt, col = "slategray", lwd=1) +
scale_colour_manual(values = cols)

tt_vel_paper = tt.df %>%
  filter(rowind == j) %>%
   mutate(velocity = 
            500.5 * 
            c(stats::filter(alvUS, c(.5, 0, -.5)))) %>%
  ggplot +
  aes(y = velocity, x = time_in_sequence) +
  geom_line(linewidth = lwd, col="slategray") +
  coord_cartesian(xlim = xlim.tt) +
  geom_segment(aes(y = velocity, x = time_in_sequence), 
               x = velon.tt, y = 0, 
               xend = velon.tt, yend = .0076 * 500.5, 
               col = "slategray", lwd=1, 
               arrow=arrow(ends="last", 
                                                                                    type="closed",length = unit(0.13, "inches"))) +
geom_hline(yintercept = 0,  col = "slategray") +
geom_segment(aes(y = velocity, x = time_in_sequence), 
               x = veloff.tt, y = 0, 
               xend = veloff.tt, 
             yend = -.006 * 500.5, lwd=1, 
             col="slategray") +
  annotate("text", 
           x = 0.22, 
           y = 0.005 * 500.5, label = "2", size=24/.pt, col="slategray") +
  annotate("text", 
           x = 0.4, 
           y = 0.001 * 500.5, label = "3", size=24/.pt, col="slategray") +
geom_segment(aes(y = velocity, x = time_in_sequence), 
               x = velon.tt, y = 0, 
               xend = veloff.tt, yend = 0, 
              lwd=1, col="slategray", 
               arrow=arrow(ends="both", type="closed")) +
  xlab("Time (s)") + 
   ylab("Velocity (units/s)") +
 theme_michigan2dpaper(k = 24) +
scale_colour_manual(values = cols)

gA.tt <- ggplotGrob(tt_dis_paper)
gB.tt <- ggplotGrob(tt_vel_paper)
grid::grid.newpage()
fig.ttipaper = grid.arrange(rbind(gA.tt, gB.tt))

ggsave(filename = file.path(pfadfigs, 
                            "fig.ttipaper.pdf"), 
      plot = fig.ttipaper, 
      device=cairo_pdf,
      width = 25, 
       height = 20, 
       units = "cm")
```



## Two sets of boxplots for the four parameters
### Param 1 and 2

```{r,fig.width = 30/2.54, fig.height=20/2.54}
lwd = 1.4
ttdisplacement = tt.df %>%
     rename(Dialect = language) %>%
  # extract values in TT plateau 
  filter(time_in_sequence <= alveolarconstriction_nucleus_off) %>%
  filter(time_in_sequence >= alveolarconstriction_nucleus_on) %>%
  group_by(rowind, coda, Dialect,  onset, vowel, speaker) %>%
  # mean value of plateau
summarise(displ = mean(alvUS - alvUS_Vokal_min)) %>%
  ungroup() %>%
     filter(coda %in% c("nt", "nd")) %>%
  ggplot +
  aes(y = displ, x = Dialect, col=coda) +
  geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 3) +
          scale_colour_manual(values = cols.coda) + scale_y_continuous("Displacement", 
                     breaks = c(0, .5, 1)) +
 # coord_cartesian(ylim = c(-3, .5)) +
theme_michigan2(n=26) +
  theme(axis.title.y = element_text(size=26)) +
  theme(axis.text.y = element_text(size=20)) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
  scale_colour_manual(values = cols.coda) 


ttvelocity = 
  tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>% 
  filter(!(is.na(alvUSV_alveolarconstriction_maxvel_onset))) %>%
  rename(Dialect = language)%>%
  ggplot +
  aes(y = alvUSV_alveolarconstriction_maxvel_onset, 
      x = Dialect, col=coda) +
  geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 3) +
  ylab("Velocity (units/s)") +
  xlab("") +
 # coord_cartesian(ylim = c(-3, .5)) +
theme_michigan3(n = 26) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
  scale_colour_manual(values = cols.coda) 

gA <- ggplotGrob(ttdisplacement)
gB <- ggplotGrob(ttvelocity)

grid::grid.newpage()
figttdisplvel= grid.arrange(rbind(gA, gB))

ggsave(filename = file.path(pfadfigs2, "figttdisplvel.pdf"), 
      plot = figttdisplvel, 
      device=cairo_pdf,
       width = 30, 
       height = 20, 
      units = "cm")
```


### Params 3 and 4

```{r,fig.width = 30/2.54, fig.height=20/2.54}
lwd = 1.4
ttartdur = 
  tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>% 
  rename(Dialect = language)%>%
  ggplot +
  aes(y = alveolarconstriction_maxvel_dur, 
      x = Dialect, col=coda) +
  geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 3) +
 # coord_cartesian(ylim = c(-3, .5)) +
theme_michigan2(n = 26) +
   scale_y_continuous("Duration (s)", breaks = c(0.2, .4)) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
  scale_colour_manual(values = cols.coda) 

ttraisingdur = 
  tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>% 
  rename(Dialect = language)%>%
  mutate(param =
           (alveolarconstriction_nucleus_on + alveolarconstriction_nucleus_off)/2 -
           alveolarconstriction_maxvel_on) %>%
  ggplot +
  aes(y = param, 
      x = Dialect, col=coda) +
  geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 3) +
  ylab("Duration (s)") +
  xlab("") +
 # coord_cartesian(ylim = c(-3, .5)) +
theme_michigan3(n = 26) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
  scale_colour_manual(values = cols.coda) 

#gA <- ggplotGrob(ttdisplacement)
#gB <- ggplotGrob(ttvelocity)

#grid::grid.newpage()
#figttdisplvel= grid.arrange(rbind(gA, gB))

figttduration = 
  grid.arrange(ttartdur, ttraisingdur, nrow=2)
ggsave(filename = file.path(pfadfigs2, "figttduration.pdf"), 
      plot = figttduration, 
      device=cairo_pdf,
       width = 30, 
       height = 20, 
      units = "cm")
```



## Comparison of /nd, nt/ on peak tongue tip displacement (1)



Displacment:
`1 .       nd    BRE - USE   0.0788 0.0289  46.9    2.72  0.0363`


```{r}
displ.df = tt.df %>%
     rename(Dialect = language) %>%
  # extract values in TT plateau 
  filter(time_in_sequence <= alveolarconstriction_nucleus_off) %>%
  filter(time_in_sequence >= alveolarconstriction_nucleus_on) %>%
  group_by(rowind, coda, Dialect,  CV, onset, vowel, speaker) %>%
  # mean value of plateau
summarise(displ = mean(alvUS - alvUS_Vokal_min)) %>%
  ungroup() %>%
     filter(coda %in% c("nt", "nd")) 

displ.df %>%
  lmer(displ ~ 
         coda * Dialect + 
         (Dialect | CV) + 
         ( coda| speaker), .) %>%
 step()

# model found
displ.df %>%
  lmer(displ ~
         coda + Dialect + (Dialect | CV) + 
         (1 | speaker) + coda:Dialect, .) %>%
 anova()


displ.df %>%
  lmer(displ ~
         coda + Dialect + (Dialect | CV) + 
         (1 | speaker) + coda:Dialect, .) %>%
emmeans(.,  ~ coda|Dialect) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

```


## Comparison of /nd, nt/ on peak tongue tip velocity raising (param 2)
The results shows that /nt/ is faster than /nd/ in USE only

Statistics
```{r}
tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>% 
  filter(!(is.na(alvUSV_alveolarconstriction_maxvel_onset))) %>%
  rename(Dialect = language)%>%
  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         coda * Dialect + 
         (Dialect | CV) + 
         (coda| speaker), .) %>%
 step()

# simplified model:

tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  # remove 3 NA
  filter(!(is.na(alveolarconstriction_maxvel_on))) %>%
  rename(Dialect = language)%>%
  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         coda + Dialect + (1 | CV) + (1 | speaker) + coda:Dialect, .) %>%
 anova()
```


```{r}
tt.df %>%
  filter(coda %in% c("nt", "nd")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language)%>%
 lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         coda + Dialect + (1 | CV) + (1 | speaker) + coda:Dialect, .) %>%
   emmeans(.,  ~ coda|Dialect) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
```

## Comparison of /nd, nt/ on tongue tip articulatory duration (param 3)

The results show:  BRE > USE

```{r}
tt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  ungroup() %>%
lmer(alveolarconstriction_maxvel_dur ~ 
         coda * Dialect + 
         (Dialect|CV) + (coda|speaker), .) %>%
  step()

tt.df %>%
 filter(coda %in% c("nd", "nt")) %>%
 filter(segment == "Vokal") %>%
 group_by(rowind) %>%
 slice_tail(n=1) %>%
 rename(Dialect = language) %>%
 ungroup() %>%
  lmer(alveolarconstriction_maxvel_dur ~ Dialect + 
         (coda | speaker) + (1 | CV), .) %>%
         anova()

# get EMMs

tt.df %>%
 filter(coda %in% c("nd", "nt")) %>%
 filter(segment == "Vokal") %>%
 group_by(rowind) %>%
 slice_tail(n=1) %>%
 rename(Dialect = language) %>%
 ungroup() %>%
  lmer(alveolarconstriction_maxvel_dur ~ Dialect + 
         (coda | speaker) + (1 | CV), .) %>%
  emmeans(. , ~ Dialect)
```


## Comparison of /nd, nt/ on duration between TT-maxvel-on and TT peak displacement (param 4)

Statistics shows that only dialect is significant.

```{r}
tt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  mutate(d = 
           (alveolarconstriction_nucleus_on + alveolarconstriction_nucleus_off)/2 -
           alveolarconstriction_maxvel_on) %>%
  lmer(d ~ 
         coda * Dialect + 
         (Dialect|CV) + (coda|speaker), .) %>%
  step()

# model found
tt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  mutate(d = 
           (alveolarconstriction_nucleus_on + alveolarconstriction_nucleus_off)/2 -
           alveolarconstriction_maxvel_on) %>%
  lmer(d ~ Dialect + (coda | speaker) + (1 | CV), .) %>%
  anova()

# EMMs

tt.df %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  mutate(d = 
           (alveolarconstriction_nucleus_on + alveolarconstriction_nucleus_off)/2 -
           alveolarconstriction_maxvel_on) %>%
  lmer(d ~ Dialect + (coda | speaker) + (1 | CV), .) %>%
  emmeans(., ~Dialect)
```


# Compensatory lengthening

/nt, nd/ durations
```{r}
finals.ntnd = c("nt", "nd")
durations.df  =  dt.df %>%
  group_by(rowind) %>%
  filter(segment == "Vokal") %>%
   filter(coda %in% finals.ntnd) %>%
  slice_tail(n=1) %>%
  mutate(alldur = velumopening_maxvel_off -  timeseriesoffset, 
         Vdur = endv_time - timeseriesoffset, 
         Ndur = velumopening_maxvel_off - endv_time) %>%
  rename(Dialect = language) 
```

Ndur
```{r}
durations.df %>% 
  group_by(coda) %>%
  summarise(m = mean(Ndur, na.rm=T))
```

```{r}
durations.theme = function(n=24, m = 20)
{
  theme(
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size=m, 
                                   angle = 90, 
                                   vjust = 1, 
                                   hjust=1),
        axis.title.y = element_text(size=n),
        legend.title = element_blank(), 
        legend.position="top",
        text = element_text(size=n))
}
```


```{r}
plot1 = durations.df  %>%
  ggplot +
  aes(y  = alldur,   x = coda, col=Dialect) +
  geom_boxplot() +
  facet_wrap(~ vowel, ncol = 5) +
  ylab("VN duration (s)") +
  xlab("") +
  durations.theme() +
coord_cartesian(ylim = c(0.15, .55)) +
  scale_colour_manual(values = cols.dial) 

plot2 = 
  durations.df %>%
  ggplot +
  aes(y = Vdur, col = Dialect, x = coda) +
  geom_boxplot() +
  durations.theme() +
    theme(
        legend.position = "none",
        axis.text.x = element_text(size=20)) +
  scale_colour_manual(values = cols.dial) +
  facet_wrap(~ vowel, ncol = 5) +
ylab("V duration (s)") +
scale_colour_manual(values = cols.dial)

durplot = grid.arrange(plot1, 
                       plot2, ncol=1, nrow=2)

ggsave(filename = file.path(pfadfigs, "durplot.pdf"), 
    plot = durplot, 
    device=cairo_pdf,
      width = 25, 
       height = 20, 
   units = "cm")

```

Test on all durations
```{r}
durations.df  %>%
  lmer(alldur ~ 
         coda * Dialect + 
         (Dialect|CV) + (coda|speaker), .) %>%
  step()

# model found
# alldur ~ coda + (coda | speaker) + (1 | CV)

durations.df  %>%
  lmer(alldur ~ coda + (coda | speaker) + (1 | CV), .) %>%
  anova()

# Emmeans
durations.df  %>%
  lmer(alldur ~ coda + (coda | speaker) + (1 | CV), .) %>%
  emmeans(., ~coda)
# nd    0.358; nt    0.275
```
Test on vowel duration

```{r}
durations.df  %>%
  lmer(Vdur ~ 
         coda * Dialect + 
         (Dialect|CV) + (coda|speaker), .) %>%
  step()
```

Model found
```{r}
durations.df  %>%
  lmer(Vdur ~ coda + Dialect + (Dialect | CV) + (coda | speaker), .) %>%
  anova()

# Emmeans
durations.df  %>%
  lmer(Vdur ~ coda + Dialect + (Dialect | CV) + (coda | speaker), .) %>%
  emmeans(., ~coda)
# nd    0.207;  nt    0.157

durations.df  %>%
  lmer(Vdur ~ coda + Dialect + (Dialect | CV) + (coda | speaker), .) %>%
  emmeans(., ~Dialect)
# BRE      0.169; USE      0.195
```


# Schematic outlines
## For voiced NC̬

```{r, fig.width=30/2.54, fig.height=30/2.54}
colscheme = c("black", "red", "blue")
ltyscheme = c(1, 3, 3)

library(emuR)
ablack = cr(p=-pi, type="l", 
            N=100, axes=F, lwd=2, ylab="", xlab="", values=T)
blue.x = seq(25, 74, length=100)
blue.y = (ablack/2) 
red.x = seq(12.5, 86.5, length=100)
red.y = (ablack * .75)
black.x = 1:100
black.y = ablack+2


vel.blue = c(stats::filter((ablack/2)-.5, c(.5, 0, -.5)))
temp = vel.blue == max(vel.blue, na.rm=T)
pktime.blue = (1:100)[temp]
pktime.blue = pktime.blue[!is.na(pktime.blue)]
pktime.blue.time= seq(25, 74, length=100)[pktime.blue]
vel.red = c(stats::filter((ablack*.75)-.25, c(.5, 0, -.5)))
temp = vel.red == max(vel.red, na.rm=T)
pktime.red = (1:100)[temp]
pktime.red = pktime.red[!is.na(pktime.red)]
pktime.red.time = seq(12.5, 86.5, length=100)[pktime.red]

pktime.black.time = 57

# horizontal shifts
ktimes = -18
black.x = black.x + ktimes
times = c(blue.x, red.x, black.x)

# vertical shifts

red.y = red.y + 2
blue.y = blue.y + 1.625

y = c(blue.y, red.y, black.y)
y = y + .5
Dialect = c(rep("USE",  length(blue.y)), 
            rep("BRE",  length(red.y)),
            rep("Both", length(black.y)))
scheme.df = data.frame(y, times, Dialect)



schematic_paper = scheme.df %>%
  ggplot +
  aes(y = y, x = times, group=Dialect, col = Dialect,
      linetype = Dialect) +
  theme_light() +
   theme_michigan2hpaper(n = 30, k = 26, position="top") +
  geom_line(linewidth = 1.5) +
  ylab("Displacement") +
  xlab("Time") + 
  coord_cartesian(xlim = c(-15, 87.5), ylim = c(0.5, 3.7)) +
  
  # double arrow for red N
geom_segment(aes(y = y, x = times), 
               x = pktime.red.time, y = 1, 
               xend = pktime.black.time, yend = 1, 
               arrow=arrow(ends="both", type="closed"), 
             col = "red", show.legend = FALSE) +
  
  # double arrow for red V
geom_segment(aes(y = y, x = times), 
               x = -12.5, y = 1, 
               xend = pktime.red.time, yend = 1, 
               arrow=arrow(ends="both", type="closed"), 
             col = "red", show.legend = FALSE) +
  
   # double arrow for blue V
geom_segment(aes(y = y, x = times), 
               x = -12.5, y = .5, 
               xend = pktime.blue.time, yend = .5, 
               arrow=arrow(ends="both", type="closed"), 
             col = "blue", show.legend = FALSE) +
  
     # double arrow for blue N
geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = .5, 
               xend = pktime.black.time, yend = .5, 
               arrow=arrow(ends="both", type="closed"), 
             col = "blue", show.legend = FALSE) +
  
  # vertical blue
  geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = 1.25, 
               xend = pktime.blue.time, yend = 2.125, 
             col = "blue", show.legend = FALSE) +
    geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = 0.5, 
               xend = pktime.blue.time, yend = 0.85, 
             col = "blue", show.legend = FALSE) +
  
  # vertical red
  geom_segment(aes(y = y, x = times), 
               x = pktime.red.time, y = 1, 
               xend = pktime.red.time, yend = 2.5, 
             col = "red", show.legend = FALSE) +
  
  #  vertical black right
    geom_segment(aes(y = y, x = times), 
               x = pktime.black.time, y = .5, 
               xend = pktime.black.time, yend = 2.55, 
             col = "black", show.legend = FALSE) +
  
  # vertical black left
      geom_segment(aes(y = y, x = times), 
               x = -12.5, y = .5, 
               xend = -12.5, yend = 3.5, 
             col = "slategray", show.legend = FALSE) +
  
   annotate("text", x=31.25, y=3.6, 
           label="Velum", 
           size=28/.pt) +
  annotate("text", x=50, y=3.6, 
           label="TT", 
           size=28/.pt, col = "red") +
  annotate("text", x=50, y=2.45, 
           label="TT", 
           size=28/.pt, col = "blue") +
  
   annotate("text", x=44.09343, y=1.1, 
           label="N", 
           size=28/.pt, col = "red") +
  
     annotate("text", x=47.18687, y=.6, 
           label="N", 
           size=28/.pt, col = "blue") +
  
  annotate("text", x=9.343434, y=1.1, 
           label="V", 
           size=28/.pt, col = "red") +
     annotate("text", x=12.43687, y=.6, 
           label="V", 
           size=28/.pt, col = "blue") +

  scale_colour_manual(values = colscheme) +
  scale_linetype_manual(values = ltyscheme)

schematic_paper

ggsave(filename = file.path(pfadfigs2,
                            "fig2schematic.pdf"), 
     plot = schematic_paper, 
     device=cairo_pdf,
       width = 30, 
      height = 20, 
       units = "cm")
```




## For voiceless NC̥


```{r, fig.width=30/2.54, fig.height=30/2.54}
colscheme = c("black", "red", "blue")
ltyscheme = c(1, 3, 3)

library(emuR)
ablack = cr(p=-pi, type="l", 
            N=100, axes=F, lwd=2, ylab="", xlab="", values=T)
#blue.x = seq(25, 74, length=100)
#blue.y = (ablack/2) 
red.x = seq(12.5, 86.5, length=100)
red.y = (ablack * .75)
blue.y = red.y
blue.y.left = blue.y[1:41]
blue.y.right = blue.y[60:100]
blue.y = c(blue.y.left, blue.y.right)


blue.x = red.x[-c(1:9, 92:100)]

black.x = 1:100
black.y = ablack+2



vel.blue = c(stats::filter((ablack/2)-.5, c(.5, 0, -.5)))
temp = vel.blue == max(vel.blue, na.rm=T)
pktime.blue = (1:100)[temp]
pktime.blue = pktime.blue[!is.na(pktime.blue)]
pktime.blue.time= seq(25, 74, length=100)[pktime.blue]
vel.red = c(stats::filter((ablack*.75)-.25, c(.5, 0, -.5)))
temp = vel.red == max(vel.red, na.rm=T)
pktime.red = (1:100)[temp]
pktime.red = pktime.red[!is.na(pktime.red)]
pktime.red.time = seq(12.5, 86.5, length=100)[pktime.red]

pktime.black.time = 57

# horizontal shifts
ktimes = -18
black.x = black.x + ktimes

times = c(blue.x, red.x, black.x)




# vertical shifts

red.y = red.y + 2
blue.y = blue.y + 2
y = c(blue.y, red.y, black.y)
y = y + .5
Dialect = c(rep("USE",  length(blue.y)), 
            rep("BRE",  length(red.y)),
            rep("Both", length(black.y)))
scheme.df = data.frame(y, times, Dialect)

orange.x = seq(31.25, 56.25, length = 51) + 5
orange.y = black.y[1:51]/4 + .25
orange.dial = rep("Both", length(orange.y))
orange.df = data.frame(y = orange.y, 
                       times = orange.x, 
                       Dialect = orange.dial)



schematic_paper = scheme.df %>%
  ggplot +
  aes(y = y, x = times, group=Dialect, col = Dialect,
      linetype = Dialect) +
  theme_light() +
  theme_michigan2hpaper(n = 30, k = 26, position="top") +
  geom_line(linewidth = 1.5) +
  ylab("Displacement") +
  xlab("Time") + 
  coord_cartesian(xlim = c(-15, 87.5), ylim = c(-.5, 4.1)) +
  geom_line(aes(y = y, x= times), 
            data = orange.df, col = "orange", lty=1, linewidth=1.5) + 
  
  # horizontal orange right
  geom_segment(aes(y = y, x = times), 
               x = 61.25, y = 1, 
               xend = 75, yend = 1, 
               linewidth=1.5,
             col = "orange", show.legend = FALSE) +
  
   # horizontal orange left
  geom_segment(aes(y = y, x = times), 
               x = -12.5, y = .5, 
               xend = 36.25, yend = .5, 
               linewidth=1.5,
             col = "orange", show.legend = FALSE) +
  
  
  # double arrow for red voiced N
geom_segment(aes(y = y, x = times), 
               x = pktime.red.time, y = 0, 
               xend = 43.75+5, yend = 0, 
               arrow=arrow(ends="both", type="closed"), 
             col = "red", show.legend = FALSE) +

    # double arrow for red voiceless N
geom_segment(aes(y = y, x = times), 
               x = 43.75+5, y = 0, 
               xend = pktime.black.time, yend = 0, 
               arrow=arrow(ends="both", type="closed"), 
             col = "red", show.legend = FALSE) +
  
  # double arrow for red V
geom_segment(aes(y = y, x = times), 
               x = -12.5, y = 0, 
               xend = pktime.red.time, yend = 0, 
               arrow=arrow(ends="both", type="closed"), 
             col = "red", show.legend = FALSE) +
  
   # double arrow for blue V
geom_segment(aes(y = y, x = times), 
               x = -12.5, y = -.5, 
               xend = pktime.blue.time, yend = -.5, 
               arrow=arrow(ends="both", type="closed"), 
             col = "blue", show.legend = FALSE) +
  
     # double arrow for blue voiced N
geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = -.5, 
               xend = 43.75+5, yend = -.5, 
               arrow=arrow(ends="both", type="closed"), 
             col = "blue", show.legend = FALSE) +

      # double arrow for blue  voiceless N
geom_segment(aes(y = y, x = times), 
               x = 43.75+5, y = -.5, 
               xend = pktime.black.time, yend = -.5, 
               arrow=arrow(ends="both", type="closed"), 
             col = "blue", show.legend = FALSE) +
  
  # vertical blue
  geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = .35, 
               xend = pktime.blue.time, yend = 2.5, 
             col = "blue", show.legend = FALSE) +
  
    geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = -.5, 
               xend = pktime.blue.time, yend = -.1, 
             col = "blue", show.legend = FALSE) +
  
  # vertical red
  geom_segment(aes(y = y, x = times), 
               x = pktime.red.time, y = 2.5, 
               xend = pktime.red.time, yend = 0, 
             col = "red", show.legend = FALSE) +
  
  #  vertical black right
    geom_segment(aes(y = y, x = times), 
               x = pktime.black.time, y = -.5, 
               xend = pktime.black.time, yend = 2.55, 
             col = "black", show.legend = FALSE) +
  
  # vertical black left
      geom_segment(aes(y = y, x = times), 
               x = -12.5, y = -.5, 
               xend = -12.5, yend = 3.5, 
             col = "slategray", show.legend = FALSE) +
  # vertical orange
      geom_segment(aes(y = y, x = times), 
               x = 43.75+5, y = 0.75, 
               xend = 43.75+5, yend = -.5, 
             col = "orange", show.legend = FALSE) +
  
  # vertical double arrow orange
      geom_segment(aes(y = y, x = times), 
               x = 75, y = 0.5, 
               xend = 75, yend = 1, 
                arrow=arrow(ends="both", type="closed"),
             col = "orange", show.legend = FALSE) +
  
  
  
   annotate("text", x=31.25, y=3.8, 
           label="Velum", 
           size=28/.pt) +
  annotate("text", x=50, y=3.8, 
           label="TT", 
           size=28/.pt, col = "red") +
  annotate("text", x=50, y=2.6, 
           label="TT", 
           size=28/.pt, col = "blue") +
    annotate("text", x=83, y=1, 
           label="−voice", 
           size=28/.pt, col = "orange") +
   annotate("text", x=83, y=.5, 
           label="+voice", 
           size=28/.pt, col = "orange") +
  
   annotate("text", x=40, y=0.19, 
           label= "N̬", 
           size=28/.pt, col = "red") +
     annotate("text", x=53, y=0.19, 
           label="N̥",  
           size=28/.pt, col = "red") +
  
     annotate("text", x=43.06187, y=-.3, 
           label="N̬", 
           size=28/.pt, col = "blue") +
  annotate("text", x=53, y=-.3, 
           label="N̥",  
           size=28/.pt, col = "blue") +
  
  annotate("text", x=9.343434, y=0.19, 
           label="V", 
           size=28/.pt, col = "red") +
     annotate("text", x=12.43687, y=-.3, 
           label="V", 
           size=28/.pt, col = "blue") +

  scale_colour_manual(values = colscheme) +
  scale_linetype_manual(values = ltyscheme)

schematic_paper

ggsave(filename = file.path(pfadfigs2,
                            "figscheme_new.pdf"), 
     plot = schematic_paper, 
     device=cairo_pdf,
       width = 30, 
      height = 20, 
       units = "cm")
```

# Figure in responses to reviewers
The proportional time at which the peak velocity of velum lowering occurs relative to the acoustic vowel onset and offset


```{r,fig.width = 30/2.54, fig.height=20/2.54}
onset_norm = dt.df %>%
  rename(Dialect = language) %>%
  filter(coda %in% c("nd", "nt")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(vn2.num = 
           velumopening_maxvel_on - Vowelonset2, 
         vn2.den =
           endv_time - Vowelonset2) %>% 
  ggplot + 
  aes(y = vn2.num/vn2.den, col = coda, x = Dialect) +
  geom_boxplot(lwd=lwd) +
  facet_wrap(~ vowel, ncol = 5) +
 theme_michigan3(n=26) +
theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1)) +
 scale_y_continuous("Proportional time", breaks=c(0.5,  1)) +
 xlab("") +
scale_colour_manual(values = cols.coda)

ggsave(filename = file.path(pfadfigs2, "fig_onset.png"), 
      plot = onset_norm, 
       width = 30, 
       height = 20, 
      units = "cm")
```

The final figure with the addition of the vertical left boundary at the time of the peak velum velocity lowering


```{r, fig.width=30/2.54, fig.height=20/2.54}
colscheme = c("black", "red", "blue")
ltyscheme = c(1, 3, 3)

library(emuR)
ablack = cr(p=-pi, type="l", 
            N=100, axes=F, lwd=2, ylab="", xlab="", values=T)
blue.x = seq(25, 74, length=100)
blue.y = (ablack/2) 
red.x = seq(12.5, 86.5, length=100)
red.y = (ablack * .75)
black.x = 1:100
black.y = ablack+2


vel.blue = c(stats::filter((ablack/2)-.5, c(.5, 0, -.5)))
temp = vel.blue == max(vel.blue, na.rm=T)
pktime.blue = (1:100)[temp]
pktime.blue = pktime.blue[!is.na(pktime.blue)]
pktime.blue.time= seq(25, 74, length=100)[pktime.blue]
vel.red = c(stats::filter((ablack*.75)-.25, c(.5, 0, -.5)))
temp = vel.red == max(vel.red, na.rm=T)
pktime.red = (1:100)[temp]
pktime.red = pktime.red[!is.na(pktime.red)]
pktime.red.time = seq(12.5, 86.5, length=100)[pktime.red]

pktime.black.time = 57

# horizontal shifts
ktimes = -18
black.x = black.x + ktimes
times = c(blue.x, red.x, black.x)

# vertical shifts

red.y = red.y + 2
blue.y = blue.y + 1.625

y = c(blue.y, red.y, black.y)
y = y + .5
Dialect = c(rep("USE",  length(blue.y)), 
            rep("BRE",  length(red.y)),
            rep("Both", length(black.y)))
scheme.df = data.frame(y, times, Dialect)



model_paper = scheme.df %>%
  ggplot +
  aes(y = y, x = times, group=Dialect, col = Dialect,
      linetype = Dialect) +
  theme_light() +
   theme_michigan2xpaper(n = 30, k = 26, position="top") +

  geom_line(linewidth = 1.5) +
  ylab("Displacement") +
  xlab("Time") + 
  coord_cartesian(xlim = c(-15, 87.5), ylim = c(0.5, 3.7)) +
  
  # double arrow for red N
geom_segment(aes(y = y, x = times), 
               x = pktime.red.time, y = 1, 
               xend = pktime.black.time, yend = 1, 
               arrow=arrow(ends="both", type="closed"), 
             col = "red", show.legend = FALSE) +
  
  # double arrow for red V
geom_segment(aes(y = y, x = times), 
               x = -12.5, y = 1, 
               xend = pktime.red.time, yend = 1, 
               arrow=arrow(ends="both", type="closed"), 
             col = "red", show.legend = FALSE) +
  
   # double arrow for blue V
geom_segment(aes(y = y, x = times), 
               x = -12.5, y = .5, 
               xend = pktime.blue.time, yend = .5, 
               arrow=arrow(ends="both", type="closed"), 
             col = "blue", show.legend = FALSE) +
  
     # double arrow for blue N
geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = .5, 
               xend = pktime.black.time, yend = .5, 
               arrow=arrow(ends="both", type="closed"), 
             col = "blue", show.legend = FALSE) +
  
  # vertical blue
  geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = 1.25, 
               xend = pktime.blue.time, yend = 2.125, 
             col = "blue", show.legend = FALSE) +
    geom_segment(aes(y = y, x = times), 
               x = pktime.blue.time, y = 0.5, 
               xend = pktime.blue.time, yend = 0.85, 
             col = "blue", show.legend = FALSE) +
  
  # vertical red
  geom_segment(aes(y = y, x = times), 
               x = pktime.red.time, y = 1, 
               xend = pktime.red.time, yend = 2.5, 
             col = "red", show.legend = FALSE) +
  
  #  vertical black right
    geom_segment(aes(y = y, x = times), 
               x = pktime.black.time, y = .5, 
               xend = pktime.black.time, yend = 2.55, 
             col = "black", show.legend = FALSE) +
  
  #  vertical black left
    geom_segment(aes(y = y, x = times), 
               x = 8, y = .5, 
               xend = 8, yend = 2.5, 
             col = "black", lty=2, show.legend = FALSE) +
  
  # vertical slategray boundary left
      geom_segment(aes(y = y, x = times), 
               x = -12.5, y = .5, 
               xend = -12.5, yend = 3.5, 
             col = "slategray", show.legend = FALSE) +
  
   annotate("text", x=31.25, y=3.6, 
           label="Velum", 
           size=28/.pt) +
  annotate("text", x=50, y=3.6, 
           label="TT", 
           size=28/.pt, col = "red") +
  annotate("text", x=50, y=2.45, 
           label="TT", 
           size=28/.pt, col = "blue") +
  
   annotate("text", x=44.09343, y=1.1, 
           label="N", 
           size=28/.pt, col = "red") +
  
     annotate("text", x=47.18687, y=.6, 
           label="N", 
           size=28/.pt, col = "blue") +
  
  annotate("text", x=9.343434, y=1.1, 
           label="V", 
           size=28/.pt, col = "red") +
     annotate("text", x=12.43687, y=.6, 
           label="V", 
           size=28/.pt, col = "blue") +

  scale_colour_manual(values = colscheme) +
  scale_linetype_manual(values = ltyscheme)

model_paper

ggsave(filename = file.path(pfadfigs2,
                            "model_paper.pdf"), 
     plot = model_paper, 
     device=cairo_pdf,
       width = 30, 
      height = 20, 
       units = "cm")
```
