---
title: "Vowel nasalisation in two varieties of English"
author: ""
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, echo=FALSE}
# knitr::opts_chunk$set(echo = TRUE, fig.width = 10)
```


Load libraries and data:
```{r, message=FALSE}
# for time_in_sequence() function
library(tidyverse)
library(gridExtra)
library(lmerTest)
library(emmeans)
pfad = "/Volumes/vdata/ERC2/Nasals23/analysis"
#pfad = "/Users/jmh/b/current/nasals23/analysis"

# directory for saving figures
# pfadfigs = "/Users/jmh/b/current/nasals23/mar24figs"
pfadfigs = "/Volumes/vdata/ERC2/Nasals23/final/figs"
pfadfigs2 = pfadfigs


load(file = file.path(pfad,       
"erc_mrinasals_concatsigs_prenasalraising_v3_prevowelVokalCodapostcoda_A_tdn__IEAVe__i_pbfszS_01.Rda"))
```
# Figures

Figures in the paper and filenames in this file as follows. Search on the names in brackets e.g. (fig1paper.png) for fig2.png to find the figures. (NB: fig1 was not derived from R commands).

- fig1 	Fig. 1.1 (fig1a.pdf)				
- fig2 	Fig. 3.1. (fig1paper.png)			
- fig3 	Fig. 3.2. (fig5paper.png)		
- fig4  	Fig.4.1 (fig3paper.png)
- fig5 	Fig. 4.2. (fig4paper.png)
- fig6 	Fig.5.1. (fig.ttipaper.png)
- fig7		Fig. 5.2. (figtt_gesture_paper.png)
- fig8		Fig. 6.1. (fig.trade1.paper.png)
- fig9		Fig. 6.2 (tradeextra_paper.png)
- fig10	Fig. 6.3. (tradeall.paper.png)
- fig11	Fig. 7.1. fig.asynch.paper.png
- fig12	Fig. 7.2. (fig.velum_tongue.paper.png)
- fig13	Fig. 7.3. (scheme.png)
- figC.1	Fig. C1 (figUSE_dial.png)
- figD.1	Fig. D1 (figappb.png)
- figE.1	Fig. E1 (figVN.png)
- figF.1	Fig. F1 (velum_tongue2.png)
- figF.2	Fig. F2 (veltongueint.png)


# Preliminaries
## 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_rs_addlang.Rda"))
dim(df_erc_prenasalraising)

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

#!! warning : very bad solution with timeseriesoffset!!

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_rs_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,
# position at gesture onset
alvUS_alveolarconstriction_gesture_onset,
# 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,
         # this works because acoustic onset
         # of the vowel is zero for all segments
         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")
```


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) 
```
## 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: 1100
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)
```

## Define 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_michigan1k <-
function(n = 24){
  theme(
axis.text.x = element_text(size=n), 
axis.title.x = element_text(size=n), 
axis.text.y = element_blank(), 
axis.title.y = element_blank(), 
text = element_text(size=n))
}
theme_michigan2 <-
function(n = 24){
    theme(
        axis.text.x = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_text(size=24),
        plot.title = element_text(size = 32),
        legend.position="top",
        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_michigan2bpaper <-
function(n = 24, k = 18){
    theme(
        axis.text.x = element_text(size=k),
        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="bottom",
legend.title = element_blank(),
        text = element_text(size=n))
}
theme_michigan2cpaper <-
function(n = 24, k = 18, 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),
        plot.title = element_text(size = n),
        legend.position=position,
        text = element_text(size=n))
}
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="top",
legend.title = element_blank(),
        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_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_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_michigan2ipaper =
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))
    }
theme_michigan3 <-
function(n = 24){
    theme(
        legend.position = "none",
        axis.text.y = element_text(size=24),
        plot.title = element_text(size = 32),
        text = element_text(size=n))
}
theme_michigan4 <-
function(n = 24){
  theme(
legend.position = "none",
plot.title = element_text(size = 32),
axis.text = element_blank(),
axis.title = element_blank())
}
theme_michigan5 <-
function(n = 32){
  theme(
legend.position = "none",
plot.title = element_text(size = 32),
axis.text.y = element_blank(),
axis.title.y = element_blank(), 
text = element_text(size=n))
}
theme_michigan6 <-
function(n = 32){
    theme(
        legend.position = "none",
        strip.text.x = element_text(size = 32),
        plot.title = element_text(size = 36),
        axis.title.x = element_text(size=n),
        axis.text.x = element_text(size=24),
        axis.text.y = element_blank(), 
        text = element_text(size=n))
}
theme_michigan6x =
function(n = 32){
    theme(
        legend.position = "top",
legend.title = element_blank(),
        strip.text.x = element_text(size = 32),
        plot.title = element_text(size = 36),
        axis.title.x = element_text(size=n),
        axis.text.x = element_text(size=24),
        axis.text.y = element_text(size=24), 
        text = element_text(size=n))
}
theme_michigan6dial =
function(n = 32){
        theme(
            legend.position = "top",
            legend.title = element_blank(),
            legend.text = element_text(size = n+4),
            strip.text.x = element_text(size = n),
            plot.title = element_text(size = 36),
            axis.title.x = element_text(size=n),
            axis.text.x = element_text(size=24),
            axis.text.y = element_blank(), 
            text = element_text(size=n))
    }

theme_michigan7 <-
function(n = 32){
  theme(
legend.position = "none",
plot.title = element_text(size = 32),
axis.text.y = element_blank(), 
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_michigan11 = 
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.title = element_text(size = n),
        legend.position=position,
        text = element_text(size=n))
}
```
# Nasalization relative to the acoustic VN boundary (section 3)

The following shows the boundary of the acoustic vowel offset, the time of the peak displacement of the velum, and the peak velocity times also of the velum

```{r, fig.width = 30/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()
# Times of pk vel. displacement
velon = dt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxvel_on) %>%
  unique()
veloff = dt.df %>%
  filter(rowind == j) %>%
  pull(velumopening_maxvel_off) %>%
  unique()

# get acoustic vowel offset
endtime = dt.df %>%
  filter(rowind == j) %>%
  pull(endv_time) %>%
  unique()

# get data to draw the area
data1 = dt.df %>%
  filter(rowind == j) %>%
  filter(time_in_sequence >= endv_time) %>%
  filter(time_in_sequence <= velumopening_maxvel_off)

data2 = dt.df %>%
  filter(rowind == j) %>%
  filter(time_in_sequence <= endv_time) %>%
  filter(time_in_sequence >= velumopening_maxvel_on)
```

## fig2 [Fig. 3.1. (fig1paper.png)].
```{r,fig.width = 30/2.54, fig.height=20/2.54}
lwd = 1.5
j = "24.63"
fig1paper = dt.df %>%
  filter(rowind == j) %>%
  ggplot +
  aes(y = velum2US, x = time_in_sequence) +
geom_segment(aes(y = velum2US, 
                 x = time_in_sequence), 
             x = velon, y = 0, xend = velon,
             yend = .8, lty=2, lwd=lwd, col="black") +
geom_segment(aes(y = velum2US, 
                 x = time_in_sequence), 
             x = endtime, y = 0, xend = endtime,
             yend = .8,  lwd=lwd, col="black") +
geom_segment(aes(y = velum2US, 
                 x = time_in_sequence), 
             x = veloff, y = 0, xend = veloff,
             yend = .8, lty=2, lwd=lwd, col="black") +
geom_segment(aes(y = velum2US, 
                 x = time_in_sequence), 
             x = maxcon, y = 0, xend = maxcon,
             yend = .8, lty=2, lwd=lwd, col="black") +
annotate("text", x=0.3, y=0.25, 
           label="italic(i[1])", 
           parse=TRUE, 
           size=32/.pt) +
  annotate("text", x=0.28, y=0.13, 
           label="italic(i[2])", 
           parse=TRUE, 
           size=32/.pt) +
  # ggtitle("Fig. 0.0") +
  xlab("Time (seconds)") +
  ylab("Size of velum lowering") + 
  geom_area(
    aes(y = velum2US, x = time_in_sequence), 
        data = data1, 
    fill = 4,
            alpha = 0.5,
            color = 1,    
            lwd = 0.5,    
            linetype = 1) +
  geom_area(
    aes(y = velum2US, x = time_in_sequence), 
        data = data2, 
    fill = "orange",
            alpha = 0.5,
            color = 1,    
            lwd = 0.5,    
            linetype = 1) +
geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = endtime, y = .2, 
               xend = maxcon, yend = .2, 
               arrow=arrow(ends="both", type="closed")) +
geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = velon, y = .1, 
               xend = veloff, yend = .1, lty=2, 
               arrow=arrow(ends="both", type="closed")) +
  annotate("text", x=c(velon,  endtime, maxcon, veloff), 
           y = rep(.86, 4),
           label=c("italic(t[1])", "italic(t[2])", 
                   "italic(t[3])", "italic(t[4])"),
           parse=TRUE, 
           size=32/.pt) +

theme_michigan1(n = 32) +
scale_colour_manual(values = cols)
fig1paper
ggsave(filename = file.path(pfadfigs2, "fig2.png"), 
       plot = fig1paper, 
      width = 30, 
       height = 20, 
       units = "cm")
```


## fig3 (Fig. 3.2. fig5paper.png)

An analysis of nasal peak delay (`peak_prop`) and of area.

Sum of nasalisation in coda

```{r}
propc.df = dt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
    filter(time_in_sequence >= endv_time) %>%
   filter(time_in_sequence <= velumopening_maxvel_off) %>%
  group_by(rowind, vowel, coda, language, onset, speaker) %>%
  summarise(cprop = sum(velum2US)) %>%
  ungroup()
```

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

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

Fig.3 for the paper
```{r,fig.width = 30/2.54, fig.height=20/2.54}
cols = c("red", "blue")

fig2apaper = dt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval = 
           velumopening_maxvel_off - velumopening_maxvel_on,
         peak_time = velumopening_maxcon_on - endv_time,
         peak_prop = peak_time/intval) %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  ggplot + 
  aes(y = peak_prop, col = Dialect, x = coda) +
  geom_boxplot() +
 scale_y_continuous(breaks=c(-1, 0, 1)) +
  facet_wrap(~ vowel, ncol = 5) +
  ylab("") +
  xlab("") +
  ggtitle("Proportional alignment of peak velum lowering") +
  coord_cartesian(ylim = c(-1, 1)) + 
  geom_hline(yintercept=0, lty=2, linewidth=1.1) +
 theme_michigan2apaper(n = 20) +
  ylab("Offset (proportion)") +
scale_colour_manual(values = cols)


fig2bpaper = propcnas.df %>%
  rename(Dialect = language) %>% 
  ggplot +
  aes(y = log(ratio), x = coda, col=Dialect) +
  geom_boxplot() +
  geom_hline(yintercept=log(.5), lty=2, linewidth=1.1) + 
  facet_wrap(~ vowel, ncol = 5) +
  ylab("Area (proportion)") +
  xlab("") +
  coord_cartesian(ylim = c(-3, .5)) +
ggtitle("Log. proportion  of nasalization in N") +
  theme_michigan2bpaper(n = 20) +
  scale_colour_manual(values = cols) 

fig5paper = 
  grid.arrange(fig2apaper, fig2bpaper, nrow=2)
fig5paper
ggsave(filename = file.path(pfadfigs2, "fig3.png"), 
  plot = fig5paper, 
      width = 30, 
      height = 20, 
     units = "cm")
```


### Statistics associated with fig2

```{r}
p.df = dt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval = 
           velumopening_maxvel_off - velumopening_maxvel_on,
         peak_time = velumopening_maxcon_on - endv_time,
         peak_prop = peak_time/intval) %>%
  rename(Dialect = language) %>%
  ungroup() 
#p.df %>%
#lmer(peak_prop ~ Dialect * vowel * coda + 
#       (vowel+coda|speaker) +
#       (1|onset), .) %>% step()
# simplified model also won't converge

#p.df %>%
#lmer(peak_prop ~ Dialect + 
#       vowel + coda + 
#       (vowel + coda | speaker) + 
#       vowel:coda, .) 

p.df %>%
lmer(peak_prop ~ Dialect + 
       vowel + coda + 
       (1 | speaker) + 
       vowel:coda, .) %>% anova()

p.em.df = p.df %>%
lmer(peak_prop ~ Dialect + 
       vowel + coda + 
       (1 | speaker) + 
       vowel:coda, .) %>% 
  emmeans(.,  ~ vowel | coda) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
print(p.em.df, n = nrow(p.em.df))
```
EMM for dialect
```{r}
p.df %>%
lmer(peak_prop ~ Dialect + 
       vowel + coda + 
       (1 | speaker) + 
       vowel:coda, .) %>% 
  emmeans(.,  ~ Dialect)
```
EMM for vowel
```{r}
p.df %>%
lmer(peak_prop ~ Dialect + 
       vowel + coda + 
       (1 | speaker) + 
       vowel:coda, .) %>% 
  emmeans(.,  ~ vowel)
```
EMM for coda
```{r}
p.df %>%
lmer(peak_prop ~ Dialect + 
       vowel + coda + 
       (1 | speaker) + 
       vowel:coda, .) %>% 
  emmeans(.,  ~ coda)
```

Statistics Fig. 2 bottom row
```{r}
# Identify and remove NA (there are 2 of them)
rowind.rem = propcnas.df %>%
  filter(is.na(ratio)) %>%
  pull(rowind) %>% unique()
rowind.rem

propcnas.df = propcnas.df %>%
  filter(!(rowind %in% rowind.rem))
```
does not converge
```{r}
#propcnas.df %>% 
#rename(Dialect = language) %>%
#lmer(log(ratio) ~ vowel + coda + Dialect + 
#       vowel:coda + vowel:Dialect + coda:Dialect+
#  (vowel + coda|speaker) +  
#(Dialect|onset), .) %>%
#step()
```

Convergence in next code snippet

Main effects:

- Vowel F[4, 928] = 591.59, p < 0.001
- Coda: F[2, 1857.0] = 21.91, p < 0.001
- Dialect F[1, 41.02] = 28.73, p < 0.001
- vowel x coda interaction F[8, 1853.76] = 8.38, p < 0.001

- area of coda nasal  influenced by vowel: æ < eɪ < /ʌ, ɛ, ɪ/ 
- additionally area of coda nasal for /n, nd/: /ʌ, ɛ < ɪ
- nd < nz < n in context of æ
- nd < /nz, n/ in context of eɪ

```{r}
propcnas.df %>% 
rename(Dialect = language) %>%
lmer(log(ratio) ~ vowel + coda +
     Dialect + 
       (1  | speaker) + 
       (1 | onset) + 
       vowel:coda, .) %>% anova()

propc.em.df = propcnas.df %>% 
rename(Dialect = language) %>%
lmer(log(ratio) ~ vowel + coda + Dialect + (1  | speaker) + (1 | onset) + vowel:coda, .) %>% 
emmeans(.,  ~ vowel | coda) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4)) 
print(propc.em.df, n = nrow(propc.em.df))
```

Estimated marginal means for dialect
```{r}
propcnas.df %>% 
rename(Dialect = language) %>%
lmer(log(ratio) ~ vowel + coda + Dialect + (1  | speaker) + (1 | onset) + vowel:coda, .) %>% 
emmeans(.,  ~ Dialect)
```
Estimated marginal means for vowel

```{r}
propcnas.df %>% 
rename(Dialect = language) %>%
lmer(log(ratio) ~ vowel + coda + Dialect + (1  | speaker) + (1 | onset) + vowel:coda, .) %>% 
emmeans(.,  ~ vowel)
```
Estimated marginal means for coda

```{r}
propcnas.df %>% 
rename(Dialect = language) %>%
lmer(log(ratio) ~ vowel + coda + Dialect + (1  | speaker) + (1 | onset) + vowel:coda, .) %>% 
emmeans(.,  ~ coda)
```

# Intra-gestural velum analysis (section 4)

```{r fig.width=20/2.54, fig.height=30/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 = c(stats::filter(velum2US, c(.5, 0, -.5)))) %>%
  pull(velocity) %>% max(na.rm=T)
```

## fig4 (Fig.4.1 (fig3paper.png)).

```{r, fig.width=25/2.54, fig.height=20/2.54}
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 = maxcon, y = plateau, 
               xend = maxcon, 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 = 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_hline(yintercept = 0,  col = "slategray") +
  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") +
  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()
fig3paper = grid.arrange(rbind(gA, gB))
fig3paper
ggsave(filename = file.path(pfadfigs, "fig4.png"), 
    plot = fig3paper, 
      width = 25, 
       height = 20, 
   units = "cm")
```

## fig5 (Fig. 4.2. (fig4paper.png))

```{r, fig.width=30/2.54, fig.height=30/2.54}
z.df = dt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  mutate(reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  endv_time), 5), 
         coda = factor(coda))

# remove 5 segments that have NA at displacement gesture onset
z.df = z.df %>% 
  filter(!is.na(velum2US_velumopening_gesture_onset)) 
z.df = z.df %>%
  mutate(velum2US = 
           velum2US - velum2US_velumopening_gesture_onset)
         
fig4paper =  z.df %>%
group_by(language, coda, vowel, reltime)  %>%
    summarise(vel = mean(velum2US)) %>%
  ungroup() %>%
  rename(Dialect = language) %>%
    ggplot + 
    aes(y = vel, x = reltime, col = Dialect, 
        group = interaction(Dialect, coda, vowel)) +
  geom_line(linewidth=1.5) +
  scale_x_continuous("Time (ms)", breaks = c(-200, 0, 200)) +
    ylab("Displacement") +
  xlab("Time (ms)") +
    facet_wrap(coda ~ vowel, ncol=5) +
    coord_cartesian(xlim = c(-250, 250)) +
    geom_vline(xintercept=0, lty=2) +
   theme_michigan2dpaper(n = 32, k=28) +
  scale_colour_manual(values = cols) 
 fig4paper
ggsave(filename = file.path(pfadfigs2, "fig5.png"), 
      plot = fig4paper, 
       width = 45, 
       height = 30, 
      units = "cm")


```



### Statistics indirectly related to Fig. 4

#### Velum opening displacement maximum: Velumdmax
Displacement is averaged between the peak velocities of 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("n", "nd", "nz")) %>%
  filter(time_in_sequence >= velumopening_nucleus_on) %>%
   filter(time_in_sequence <= velumopening_nucleus_off) %>%
  group_by(rowind, language, onset,
           coda, vowel, speaker) %>%
   summarise(displ = mean(velum2US)) %>%
  rename(Dialect = language) %>%
  ungroup() 

#displ.df %>%
#     lmer(displ ~ 
#              coda * vowel * Dialect   
#              (1|onset)  
#              (vowel + coda| speaker), .) %>%
#     step()

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

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

EMM for vowel
```{r}
displ.df %>%
lmer(displ ~ coda + vowel + 
       Dialect + (vowel + coda | speaker) + vowel:Dialect, .) %>%
  emmeans(.,  ~ vowel)
```

EMM for coda

```{r}
displ.df %>%
lmer(displ ~ coda + vowel + 
       Dialect + (vowel + coda | speaker) + vowel:Dialect, .) %>%
  emmeans(.,  ~ coda)
```

Test for coda differences: n > nd; no nz > nd, no {n, nz} differences
```{r}
displ.em.df = displ.df %>%
lmer(displ ~ coda + vowel + 
       Dialect + (vowel + coda | speaker) + vowel:Dialect, .) %>%
  emmeans(.,  ~ coda) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
print(displ.em.df, n= nrow(displ.em.df))
```

#### peak velocity of opening: Velumvmaxo
```{r}

pklower.df = dt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  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 * vowel * Dialect +
#       (1|onset) +
#     (vowel + coda | speaker), .) %>%
#step()
# model found, won't converge
#pklower.df %>%
#lmer(velum2USV_velumopening_maxvel_onset ~ 
#       coda + vowel + Dialect + (1 | onset) + (vowel | speaker) + #vowel:Dialect, .) %>%
#anova()

pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ 
       coda + vowel + Dialect +  (1|onset) + (1| speaker) + vowel:Dialect, .) %>%
anova()

pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ 
       coda + vowel + Dialect +  (1|onset) + (1| speaker) + vowel:Dialect, .) %>%
emmeans(.,  ~ Dialect|vowel) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
```
EMM for vowel
```{r}
pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ 
       coda + vowel + Dialect +  (1|onset) + (1| speaker) + vowel:Dialect, .) %>%
emmeans(.,  ~ vowel) 
```

EMM for coda
```{r}
pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ 
       coda + vowel + Dialect +  (1|onset) + (1| speaker) + vowel:Dialect, .) %>%
emmeans(.,  ~ coda) 
```
Test for coda differences: nd, nz > n
```{r}
pklower.df %>%
lmer(velum2USV_velumopening_maxvel_onset ~ 
       coda + vowel + Dialect +  (1|onset) + (1| speaker) + vowel:Dialect, .) %>%
emmeans(.,  ~ coda) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
```



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

#artdur.df %>%
#  lmer(art_time ~ vowel * coda * Dialect +
#         (1|onset) +
#         (vowel + coda|speaker), .) %>%
# step()
# model found: gives singular
#artdur.df %>%
#  lmer(art_time ~ vowel + coda + 
#         (1 | onset) + 
#         (vowel + coda | speaker) + 
#         vowel:coda, .) %>% anova()

artdur.df %>%
  lmer(art_time ~ vowel + coda + 
         (1 | onset) + 
         (1  | speaker) + 
         vowel:coda, .) %>% anova()

# Test significant `vowel * coda` interaction:

artdur.em.df = artdur.df %>%
  lmer(art_time ~ vowel + coda + 
         (1 | onset) + 
         (1  | speaker) + 
         vowel:coda, .) %>% 
  emmeans(.,  ~ coda | vowel) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
print(artdur.em.df, n = nrow(artdur.em.df))

```
EMM means: vowel
```{r}
 artdur.df %>%
  lmer(art_time ~ vowel + coda + 
         (1 | onset) + 
         (1  | speaker) + 
         vowel:coda, .) %>% 
  emmeans(.,  ~ vowel)
```
EMM means: coda
```{r}
 artdur.df %>%
  lmer(art_time ~ vowel + coda + 
         (1 | onset) + 
         (1  | speaker) + 
         vowel:coda, .) %>% 
  emmeans(.,  ~ coda)
```

# Intra-gestural tongue tip analysis (section 5)

```{r, fig.width=20/2.54, fig.height=30/2.54}
j = 2.5
# TT displacement maximum
ttmax = tt.df %>%
  filter(rowind == j) %>%
  pull(alvUS) %>% max()

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

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

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

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

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

## fig6  (Fig.5.1. (fig.ttipaper.png))
```{r, fig.width=25/2.54, fig.height=20/2.54}

tongue_tip_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 = velon.tt, y = 0.4, 
               xend = ttmax.t, yend = 0.4, 
             col="black", lwd=lwd,
               arrow=arrow(ends="both", type="closed")) +
geom_segment(aes(y = alvUS, x = time_in_sequence), 
               x = ttmax.t, y = 0, lwd=lwd,
               xend = ttmax.t, yend = ttmax, 
               arrow=arrow(ends="last", type="closed")) +
  geom_vline(xintercept = c(velon.tt), 
             lty=2, linewidth=lwd,
             col = c("black")) +
  annotate("text", x = .4, y = .6, label = "1", size=24/.pt) +
  annotate("text", 
           x = mean(c(velon.tt, ttmax.t)), 
           y = 0.45, label = "2", size=24/.pt) +
  coord_cartesian(xlim = xlim.tt) +
 theme_michigan2apaper(k = 24) +
scale_colour_manual(values = cols)

tt_vel_paper = tt.df %>%
  filter(rowind == j) %>%
   mutate(velocity = c(stats::filter(alvUS, c(.5, 0, -.5)))) %>%
  ggplot +
  aes(y = velocity, x = time_in_sequence) +
  geom_line(linewidth = lwd) +
  coord_cartesian(xlim = xlim.tt) +
  geom_vline(xintercept = c(velon.tt, veloff.tt), 
             lty=2, linewidth=lwd,
             col = c("black", "black")) +
geom_hline(yintercept = 0,  col = "slategray") +
geom_segment(aes(y = velocity, x = time_in_sequence), 
               x = velon.tt, y = ttvel, 
               xend = velon.tt, yend = 0, lwd=lwd,
               arrow=arrow(ends="first", type="closed")) +
  annotate("text", 
           x = 0.31, 
           y = 0.0035, label = "4", size=24/.pt) +
  annotate("text", 
           x = 0.4, 
           y = 0.001, label = "3", size=24/.pt) +
geom_segment(aes(y = velocity, x = time_in_sequence), 
               x = velon.tt, y = 0, 
               xend = veloff.tt, yend = 0, 
             col="black", lwd=lwd,
               arrow=arrow(ends="both", type="closed")) +
  ylab("Velocity") +
  xlab("Time (s)") + 
 theme_michigan2dpaper(k = 24) +
scale_colour_manual(values = cols)

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

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

## fig7 (Fig. 5.2. (figtt_gesture_paper.png))

```{r, fig.width=30/2.54, fig.height=30/2.54}
figtt_gesture_paper = tt.df %>%
filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>% 
  rename(Dialect = language) %>%
  ggplot + 
  aes(x = alveolarconstriction_maxvel_dur,  
      y = alvUS_alveolarconstriction_maxcon_onset -
        alvUS_alveolarconstriction_gesture_onset, 
      col=Dialect) +
  theme_michigan2epaper(position="top") + 
  xlab("Duration (s) between TT peak velocities") +
  ylab("TT displacement magnitude") +
  geom_point(size = 3) +
  facet_wrap(coda ~ vowel) +
scale_colour_manual(values = cols)
figtt_gesture_paper

ggsave(filename = file.path(pfadfigs2, "fig7.png"), 
plot = figtt_gesture_paper, 
       width = 30, 
       height = 30, 
      units = "cm")
```


### Statistics
#### peak tongue tip displacement
```{r}
tonguetip.df = tt.df %>%
filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>% 
  rename(Dialect = language) %>%
  mutate(displ = alvUS_alveolarconstriction_maxcon_onset -
           alvUS_alveolarconstriction_gesture_onset)

#tonguetip.df %>%
#  lmer(displ ~ vowel * coda * Dialect + 
#         (1|onset) + (vowel + coda|speaker), .) %>%
#  step()

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

EMM means: dialect
```{r}
tonguetip.df %>%
  lmer(displ ~ vowel + 
         coda + 
         Dialect + 
         (1 | onset) + 
         (vowel + coda | speaker) + 
         vowel:coda + vowel:Dialect, .) %>%
  emmeans(.,  ~ Dialect)
```
EMM means: vowel
```{r}
tonguetip.df %>%
  lmer(displ ~ vowel + 
         coda + 
         Dialect + 
         (1 | onset) + 
         (vowel + coda | speaker) + 
         vowel:coda + vowel:Dialect, .) %>%
  emmeans(.,  ~ vowel)
```


EMM means: coda
```{r}
tonguetip.df %>%
  lmer(displ ~ vowel + 
         coda + 
         Dialect + 
         (1 | onset) + 
         (vowel + coda | speaker) + 
         vowel:coda + vowel:Dialect, .) %>%
  emmeans(.,  ~ coda)
```

Post-hoc tests
```{r}
tonguetip.df %>%
  lmer(displ ~ vowel + 
         coda + 
         Dialect + 
         (1 | onset) + 
         (vowel + coda | speaker) + 
         vowel:coda + vowel:Dialect, .) %>%
  emmeans(.,  ~ Dialect|vowel) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

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

```

#### peak tongue tip velocity

```{r}
#tonguetip.df %>%
#  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
#         vowel * coda * Dialect + 
#         (1|onset) + (vowel + coda|speaker), .) %>%
#  step()
# model found - won't converge
#tonguetip.df %>%
#  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
#         vowel + 
#         coda + 
#         (1 | onset) + 
 #        (vowel + coda | speaker), .) %>%
 # anova()

tonguetip.df %>%
  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         vowel + 
         coda + 
         (1 | onset) + 
         (1 | speaker), .) %>%
  anova()
```
EMM means: vowel
```{r}
tonguetip.df %>%
  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         vowel + 
         coda + 
         (1 | onset) + 
         (1 | speaker), .) %>%
     emmeans(.,  ~ vowel)

```
EMM means: coda

```{r}
tonguetip.df %>%
  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         vowel + 
         coda + 
         (1 | onset) + 
         (1 | speaker), .) %>%
     emmeans(.,  ~ coda)

```
Post-hoc tests
```{r}  
# test for vowel
tonguetip.df %>%
  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         vowel + 
         coda + 
         (1 | onset) + 
         (1 | speaker), .) %>%
     emmeans(.,  ~ vowel) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

# test for coda
tonguetip.df %>%
  lmer(alvUSV_alveolarconstriction_maxvel_onset ~ 
         vowel + 
         coda + 
         (1 | onset) + 
         (1 | speaker), .) %>%
  emmeans(.,  ~ coda) %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))
```

#### duration between tongue tip peak velocities
alveolarconstriction_maxvel_dur
```{r}
#tonguetip.df %>%
#  lmer(alveolarconstriction_maxvel_dur ~ 
#         vowel * coda * Dialect + 
#         (1|onset) + (vowel + coda|speaker), .) %>%
#  step()
# model found
tonguetip.df %>%
  lmer(alveolarconstriction_maxvel_dur ~ 
         coda + 
         Dialect + 
         (1 | onset) + 
         (coda | speaker), .) %>%
         anova()
```
EMM mean dialect:
```{r}
tonguetip.df %>%
  lmer(alveolarconstriction_maxvel_dur ~ 
         coda + 
         Dialect + 
         (1 | onset) + 
         (coda | speaker), .) %>%
  emmeans(.,  ~ Dialect)

```
EMM means coda

```{r}
tonguetip.df %>%
  lmer(alveolarconstriction_maxvel_dur ~ 
         coda + 
         Dialect + 
         (1 | onset) + 
         (coda | speaker), .) %>%
  emmeans(.,  ~ coda) 
```


#### Interval between time of maximum peak velocity and peak displacement

```{r}
#tt.df %>%
#  filter(coda %in% c("n", "nd", "nz")) %>%
#  filter(segment == "Vokal") %>%
#  group_by(rowind) %>%
#  slice_tail(n=1) %>%
#  mutate(ttint = 
#           alveolarconstriction_maxcon_on - 
#           alveolarconstriction_maxvel_on)    %>%
#  rename(Dialect = language) %>%
#  ungroup() %>% 
#  lmer(ttint ~ vowel * coda * Dialect +
#         (Dialect|stem) + (vowel+coda|speaker), .) %>%
#  step()
# 
tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(ttint = 
           alveolarconstriction_maxcon_on - 
           alveolarconstriction_maxvel_on)    %>%
  rename(Dialect = language) %>%
  ungroup() %>% 
  lmer(ttint ~ vowel + coda + Dialect + (coda | speaker) + (1 | stem), .) %>%
  anova()
```
Vowel, coda, dialect all sig in the above.

Analysis of vowel and coda
coda: n < /nd, nz/. vowel: ʌ >  ɪ
```{r}
# estimated marginal means for coda
tt.df.i3.coda.emmeans =
tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(ttint = 
           alveolarconstriction_maxcon_on - 
           alveolarconstriction_maxvel_on)    %>%
  rename(Dialect = language) %>%
  ungroup() %>% 
  lmer(ttint ~ vowel + coda + 
         Dialect + (coda | speaker) + (1 | stem), .) %>%
  emmeans(.,  ~ coda)
tt.df.i3.coda.emmeans
# signficance of coda
tt.df.i3.coda.emmeans %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

# estimated marginal means for vowel
tt.df.i3.vowel.emmeans =
tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(ttint = 
           alveolarconstriction_maxcon_on - 
           alveolarconstriction_maxvel_on)    %>%
  rename(Dialect = language) %>%
  ungroup() %>% 
  lmer(ttint ~ vowel + coda + 
         Dialect + (coda | speaker) + (1 | stem), .) %>%
  emmeans(.,  ~ vowel) 
tt.df.i3.vowel.emmeans 

# significance of vowel
tt.df.i3.vowel.emmeans %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

# estimated marginal means for Dialect
tt.df.i3.dialect.emmeans =
tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(ttint = 
           alveolarconstriction_maxcon_on - 
           alveolarconstriction_maxvel_on)    %>%
  rename(Dialect = language) %>%
  ungroup() %>% 
  lmer(ttint ~ vowel + coda + 
         Dialect + (coda | speaker) + (1 | stem), .) %>%
  emmeans(.,  ~ Dialect) 
tt.df.i3.dialect.emmeans 
```


# Vowel nasalization and tongue tip lenition (section 6)

```{r}
# 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()

# set the time scale to be the same for
# TT and velum
 xlim.trade  = tt.df %>%
  filter(rowind == j) %>%
  pull(time_in_sequence) %>%
  range()
```

## fig8  (fig.trade1.paper.png)
```{r, fig.width=30/2.54, fig.height=30/2.54}
j = 2.5
velum_trade_paper = tt.df %>%
  filter(rowind == j) %>%
  ggplot +
  aes(y = velum2US, x = time_in_sequence) +
  geom_line(linewidth = lwd) +
    coord_cartesian(xlim = xlim.trade) +
  annotate("text", 
           x=mean(c(velon.tt, maxcon.vel)), y=0.57, 
           label="italic(i[1])", 
           parse=TRUE, 
           size=24/.pt) +
  annotate("text", x=.35, y=0.25, 
           label="italic(i[2])", 
           parse=TRUE, 
           size=23/.pt) +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = velon.vel, y = .15, 
               xend = velon.vel, yend = .85, lty=2, lwd=lwd) +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = velon.tt, y = .15, 
               xend = velon.tt, yend = .85, col = "slategray", 
               lty=2, lwd=lwd) +
geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = maxcon.vel, y = .15, 
               xend = maxcon.vel, yend = .85, lty=2, lwd=lwd) +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = veloff.vel, y = .15, 
               xend = veloff.vel, yend = .85, lty=2, lwd=lwd) +
annotate("text", 
         x=c(velon.vel, velon.tt, maxcon.vel, veloff.vel), 
           y = rep(.9, 4),
           label=c("italic(t[1])", "italic(t[2])", 
"italic(t[3])", "italic(t[4])"),
           parse=TRUE, 
           size=24/.pt) +
  ylab("Displacement") +

ggtitle("Velum") +
  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 = velon.vel, y = .2, 
               xend = veloff.vel, yend = .2, col="black",
               arrow=arrow(ends="both", type="closed")) +
theme_michigan2fpaper(n= 28, k = 24)

# draw tongue tip with kinematic markers
tongue_trade_tip_paper = tt.df %>%
  filter(rowind == j) %>%
  ggplot +
  aes(y = alvUS, x = time_in_sequence) +
  geom_line(col="slategray",  linewidth = lwd) +
  coord_cartesian(xlim = xlim.trade) +
  scale_x_continuous("Time (s)", breaks = c(0.1, 0.3, 0.5)) +
  geom_vline(xintercept = c(velon.tt), 
             lty=2, linewidth = lwd, 
             col = c("slategray")) +
geom_segment(aes(y = alvUS, x = time_in_sequence), 
               x = ttmax.t, y = 0, col = c("slategray"),
               xend = ttmax.t, yend = ttmax,
               arrow=arrow(ends="last", type="closed")) +
 # annotate("text", x=0.01, y=0.79, 
 #         label="italic(TT[pd])", 
#         parse=TRUE, 
 #          size=24/.pt) +
   ylab("Displacement") + 
  ggtitle("Tongue tip") +
  theme_michigan2gpaper(n = 28, k = 24)

gA <- ggplotGrob(velum_trade_paper)
gB <- ggplotGrob(tongue_trade_tip_paper)
grid::grid.newpage()
fig.trade1.paper = grid.arrange(rbind(gA, gB))
fig.trade1.paper


ggsave(filename = file.path(pfadfigs,
                            "fig8.png"), 
     plot = fig.trade1.paper, 
       width = 25, 
      height = 20, 
       units = "cm")

```

## fig9 (Fig. 6.2 (tradeextra_paper.png))
```{r, fig.width=30/2.54, fig.height=30/2.54}
trade.df = tt.df %>%
  rename(Dialect = language) %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval = 
           velumopening_maxvel_off - velumopening_maxvel_on,
         peak_time = velumopening_maxcon_on - alveolarconstriction_maxvel_on,
         peak_prop = peak_time/intval, 
         displacement = alvUS_alveolarconstriction_maxcon_onset -
        alvUS_alveolarconstriction_gesture_onset) %>%
  ungroup() 

# remove 9 NAs
trade.df = trade.df %>%
  filter(!(is.na(displacement))) %>%
  filter(!(is.na(peak_prop))) 

dataBRE = trade.df %>%
  filter(Dialect == "BRE")
dataUSE = trade.df %>%
  filter(Dialect == "USE")

fig.tradeextra_paper = trade.df %>%
  ggplot + 
  aes(y = peak_prop, x = displacement, col = Dialect) +
  geom_point(size = 1.2) + 
facet_wrap(coda~vowel, ncol=3) +
  ylab(expression(Proportional~alignment~of~peak~velum~lowering~(italic(i[1]/i[2])))) +
  xlab("Magnitude of peak tongue tip raising") +
  geom_smooth(
    aes(y = peak_prop, x = displacement), 
    data = dataBRE, 
    method = "lm", color = "red") + 
  geom_smooth(
    aes(y = peak_prop, x = displacement), 
    data = dataUSE, 
    method = "lm", color = "blue") + 
    theme_michigan2epaper(position="top") + 
  theme(legend.title = element_blank()) +
  scale_colour_manual(values = cols) 
 fig.tradeextra_paper
 
 ggsave(filename = file.path(pfadfigs2, 
                             "fig9.png"), 
plot = fig.tradeextra_paper, 
width = 30, 
height = 30, 
units = "cm")
```

## fig10 (Fig. 6.3. (tradeall.paper.png))
```{r, fig.width=30/2.54, fig.height=30/2.54}
trade.lmer = trade.df %>%
  mutate(vowel_coda = interaction(vowel, coda)) %>%
  lmer(peak_prop ~ displacement * Dialect +
       (displacement  | vowel_coda) + 
         (displacement | speaker), . )

anova(trade.lmer)

emtrends(trade.lmer, pairwise ~ Dialect,  
         var = 'displacement')

emm.df = emmip(trade.lmer, Dialect ~ displacement, 
      cov.reduce = range, plotit=F)

g = emmip_ggplot(emm.df, Dialect ~ displacement, 
      cov.reduce = range)

fig.tradeall.paper = g + 
 theme_michigan10(n = 32) +
  theme(legend.title=element_blank())+
  geom_line(linewidth=1.5, lty=1) +
  ylab(expression(Proportional~alignment~of~peak~velum~lowering~(italic(i[1]/i[2])))) +
  xlab("Magnitude of peak tongue tip raising") +
scale_colour_manual(values = cols)
fig.tradeall.paper

ggsave(filename = file.path(pfadfigs2, "fig10.png"),  
       plot = fig.tradeall.paper, 
       width = 30, 
       height = 30, 
       units = "cm")
```

# Inter-gestural analysis of tongue tip and velum (section 7)
## fig11 (Fig. 7.1. fig.asynch.paper.png)

```{r, fig.width=45/2.54, fig.height=30/2.54}
j = 2.5
n = 54; k = 40
velum_asynch_paper = tt.df %>%
  filter(rowind == j) %>%
  ggplot +
  aes(y = velum2US, x = time_in_sequence) +
  geom_line(linewidth = lwd) +
  coord_cartesian(xlim = c(0, 0.5), ylim = c(0.18, 1.0)) +
  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 = .9, lty=2, lwd=lwd, col="black") +
   geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = velon.tt, y = .1, 
               xend = velon.tt, yend = .9, lty=2, lwd=lwd, col="slategray") +
  geom_segment(aes(y = velum2US, x = time_in_sequence), 
               x = ttmax.t, y = .1, 
               xend = ttmax.t, yend = .9, 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", 
                           length = unit(0.1, "inches"))) +
   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(.95, 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) %>%
  ggplot +
  aes(y = alvUS, x = time_in_sequence) +
  geom_line(col="slategray",  linewidth = lwd) +
  coord_cartesian(xlim = c(0, 0.5)) +
  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,
                            "fig11.png"), 
 plot = fig.asynch.paper, 
       width = 45, 
      height = 30, 
       units = "cm")

```

## fig12 (Fig. 7.2. (fig.velum_tongue.paper.png))
```{r}
x.df = tt.df %>%
  filter(coda !="nt") %>%
  mutate(reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  velumopening_maxcon_on), 5), 
         coda = factor(coda))

# remove 9 segments that have NA at displacement gesture onset
x.df = x.df %>% 
  filter(!is.na(alvUS_alveolarconstriction_gesture_onset)) 
x.df = x.df %>%
  mutate(alvUSn = alvUS - alvUS_alveolarconstriction_gesture_onset)

# get the mean relative times of TT peak opening velocitie

tt.maxr.t = x.df %>%
  rename(Dialect = language) %>%
  group_by(vowel, coda, Dialect) %>%
  summarise(m1 = mean(alveolarconstriction_maxvel_on -  
                       velumopening_maxcon_on, na.rm=T) * 1000, 
            m2 = mean(endv_time -  
                       velumopening_maxcon_on, na.rm=T) * 1000) %>%  ungroup()
```


```{r,fig.width=30/2.54, fig.height=30/2.54}
 cols = c("red", "blue")
 velum_tongue.paper = x.df %>%
group_by(language, coda, vowel, reltime)  %>%
    summarise(velum = mean(velum2US), tt = mean(alvUSn)) %>%
  ungroup() %>%
  rename(Dialect = language) %>%
    ggplot + 
    aes(y = tt, x = reltime, col = Dialect, 
        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(coda ~ vowel, ncol=3) +
    coord_cartesian(xlim = c(-200, 200)) +
    geom_vline(xintercept=0, lty=2) +
   theme_michigan2epaper(position="top") + 
   theme(legend.title = element_blank())+
   geom_vline(data = tt.maxr.t, aes(xintercept=m2,  col = Dialect),  lty=3) + 
   geom_vline(data = tt.maxr.t, aes(xintercept=m1,  col = Dialect)) + 
 #  ggtitle("Tongue tip height", 
 #          subtitle = "Aligned at peak velum opening") +
  scale_colour_manual(values = cols) 
 velum_tongue.paper
ggsave(filename = file.path(pfadfigs2,
                            "fig12.png"), 
     plot = velum_tongue.paper, 
       width = 30, 
      height = 30, 
       units = "cm")
```

## Statistics
### Interval i1: Interval between time of  maximum TT peak velocity and maximum velum opening

Statistics: both coda and dialect sig. with no interaction.
```{r}
#tt.df %>%
#  filter(coda %in% c("n", "nd", "nz")) %>%
#  filter(segment == "Vokal") %>%
#  group_by(rowind) %>%
#  slice_tail(n=1) %>%
# mutate(intval = 
 #          velumopening_maxcon_on - 
 #          alveolarconstriction_maxvel_on)     %>%
 # rename(Dialect = language) %>%
 # ungroup() %>%
#  lmer(intval ~ vowel * coda * Dialect +
#         (Dialect|stem) + (vowel+coda|speaker), .) %>%
#  step()

tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval = 
           velumopening_maxcon_on - 
           alveolarconstriction_maxvel_on)     %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  lmer(intval ~ coda + Dialect + 
         (1 | stem) + (1 | speaker), .) %>%
  anova()

# estimated marginal means for coda
tt.df.i1.coda.emmeans = tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval = 
           velumopening_maxcon_on - 
           alveolarconstriction_maxvel_on)     %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  lmer(intval ~ coda + Dialect + 
         (1 | stem) + (1 | speaker), .) %>%
  emmeans(.,  ~ coda) 

tt.df.i1.coda.emmeans

# significance of coda
tt.df.i1.coda.emmeans %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

# estimated marginal means for dialect
tt.df.i1.dialect.emmeans = tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(intval = 
           velumopening_maxcon_on - 
           alveolarconstriction_maxvel_on)     %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  lmer(intval ~ coda + Dialect + 
         (1 | stem) + (1 | speaker), .) %>%
  emmeans(.,  ~ Dialect) 
tt.df.i1.dialect.emmeans
```

### i2:  Interval between tongue tip maximum and velum opening maximum

Statistics: both vowel and coda, not dialect sig.
```{r}
#tt.df %>%
#  filter(coda %in% c("n", "nd", "nz")) %>%
#  filter(segment == "Vokal") %>%
#  group_by(rowind) %>%
#  slice_tail(n=1) %>%
#  mutate(peakintval = 
#           velumopening_maxcon_on - 
#           alveolarconstriction_maxcon_on)     %>%
#  rename(Dialect = language) %>%
#  ungroup() %>%
# lmer(peakintval ~ vowel * coda * Dialect + 
#       (Dialect|stem) + (vowel+coda|speaker), .) %>%
#  step()

tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(peakintval = 
           velumopening_maxcon_on - 
           alveolarconstriction_maxcon_on)     %>%
  rename(Dialect = language) %>%
  ungroup() %>%
lmer(peakintval ~ vowel + coda + 
       (coda | speaker), .) %>%
  anova()

# estimated marginal means for coda
tt.df.i2.coda.emmeans = 
tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(peakintval = 
           velumopening_maxcon_on - 
           alveolarconstriction_maxcon_on)     %>%
  rename(Dialect = language) %>%
  ungroup() %>%
lmer(peakintval ~ vowel + coda + 
       (coda | speaker), .) %>%
emmeans(.,  ~ coda) 
tt.df.i2.coda.emmeans

# significance for coda
tt.df.i2.coda.emmeans %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

# estimated marginal means for vowel
tt.df.i2.vowel.emmeans = 
tt.df %>%
  filter(coda %in% c("n", "nd", "nz")) %>%
  filter(segment == "Vokal") %>%
  group_by(rowind) %>%
  slice_tail(n=1) %>%
  mutate(peakintval = 
           velumopening_maxcon_on - 
           alveolarconstriction_maxcon_on)     %>%
  rename(Dialect = language) %>%
  ungroup() %>%
lmer(peakintval ~ vowel + coda + 
       (coda | speaker), .) %>%
emmeans(.,  ~ vowel) 
tt.df.i2.vowel.emmeans

# significance for vowel
tt.df.i2.vowel.emmeans %>%
pairs(., simple = "each", combine=T) %>%
  as.data.frame() %>%
  filter(p.value <= 0.05) %>%
  as_tibble() %>%
  mutate(p.value = round(p.value, 4))

```


# Discussion (section 7.3)

Schematic outline in final conclusion for the paper
```{r}
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) -.5
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]
times = c(blue.x, red.x, black.x)
y = c(blue.y, red.y, black.y)
Dialect = c(rep("USE",  length(blue.y)), 
            rep("BRE",  length(red.y)),
            rep("Both", length(black.y)))
scheme.df = data.frame(y, times, Dialect)
```

## fig13 (Fig. 7.3	(scheme.png))
```{r, fig.width=30/2.54, fig.height=30/2.54}
colscheme = c("black", "red", "blue")
ltyscheme = c(1, 2, 2)

schematic_paper = scheme.df %>%
  ggplot +
  aes(y = y, x = times, group=Dialect, col = Dialect,
      linetype = Dialect) +
   theme_michigan2hpaper(n = 30, k = 26, position="top") +
  geom_line(linewidth = 1.5) +
  ylab("Displacement") +
  xlab("Time") + 
  coord_cartesian(xlim = c(-15, 100)) +
   annotate("text", x=-5, y=2, 
           label="Velum", 
           size=28/.pt) +
  annotate("text", x=-5, y=-.5, 
           label="Tongue tip", 
           size=28/.pt) +
  geom_vline(xintercept = c(pktime.blue.time, pktime.red.time), 
              col=c("blue", "red")) +
  scale_colour_manual(values = colscheme) +
  scale_linetype_manual(values = ltyscheme)

schematic_paper

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

# Appendices

## Appendix C. fig.C.1 (Fig. C1		(figUSE_dial.png))
```{r}
a1 = dt.df %>%
  # change the USE dialect names
  mutate(Region =
           case_when(region4 == "US_Atlantic" ~ "Atlantic", 
                     region4 == "US_south"~ "South",
                     region4 == "US_midland" ~ "Midland", 
                     region4 == "US_west"~ "West",
                      TRUE ~ language), 
         Region2 = paste(Region, speaker), 
         # align at the acoustic vowel midpoint
         midv_time = 
    (endv_time+timeseriesoffset)/2, 
    reltime = 
           plyr::round_any(1000 * (time_in_sequence - midv_time), 5)) %>%
  # choose only fed words
   filter(coda %in% c("d") & vowel == c("ɛ")) %>%
  filter(language == "USE") %>%
  # filter between acoustic vowel onset and offset
  filter(time_in_sequence >= timeseriesoffset) %>%
  filter(time_in_sequence <= endv_time) %>%
  # aggregate
group_by(Region2, coda, vowel, reltime)  %>%
    summarise(pal = mean(palUS)) %>%
  ungroup() 
```

This is exactly the same as `a1` except for /æ/ before coda /n, nd, nz/
```{r}
a2 = dt.df %>%
  mutate(Region =
           case_when(region4 == "US_Atlantic" ~ "Atlantic", 
                     region4 == "US_south"~ "South",
                     region4 == "US_midland" ~ "Midland", 
                     region4 == "US_west"~ "West",
                      TRUE ~ language), 
         Region2 = paste(Region, speaker), 
         midv_time = 
    (endv_time+timeseriesoffset)/2, 
    reltime = 
           plyr::round_any(1000 * (time_in_sequence - midv_time), 5)) %>%
    filter(coda %in% c("n", "nd", "nz") & vowel == c("æ")) %>%
  filter(language == "USE") %>%
  filter(time_in_sequence >= timeseriesoffset) %>%
  filter(time_in_sequence <= endv_time) %>%
group_by(Region2, coda, vowel, reltime)  %>%
    summarise(pal = mean(palUS)) %>%
  ungroup() 
```

Bind `a1` and `a2` together and rename /æ/ as /æ̃/
```{r}
a3 = rbind(a1, a2)
a3 = a3 %>% 
  mutate(vowel = 
           case_when(vowel == "æ" ~ "æ̃", 
                     TRUE ~ vowel))
```

Aggregated plots by speaker separately for the two categories above
```{r, fig.width=30/2.54, fig.height=30/2.54}
cols.vowel = c("black", "orange")
USE_dialect = a3 %>%
    ggplot + 
    aes(y = pal, x = reltime, col = vowel, 
        group = interaction(Region2, coda, vowel)) +
  geom_line(linewidth=1.1) +
  theme_michigan6dial(n = 28) + 
  scale_x_continuous("Time (ms)", breaks = c(-200, 0, 200)) +
  facet_wrap(~Region2) +
    ylab("Tongue dorsum height") +
    coord_cartesian(xlim = c(-250, 250)) +
    geom_vline(xintercept=0, lty=2, linewidth=1) +
  scale_colour_manual(values = cols.vowel) 
USE_dialect
ggsave(filename = file.path(pfadfigs, "figC.1.png"), 
       plot = USE_dialect, 
       width = 30, 
       height = 30, 
       units = "cm")
```

## Appendix D. fig.D.1 (Fig. D1		(figappb.png))

```{r, fig.width=45/2.54, fig.height=30/2.54}
 appb.df = dt.df %>%
    filter(coda %in% c("n", "nd", "nz")) %>%
  mutate(reltime = 
           plyr::round_any(1000 * (time_in_sequence - timeseriesoffset), 5))

figappb =  appb.df %>%
group_by(language, coda, vowel, reltime)  %>%
    summarise(vel = mean(velum2US)) %>%
  ungroup() %>%
  rename(Dialect = language) %>%
    ggplot + 
    aes(y = vel, x = reltime, col = Dialect, 
        group = interaction(Dialect, coda, vowel)) +
  geom_line(linewidth=1.5) +
  scale_x_continuous("Time (ms)", breaks = c(-0, 50, 100)) +
    ylab("Velum height") +
    facet_wrap(coda ~ vowel, ncol=5) +
    coord_cartesian(xlim = c(-50, 150)) +
    geom_vline(xintercept=0, lty=2, linewidth=1.5) +
   theme_michigan6x() +
  scale_colour_manual(values = cols) 
 figappb
 
ggsave(filename = file.path(pfadfigs, "figD.1.png"), 
       plot = figappb, 
       width = 45, 
       height = 30, 
       units = "cm")
```


## Appendix E. fig.E.1 (Fig. E1		(figVN.png))
```{r,fig.width = 30/2.54, fig.height=20/2.54}

cols.dial = c("red", "blue")
vowel_duration = dt.df %>%
  group_by(rowind) %>%
  filter(segment == "Vokal") %>%
   filter(coda %in% c("n", "nd", "nz")) %>%
  slice_tail(n=1) %>%
  mutate(Vdur = endv_time - timeseriesoffset) %>%
  rename(Dialect = language) %>%
  ungroup() %>%
  ggplot +
  aes(y = Vdur, col = Dialect, x = coda) +
  geom_boxplot() +
  facet_wrap(~ vowel, ncol = 5) +
  theme_michigan2(n=32) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1), 
axis.title.y = element_text(size=28), 
legend.title = element_blank()) +
 #  scale_y_continuous("Proportion", breaks=c(-1, 0, 1)) +
  xlab("") +
ylab("V duration (s)") +
  # ggtitle("Proportional alignment of peak velum lowering") +
 # coord_cartesian(ylim = c(-1, 1)) + 
 # geom_hline(yintercept=0, lty=2) +
scale_colour_manual(values = cols.dial)

# duration from the acoustic vowel offset to time of peak velocity velum raising
N_duration =  dt.df %>%
  group_by(rowind) %>%
  filter(segment == "Vokal") %>%
   filter(coda %in% c("n", "nd", "nz")) %>%
  slice_tail(n=1) %>%
  mutate(Ndur = velumopening_maxvel_off - endv_time) %>%
  rename(Dialect = language) %>%
  ggplot +
  aes(y  = Ndur,   x = coda, col=Dialect) +
  geom_boxplot() +
  facet_wrap(~ vowel, ncol = 5) +
  ylab("N duration (s)") +
  xlab("") +
coord_cartesian(ylim = c(0, .6)) +
# ggtitle("Log. proportion  of nasalization in coda-/n/") +
theme_michigan3(n = 32) +
  theme(axis.text.y = 
element_text(angle = 90, vjust = 1, hjust=1), 
axis.title.y = element_text(size=28)) +
  scale_colour_manual(values = cols.dial) 

figVN = 
  grid.arrange(vowel_duration, N_duration, nrow=2)

ggsave(filename = file.path(pfadfigs, "figE.1.png"), 
      plot = figVN, 
       width = 30, 
       height = 20, 
      units = "cm")

```

## Appendix F
### fig.F.1 (Fig. F1		(velum_tongue2.png))

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

z.df = tt.df %>%
  filter(coda !="nt") %>%
  mutate(reltime = 
           plyr::round_any(1000 * 
                             (time_in_sequence -  velumopening_maxvel_off), 5), 
         coda = factor(coda))

# remove 9 segments that have NA at displacement gesture onset
z.df = z.df %>% 
  filter(!is.na(alvUS_alveolarconstriction_gesture_onset)) 
z.df = z.df %>%
  mutate(alvUSn = alvUS - alvUS_alveolarconstriction_gesture_onset)

# get the mean relative times of TT peak opening velocity

tt.maxr.t = z.df %>%
  rename(Dialect = language) %>%
  group_by(vowel, coda, Dialect) %>%
  summarise(m1 = mean(alveolarconstriction_maxcon_on -  
                       velumopening_maxvel_off, na.rm=T) * 1000, 
            m2 = mean(alveolarconstriction_maxvel_off -  
                       velumopening_maxvel_off, na.rm=T) * 1000) %>%  ungroup()


 cols.dial = c("red", "blue")
lwd = 1.2
velum_tongue2 = z.df %>%
group_by(language, coda, vowel, reltime)  %>%
    summarise(velum = mean(velum2US), tt = mean(alvUSn)) %>%
  ungroup() %>%
  rename(Dialect = language) %>%
    ggplot + 
    aes(y = tt, x = reltime, col = Dialect, 
        group = interaction(Dialect, coda, vowel)) +
     theme_michigan2epaper(position="top") +
  theme(legend.title = element_blank()) + 
  geom_line(linewidth=lwd) +
  scale_x_continuous("Time (ms)", breaks = c(-200, 0, 200)) +
    ylab("Tongue tip displacement") +
    facet_wrap(coda ~ vowel, ncol=3) +
    coord_cartesian(xlim = c(-250, 250)) +
    geom_vline(xintercept=0, lty=2, linewidth=lwd) +
 #  theme_michigan11(54, 36, position="top") + 
   geom_vline(data = tt.maxr.t, 
             aes(xintercept=m2,  col = Dialect),  
              lty=3,linewidth=lwd) + 
   geom_vline(data = tt.maxr.t, 
              aes(xintercept=m1,  col = Dialect),
              linewidth = lwd) + 
  scale_colour_manual(values = cols.dial) 
velum_tongue2
ggsave(filename = file.path(pfadfigs,
                            "figF.1.png"), 
     plot = velum_tongue2, 
       width = 30, 
      height = 30, 
       units = "cm")
```


### fig.F.2 (Fig. F2. veltongueint.png)
```{r, fig.width=30/2.54, fig.height=30/2.54}
cols.dial = c("red", "blue")
veltongueint = tt.df %>%
  filter(coda != "nt") %>%
  mutate(veltongueint =
           alveolarconstriction_maxcon_on -  
                       velumopening_maxvel_off) %>%
  rename(Dialect = language) %>%
  ggplot +
  aes(y = veltongueint,x = coda, col=Dialect) +
  geom_boxplot(lwd=1.4) +
  theme_michigan2(n = 32) +
  theme(legend.title = element_blank(), 
        axis.text.x = element_text(size=24)) +
  geom_hline(yintercept=0, lty=2, linewidth=1.1) + 
  facet_wrap(~ vowel, ncol = 5) +
  ylab("Duration (s)") +
  xlab("") +
scale_colour_manual(values = cols.dial)
veltongueint

ggsave(filename = file.path(pfadfigs, "figF.2.png"), 
   plot = veltongueint, 
      width = 30, 
      height = 20, 
      units = "cm")
```
