Vowel nasalisation in two varieties of English

Load libraries and data:

# 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"))

1 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)

2 Preliminaries

2.1 keep only monosyllables

df_erc_prenasalraising = df_erc_prenasalraising %>%
  filter(secvow == "_")

2.2 Create a segment identifier

This is just speaker and trial number pasted together

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

2.3 Change speaker id to a factor

df_erc_prenasalraising = df_erc_prenasalraising  %>%
  mutate(speaker = factor(paste0("S", speaker)))

2.4 remove S41

df_erc_prenasalraising = df_erc_prenasalraising  %>%
  filter(speaker!="S41") %>%
  mutate(speaker = factor(speaker))

2.5 Identify nasality and change vowel labels

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)

## 
##      n      o 
## 724957 333591

table(dt.df$vowel)

## 
##      æ     eɪ      ʌ      ɛ      ɪ 
## 214116 235598 190151 211110 207573

2.6 give orthographic labels

# 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))

## [1] FALSE

# 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)

## 
##    S1   S10   S11   S12   S13   S14   S15   S16   S17   S18    S2   S20   S21 
## 23883 21633 25510 25911 25595 25923 21346 24449 22565 18280 19369 22755 23681 
##   S22   S23   S24   S25   S26   S27   S28   S29    S3   S30   S31   S32   S33 
## 19433 24764 27599 28343 23141 22461 25327 28109 25393 25086 22805 26681 29153 
##   S34   S35   S36   S37   S38   S39    S4   S40   S42   S43   S44   S45    S5 
## 26516 23502 22725 24958 32093 28634 23480 27368 19567 27020 23442 22925 29230 
##    S6    S7    S8    S9 
## 27438 21233 23672 25550

# 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)

## 
##           BR US_northeast     US_south      US_west 
##       665276       178869        70089       144314

   table(dt.df$region4)

## 
##          BR US_Atlantic    US_south  US_midland     US_west 
##      665276       77525       70089      101344      144314

2.7 Read in kinematic markers for nasalisation

load(file.path(pfad, 
"erc_mrinasals_allsubjects_artaku_allsegandsig_table_v2_rs_addlang.Rda"))
dim(df_erc_prenasalraising)

## [1] 13577   156

# 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)

2.8 Read in kinematic markers for tongue tip

load(file.path(pfad, 
"erc_mrinasals_allsubjects_artaku_allsegandsig_table_v2_rs_addlang_addtt_sig.Rda"))
dim(df_erc_prenasalraising)

## [1] 13577   172

# 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) 

2.9 Acoustic times

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:

  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

##  [1] "14.259" "29.22"  "3.299"  "30.261" "34.131" "38.147" "38.216" "38.257"
##  [9] "39.136" "40.49"

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

2.10 data-frame of tongue tip data

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()

## [1] 1337

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)

2.11 Define themes

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))
}

3 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

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)

3.1 fig2 [Fig. 3.1. (fig1paper.png)].

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")

3.2 fig3 (Fig. 3.2. fig5paper.png)

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

Sum of nasalisation in coda

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()

## `summarise()` has grouped output by 'rowind', 'vowel', 'coda', 'language',
## 'onset'. You can override using the `.groups` argument.

Sum of nasalisation between points of maximum velocity.

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

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

Fig.3 for the paper

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

## TableGrob (2 x 1) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (2-2,1-1) arrange gtable[layout]

ggsave(filename = file.path(pfadfigs2, "fig3.png"), 
  plot = fig5paper, 
      width = 30, 
      height = 20, 
     units = "cm")

3.2.1 Statistics associated with fig2

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()

## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
## Dialect     0.684  0.6836     1   41.05  26.9030 6.121e-06 ***
## vowel      34.858  8.7145     4 1969.14 342.9443 < 2.2e-16 ***
## coda        0.813  0.4067     2 1969.10  16.0045 1.274e-07 ***
## vowel:coda  0.401  0.0502     8 1969.14   1.9749   0.04594 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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))

## # A tibble: 25 × 8
##    coda  vowel contrast estimate     SE    df t.ratio p.value
##    <chr> <chr> <chr>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>
##  1 n     .     æ - ʌ     -0.234  0.0186 1969.  -12.6   0     
##  2 n     .     æ - ɛ     -0.255  0.0185 1969.  -13.8   0     
##  3 n     .     æ - ɪ     -0.326  0.0186 1969.  -17.6   0     
##  4 n     .     eɪ - ʌ    -0.184  0.0172 1969.  -10.7   0     
##  5 n     .     eɪ - ɛ    -0.205  0.0171 1969.  -12.0   0     
##  6 n     .     eɪ - ɪ    -0.276  0.0172 1969.  -16.1   0     
##  7 n     .     ʌ - ɪ     -0.0927 0.0171 1969.   -5.41  0     
##  8 n     .     ɛ - ɪ     -0.0712 0.0171 1969.   -4.17  0.0014
##  9 nd    .     æ - eɪ    -0.0782 0.0211 1969.   -3.71  0.0094
## 10 nd    .     æ - ʌ     -0.308  0.0186 1969.  -16.6   0     
## 11 nd    .     æ - ɛ     -0.314  0.0172 1969.  -18.2   0     
## 12 nd    .     æ - ɪ     -0.397  0.0172 1969.  -23.1   0     
## 13 nd    .     eɪ - ʌ    -0.230  0.0222 1969.  -10.4   0     
## 14 nd    .     eɪ - ɛ    -0.236  0.0211 1969.  -11.2   0     
## 15 nd    .     eɪ - ɪ    -0.319  0.0211 1969.  -15.2   0     
## 16 nd    .     ʌ - ɪ     -0.0887 0.0186 1969.   -4.78  0.0001
## 17 nd    .     ɛ - ɪ     -0.0832 0.0172 1969.   -4.83  0.0001
## 18 nz    .     æ - eɪ    -0.0963 0.0242 1969.   -3.97  0.0034
## 19 nz    .     æ - ʌ     -0.276  0.0241 1969.  -11.5   0     
## 20 nz    .     æ - ɛ     -0.287  0.0242 1969.  -11.9   0     
## 21 nz    .     æ - ɪ     -0.340  0.0220 1969.  -15.4   0     
## 22 nz    .     eɪ - ʌ    -0.180  0.0243 1969.   -7.41  0     
## 23 nz    .     eɪ - ɛ    -0.191  0.0244 1969.   -7.82  0     
## 24 nz    .     eɪ - ɪ    -0.243  0.0222 1969.  -10.9   0     
## 25 .     æ     n - nd     0.0934 0.0186 1969.    5.02  0

EMM for dialect

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

##  Dialect emmean     SE   df lower.CL upper.CL
##  BRE      0.280 0.0183 41.3   0.2435    0.317
##  USE      0.125 0.0237 41.1   0.0773    0.173
## 
## Results are averaged over the levels of: vowel, coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM for vowel

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

## NOTE: Results may be misleading due to involvement in interactions

##  vowel  emmean     SE   df lower.CL upper.CL
##  æ     0.00527 0.0168 64.9  -0.0283   0.0388
##  eɪ    0.08011 0.0171 70.5   0.0459   0.1143
##  ʌ     0.27816 0.0168 64.9   0.2446   0.3117
##  ɛ     0.29068 0.0166 62.6   0.2574   0.3239
##  ɪ     0.35969 0.0163 57.7   0.3270   0.3923
## 
## Results are averaged over the levels of: Dialect, coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM for coda

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

## NOTE: Results may be misleading due to involvement in interactions

##  coda emmean     SE   df lower.CL upper.CL
##  n     0.226 0.0156 48.2    0.194    0.257
##  nd    0.179 0.0158 50.6    0.147    0.210
##  nz    0.204 0.0163 57.9    0.171    0.237
## 
## Results are averaged over the levels of: Dialect, vowel 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Statistics Fig. 2 bottom row

# Identify and remove NA (there are 2 of them)
rowind.rem = propcnas.df %>%
  filter(is.na(ratio)) %>%
  pull(rowind) %>% unique()
rowind.rem

## character(0)

propcnas.df = propcnas.df %>%
  filter(!(rowind %in% rowind.rem))

does not converge

#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ɪ

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

## Type III Analysis of Variance Table with Satterthwaite's method
##             Sum Sq Mean Sq NumDF   DenDF  F value    Pr(>F)    
## vowel      144.942  36.235     4  928.00 591.5946 < 2.2e-16 ***
## coda         2.684   1.342     2 1856.99  21.9107 3.934e-10 ***
## Dialect      1.760   1.760     1   41.02  28.7290 3.496e-06 ***
## vowel:coda   4.107   0.513     8 1853.76   8.3824 3.086e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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))

## # A tibble: 30 × 8
##    coda  vowel contrast estimate     SE    df t.ratio p.value
##    <chr> <chr> <chr>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>
##  1 n     .     æ - eɪ    -0.167  0.0296 1933.   -5.65  0     
##  2 n     .     æ - ʌ     -0.488  0.0295 1736.  -16.6   0     
##  3 n     .     æ - ɛ     -0.522  0.0296 1287.  -17.6   0     
##  4 n     .     æ - ɪ     -0.615  0.0293 1938.  -21.0   0     
##  5 n     .     eɪ - ʌ    -0.321  0.0278 1404.  -11.5   0     
##  6 n     .     eɪ - ɛ    -0.354  0.0279  967.  -12.7   0     
##  7 n     .     eɪ - ɪ    -0.448  0.0269 1956.  -16.6   0     
##  8 n     .     ʌ - ɪ     -0.127  0.0275 1405.   -4.63  0.0002
##  9 n     .     ɛ - ɪ     -0.0939 0.0275  958.   -3.41  0.0303
## 10 nd    .     æ - eɪ    -0.271  0.0334 1894.   -8.12  0     
## 11 nd    .     æ - ʌ     -0.751  0.0303 1027.  -24.8   0     
## 12 nd    .     æ - ɛ     -0.764  0.0268 1956.  -28.5   0     
## 13 nd    .     æ - ɪ     -0.865  0.0267 1956.  -32.4   0     
## 14 nd    .     eɪ - ʌ    -0.480  0.0353 1625.  -13.6   0     
## 15 nd    .     eɪ - ɛ    -0.493  0.0334 1894.  -14.8   0     
## 16 nd    .     eɪ - ɪ    -0.594  0.0334 1893.  -17.8   0     
## 17 nd    .     ʌ - ɪ     -0.114  0.0303 1022.   -3.77  0.0078
## 18 nd    .     ɛ - ɪ     -0.102  0.0267 1956.   -3.80  0.0066
## 19 nz    .     æ - eɪ    -0.294  0.0387 1879.   -7.58  0     
## 20 nz    .     æ - ʌ     -0.630  0.0395  971.  -15.9   0     
## 21 nz    .     æ - ɛ     -0.660  0.0375 1956.  -17.6   0     
## 22 nz    .     æ - ɪ     -0.734  0.0347 1920.  -21.2   0     
## 23 nz    .     eɪ - ʌ    -0.336  0.0398 1056.   -8.43  0     
## 24 nz    .     eɪ - ɛ    -0.366  0.0390 1877.   -9.40  0     
## 25 nz    .     eɪ - ɪ    -0.440  0.0359 1735.  -12.3   0     
## 26 .     æ     n - nd     0.256  0.0293 1940.    8.73  0     
## 27 .     æ     n - nz     0.123  0.0347 1956.    3.56  0.017 
## 28 .     æ     nd - nz   -0.133  0.0331 1909.   -4.00  0.0029
## 29 .     eɪ    n - nd     0.152  0.0336 1889.    4.53  0.0003
## 30 .     eɪ    nd - nz   -0.155  0.0380 1956.   -4.09  0.002

Estimated marginal means for dialect

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

##  Dialect emmean     SE   df lower.CL upper.CL
##  BRE     -0.372 0.0330 42.1   -0.439   -0.306
##  USE     -0.646 0.0419 43.5   -0.731   -0.562
## 
## Results are averaged over the levels of: vowel, coda 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Estimated marginal means for vowel

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

## NOTE: Results may be misleading due to involvement in interactions

##  vowel emmean     SE   df lower.CL upper.CL
##  æ     -0.960 0.0305 52.9   -1.021   -0.899
##  eɪ    -0.716 0.0311 56.9   -0.778   -0.654
##  ʌ     -0.337 0.0305 52.0   -0.398   -0.276
##  ɛ     -0.312 0.0300 52.5   -0.372   -0.252
##  ɪ     -0.222 0.0297 47.8   -0.282   -0.162
## 
## Results are averaged over the levels of: coda, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

Estimated marginal means for coda

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

## NOTE: Results may be misleading due to involvement in interactions

##  coda emmean     SE   df lower.CL upper.CL
##  n    -0.475 0.0285 43.1   -0.533   -0.418
##  nd   -0.559 0.0289 44.0   -0.618   -0.501
##  nz   -0.494 0.0297 49.2   -0.553   -0.434
## 
## Results are averaged over the levels of: vowel, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

4 Intra-gestural velum analysis (section 4)

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)

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

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)

## Warning: Removed 2 rows containing missing values (`geom_line()`).

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

fig3paper

## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]

ggsave(filename = file.path(pfadfigs, "fig4.png"), 
    plot = fig3paper, 
      width = 25, 
       height = 20, 
   units = "cm")

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

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) 

## `summarise()` has grouped output by 'language', 'coda', 'vowel'. You can
## override using the `.groups` argument.

 fig4paper

ggsave(filename = file.path(pfadfigs2, "fig5.png"), 
      plot = fig4paper, 
       width = 45, 
       height = 30, 
      units = "cm")

4.2.1 Statistics indirectly related to Fig. 4

4.2.1.1 Velum opening displacement maximum: Velumdmax

Displacement is averaged between the peak velocities of velum movement

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() 

## `summarise()` has grouped output by 'rowind', 'language', 'onset', 'coda',
## 'vowel'. You can override using the `.groups` argument.

#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()

## boundary (singular) fit: see help('isSingular')

## Type III Analysis of Variance Table with Satterthwaite's method
##                 Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## coda          0.107072 0.053536     2 52.461  8.9444 0.0004545 ***
## vowel         0.280118 0.070030     4 49.350 11.7000 9.265e-07 ***
## Dialect       0.002574 0.002574     1 40.969  0.4300 0.5156372    
## vowel:Dialect 0.083992 0.020998     4 51.463  3.5082 0.0131592 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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))

## boundary (singular) fit: see help('isSingular')

print(displ.em.df, n= nrow(displ.em.df))

## # A tibble: 8 × 8
##   vowel Dialect contrast estimate      SE    df t.ratio p.value
##   <chr> <chr>   <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 .     BRE     æ - eɪ     0.0866 0.0150   43.3    5.76  0     
## 2 .     BRE     æ - ɛ      0.0348 0.00965  42.3    3.60  0.0206
## 3 .     BRE     æ - ɪ      0.0646 0.0124   42.8    5.22  0.0001
## 4 .     BRE     eɪ - ʌ    -0.101  0.0153   42.4   -6.56  0     
## 5 .     BRE     eɪ - ɛ    -0.0519 0.0130   43.3   -3.99  0.0063
## 6 .     BRE     ʌ - ɛ      0.0487 0.00911  42.0    5.35  0.0001
## 7 .     BRE     ʌ - ɪ      0.0786 0.0121   42.0    6.51  0     
## 8 .     BRE     ɛ - ɪ      0.0299 0.00818  42.3    3.65  0.0179

EMM for vowel

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

## boundary (singular) fit: see help('isSingular')

## NOTE: Results may be misleading due to involvement in interactions

##  vowel emmean     SE   df lower.CL upper.CL
##  æ      0.701 0.0139 41.2    0.673    0.729
##  eɪ     0.650 0.0194 41.4    0.610    0.689
##  ʌ      0.719 0.0126 41.2    0.694    0.745
##  ɛ      0.689 0.0145 41.1    0.660    0.719
##  ɪ      0.660 0.0173 41.2    0.625    0.695
## 
## Results are averaged over the levels of: coda, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM for coda

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

## boundary (singular) fit: see help('isSingular')

##  coda emmean     SE   df lower.CL upper.CL
##  n     0.697 0.0142 41.0    0.669    0.726
##  nd    0.670 0.0154 41.8    0.639    0.701
##  nz    0.684 0.0153 41.7    0.653    0.714
## 
## Results are averaged over the levels of: vowel, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Test for coda differences: n > nd; no nz > nd, no {n, nz} differences

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))

## boundary (singular) fit: see help('isSingular')

print(displ.em.df, n= nrow(displ.em.df))

## # A tibble: 2 × 6
##   contrast estimate      SE    df t.ratio p.value
##   <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 n - nd     0.0270 0.00689  42.0    3.92   0.001
## 2 nd - nz   -0.0132 0.00494  41.9   -2.67   0.032

4.2.1.2 peak velocity of opening: Velumvmaxo

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()

## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF   DenDF F value    Pr(>F)    
## coda           22.287 11.1433     2 1940.78  9.8375 5.612e-05 ***
## vowel         111.871 27.9678     4 1576.02 24.6904 < 2.2e-16 ***
## Dialect         0.556  0.5557     1   41.01  0.4906    0.4876    
## vowel:Dialect  35.923  8.9807     4 1966.85  7.9283 2.448e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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))

## # A tibble: 10 × 8
##    vowel Dialect contrast estimate     SE    df t.ratio p.value
##    <chr> <chr>   <chr>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>
##  1 .     BRE     æ - eɪ      1.09  0.101  1970.   10.8   0     
##  2 .     BRE     æ - ɛ       0.676 0.0955 1946.    7.07  0     
##  3 .     BRE     æ - ɪ       0.542 0.0930 1971.    5.82  0     
##  4 .     BRE     eɪ - ʌ     -0.983 0.107  1626.   -9.17  0     
##  5 .     BRE     eɪ - ɛ     -0.417 0.100  1917.   -4.15  0.0009
##  6 .     BRE     eɪ - ɪ     -0.551 0.0970 1971.   -5.68  0     
##  7 .     BRE     ʌ - ɛ       0.566 0.104  1259.    5.47  0     
##  8 .     BRE     ʌ - ɪ       0.432 0.102  1381.    4.24  0.0006
##  9 .     USE     eɪ - ʌ     -0.525 0.134  1831.   -3.90  0.0024
## 10 .     USE     ʌ - ɪ       0.395 0.127  1686.    3.12  0.0463

EMM for vowel

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

## NOTE: Results may be misleading due to involvement in interactions

##  vowel emmean    SE   df lower.CL upper.CL
##  æ       6.39 0.194 42.2     6.00     6.78
##  eɪ      5.72 0.195 43.2     5.33     6.11
##  ʌ       6.47 0.195 42.0     6.08     6.87
##  ɛ       6.05 0.192 41.2     5.66     6.44
##  ɪ       6.06 0.192 40.7     5.67     6.45
## 
## Results are averaged over the levels of: coda, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM for coda

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

##  coda emmean    SE   df lower.CL upper.CL
##  n      5.99 0.188 37.9     5.61     6.37
##  nd     6.19 0.189 38.4     5.81     6.57
##  nz     6.24 0.192 40.4     5.85     6.62
## 
## Results are averaged over the levels of: vowel, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Test for coda differences: nd, nz > n

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))

## # A tibble: 2 × 6
##   contrast estimate     SE    df t.ratio p.value
##   <chr>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>
## 1 n - nd     -0.198 0.0559 1921.   -3.55  0.0012
## 2 n - nz     -0.245 0.0635 1936.   -3.86  0.0004

4.2.1.3 Articulatory duration between peak velocities: Velumduro

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()

## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## vowel      1.2553 0.31381     4 1099.8  59.1849 < 2.2e-16 ***
## coda       3.8772 1.93859     2 1904.7 365.6162 < 2.2e-16 ***
## vowel:coda 0.1175 0.01468     8 1886.7   2.7695  0.004803 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 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))

## # A tibble: 25 × 8
##    vowel coda  contrast estimate      SE    df t.ratio p.value
##    <chr> <chr> <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
##  1 æ     .     n - nd     0.0836 0.00860 1966.    9.72  0     
##  2 æ     .     n - nz     0.0653 0.0102  1969.    6.41  0     
##  3 eɪ    .     n - nd     0.0755 0.00983 1953.    7.68  0     
##  4 eɪ    .     n - nz     0.0706 0.00987 1954.    7.16  0     
##  5 ʌ     .     n - nd     0.125  0.00855 1967.   14.6   0     
##  6 ʌ     .     n - nz     0.0985 0.00985 1841.   10.0   0     
##  7 ɛ     .     n - nd     0.0946 0.00826  967.   11.5   0     
##  8 ɛ     .     n - nz     0.0886 0.00998 1405.    8.88  0     
##  9 ɪ     .     n - nd     0.106  0.00784 1964.   13.5   0     
## 10 ɪ     .     n - nz     0.0920 0.00852 1969.   10.8   0     
## 11 .     n     æ - eɪ     0.0333 0.00862 1964.    3.86  0.0052
## 12 .     n     æ - ɛ      0.0497 0.00880 1240.    5.65  0     
## 13 .     n     æ - ɪ      0.0583 0.00859 1966.    6.78  0     
## 14 .     n     ʌ - ɛ      0.0323 0.00832  771.    3.88  0.0052
## 15 .     n     ʌ - ɪ      0.0408 0.00814 1630.    5.01  0     
## 16 .     nd    æ - ʌ      0.0586 0.00899 1330.    6.52  0     
## 17 .     nd    æ - ɛ      0.0608 0.00786 1964.    7.73  0     
## 18 .     nd    æ - ɪ      0.0808 0.00785 1964.   10.3   0     
## 19 .     nd    eɪ - ɛ     0.0355 0.00982 1956.    3.62  0.0137
## 20 .     nd    eɪ - ɪ     0.0555 0.00981 1956.    5.66  0     
## 21 .     nz    æ - eɪ     0.0386 0.0114  1957.    3.40  0.0306
## 22 .     nz    æ - ʌ      0.0507 0.0117  1284.    4.32  0.0008
## 23 .     nz    æ - ɛ      0.0730 0.0110  1964.    6.61  0     
## 24 .     nz    æ - ɪ      0.0850 0.0102  1962.    8.33  0     
## 25 .     nz    eɪ - ɪ     0.0463 0.0105  1912.    4.40  0.0005

EMM means: vowel

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

## NOTE: Results may be misleading due to involvement in interactions

##  vowel emmean     SE   df lower.CL upper.CL
##  æ      0.339 0.0116 50.4    0.316    0.362
##  eɪ     0.307 0.0118 52.6    0.283    0.330
##  ʌ      0.297 0.0117 50.0    0.273    0.320
##  ɛ      0.278 0.0115 49.7    0.255    0.301
##  ɪ      0.264 0.0115 47.5    0.241    0.287
## 
## Results are averaged over the levels of: coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM means: coda

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

## NOTE: Results may be misleading due to involvement in interactions

##  coda emmean     SE   df lower.CL upper.CL
##  n     0.357 0.0112 44.1    0.334    0.379
##  nd    0.260 0.0113 44.9    0.237    0.283
##  nz    0.274 0.0114 47.9    0.251    0.297
## 
## Results are averaged over the levels of: vowel 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

5 Intra-gestural tongue tip analysis (section 5)

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()

5.1 fig6 (Fig.5.1. (fig.ttipaper.png))

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)

## Warning: Removed 2 rows containing missing values (`geom_line()`).

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")

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

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

## Warning: Removed 16 rows containing missing values (`geom_point()`).

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

## Warning: Removed 16 rows containing missing values (`geom_point()`).

5.2.1 Statistics

5.2.1.1 peak tongue tip displacement

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()

## Type III Analysis of Variance Table with Satterthwaite's method
##                Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## vowel         2.74704 1.37352     2  44.23 259.9266 < 2.2e-16 ***
## coda          0.45885 0.22943     2  47.78  43.4169 1.791e-11 ***
## Dialect       0.10989 0.10989     1  40.92  20.7951 4.573e-05 ***
## vowel:coda    0.10156 0.02539     4 882.95   4.8049 0.0007737 ***
## vowel:Dialect 0.03975 0.01987     2  40.46   3.7611 0.0317595 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

EMM means: dialect

tonguetip.df %>%
  lmer(displ ~ vowel + 
         coda + 
         Dialect + 
         (1 | onset) + 
         (vowel + coda | speaker) + 
         vowel:coda + vowel:Dialect, .) %>%
  emmeans(.,  ~ Dialect)

## NOTE: Results may be misleading due to involvement in interactions

##  Dialect emmean     SE   df lower.CL upper.CL
##  BRE      0.353 0.0197 11.2    0.309    0.396
##  USE      0.258 0.0226 17.1    0.210    0.305
## 
## Results are averaged over the levels of: vowel, coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM means: vowel

tonguetip.df %>%
  lmer(displ ~ vowel + 
         coda + 
         Dialect + 
         (1 | onset) + 
         (vowel + coda | speaker) + 
         vowel:coda + vowel:Dialect, .) %>%
  emmeans(.,  ~ vowel)

## NOTE: Results may be misleading due to involvement in interactions

##  vowel emmean     SE   df lower.CL upper.CL
##  ʌ      0.474 0.0224 17.6    0.427    0.521
##  ɛ      0.294 0.0201 12.3    0.250    0.337
##  ɪ      0.148 0.0172  6.9    0.107    0.189
## 
## Results are averaged over the levels of: coda, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM means: coda

tonguetip.df %>%
  lmer(displ ~ vowel + 
         coda + 
         Dialect + 
         (1 | onset) + 
         (vowel + coda | speaker) + 
         vowel:coda + vowel:Dialect, .) %>%
  emmeans(.,  ~ coda)

## NOTE: Results may be misleading due to involvement in interactions

##  coda emmean     SE   df lower.CL upper.CL
##  n     0.253 0.0190 10.2    0.211    0.295
##  nd    0.320 0.0190 10.2    0.278    0.363
##  nz    0.342 0.0192 10.5    0.299    0.384
## 
## Results are averaged over the levels of: vowel, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Post-hoc tests

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))

## # A tibble: 9 × 8
##   vowel Dialect contrast  estimate     SE    df t.ratio p.value
##   <chr> <chr>   <chr>        <dbl>  <dbl> <dbl>   <dbl>   <dbl>
## 1 ʌ     .       BRE - USE   0.127  0.0338  40.9    3.76  0.0048
## 2 ɛ     .       BRE - USE   0.0992 0.0286  41.1    3.47  0.0111
## 3 ɪ     .       BRE - USE   0.0592 0.0182  41.2    3.25  0.0208
## 4 .     BRE     ʌ - ɛ       0.194  0.0205  44.4    9.43  0     
## 5 .     BRE     ʌ - ɪ       0.360  0.0199  44.6   18.1   0     
## 6 .     BRE     ɛ - ɪ       0.166  0.0118  43.2   14.1   0     
## 7 .     USE     ʌ - ɛ       0.166  0.0266  43.1    6.24  0     
## 8 .     USE     ʌ - ɪ       0.292  0.0257  43.2   11.4   0     
## 9 .     USE     ɛ - ɪ       0.126  0.0155  43.7    8.13  0

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))

## # A tibble: 15 × 8
##    coda  vowel contrast estimate     SE    df t.ratio p.value
##    <chr> <chr> <chr>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>
##  1 n     .     ʌ - ɛ      0.176  0.0180  58.0    9.78  0     
##  2 n     .     ʌ - ɪ      0.297  0.0174  58.3   17.1   0     
##  3 n     .     ɛ - ɪ      0.121  0.0121 108.    10.0   0     
##  4 nd    .     ʌ - ɛ      0.194  0.0185  64.7   10.5   0     
##  5 nd    .     ʌ - ɪ      0.345  0.0180  66.8   19.2   0     
##  6 nd    .     ɛ - ɪ      0.151  0.0113  85.8   13.3   0     
##  7 nz    .     ʌ - ɛ      0.169  0.0212 110.     8.00  0     
##  8 nz    .     ʌ - ɪ      0.334  0.0203 108.    16.4   0     
##  9 nz    .     ɛ - ɪ      0.165  0.0133 155.    12.5   0     
## 10 .     ʌ     n - nd    -0.0895 0.0120 241.    -7.48  0     
## 11 .     ʌ     n - nz    -0.0990 0.0173 181.    -5.71  0     
## 12 .     ɛ     n - nd    -0.0711 0.0105 146.    -6.80  0     
## 13 .     ɛ     n - nz    -0.106  0.0147  96.8   -7.22  0     
## 14 .     ɪ     n - nd    -0.0416 0.0100 127.    -4.15  0.0011
## 15 .     ɪ     n - nz    -0.0621 0.0137  75.0   -4.53  0.0004

5.2.1.2 peak tongue tip velocity

#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()

## Type III Analysis of Variance Table with Satterthwaite's method
##        Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## vowel 2282.10 1141.05     2 1025.3 852.103 < 2.2e-16 ***
## coda   106.79   53.39     2 1029.8  39.873 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

EMM means: vowel

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

##  vowel emmean    SE   df lower.CL upper.CL
##  ʌ       6.21 0.231 17.2     5.73     6.70
##  ɛ       3.95 0.221 15.0     3.48     4.42
##  ɪ       2.16 0.223 15.1     1.68     2.63
## 
## Results are averaged over the levels of: coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM means: coda

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

##  coda emmean    SE   df lower.CL upper.CL
##  n      3.65 0.221 14.9     3.18     4.12
##  nd     4.23 0.224 15.3     3.75     4.71
##  nz     4.44 0.230 17.2     3.95     4.92
## 
## Results are averaged over the levels of: vowel 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

Post-hoc tests

# 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))

## # A tibble: 3 × 6
##   contrast estimate     SE    df t.ratio p.value
##   <chr>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>
## 1 ʌ - ɛ        2.27 0.0983 1028.    23.0       0
## 2 ʌ - ɪ        4.06 0.0991 1041.    40.9       0
## 3 ɛ - ɪ        1.79 0.0840 1013.    21.3       0

# 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))

## # A tibble: 2 × 6
##   contrast estimate     SE    df t.ratio p.value
##   <chr>       <dbl>  <dbl> <dbl>   <dbl>   <dbl>
## 1 n - nd     -0.578 0.0831 1022.   -6.95       0
## 2 n - nz     -0.785 0.0987 1032.   -7.95       0

5.2.1.3 duration between tongue tip peak velocities

alveolarconstriction_maxvel_dur

#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()

## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## coda    0.57360 0.286800     2 43.245  85.743 8.948e-16 ***
## Dialect 0.05527 0.055265     1 41.169  16.522 0.0002112 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

EMM mean dialect:

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

##  Dialect emmean     SE   df lower.CL upper.CL
##  BRE      0.296 0.0129 33.7    0.270    0.322
##  USE      0.226 0.0158 42.8    0.194    0.258
## 
## Results are averaged over the levels of: coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

EMM means coda

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

##  coda emmean     SE   df lower.CL upper.CL
##  n     0.199 0.0109 23.7    0.176    0.221
##  nd    0.254 0.0113 25.6    0.230    0.277
##  nz    0.331 0.0151 41.0    0.301    0.362
## 
## Results are averaged over the levels of: Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

5.2.1.4 Interval between time of maximum peak velocity and peak displacement

#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()

## Type III Analysis of Variance Table with Satterthwaite's method
##           Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## vowel   0.039943 0.019972     2  7.113  11.045  0.006584 ** 
## coda    0.058600 0.029300     2 44.420  16.203 5.194e-06 ***
## Dialect 0.036529 0.036529     1 40.044  20.201 5.822e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Vowel, coda, dialect all sig in the above.

Analysis of vowel and coda coda: n < /nd, nz/. vowel: ʌ > ɪ

# 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

##  coda emmean      SE   df lower.CL upper.CL
##  n    0.0934 0.00600 45.7   0.0813    0.105
##  nd   0.1165 0.00529 44.4   0.1059    0.127
##  nz   0.1243 0.00810 47.1   0.1080    0.141
## 
## Results are averaged over the levels of: vowel, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

# 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))

## # A tibble: 2 × 6
##   contrast estimate      SE    df t.ratio p.value
##   <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 n - nd    -0.0232 0.00420  44.7   -5.51  0     
## 2 n - nz    -0.0309 0.00769  43.8   -4.02  0.0007

# 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 

##  vowel emmean      SE   df lower.CL upper.CL
##  ʌ     0.1242 0.00671 33.4   0.1105    0.138
##  ɛ     0.1122 0.00608 41.5   0.0999    0.124
##  ɪ     0.0978 0.00613 36.0   0.0854    0.110
## 
## Results are averaged over the levels of: coda, Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

# 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))

## # A tibble: 1 × 6
##   contrast estimate      SE    df t.ratio p.value
##   <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 ʌ - ɪ      0.0263 0.00573  8.41    4.59  0.0047

# 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 

##  Dialect emmean      SE   df lower.CL upper.CL
##  BRE     0.1330 0.00650 47.4   0.1199    0.146
##  USE     0.0898 0.00826 47.1   0.0732    0.106
## 
## Results are averaged over the levels of: vowel, coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

6 Vowel nasalization and tongue tip lenition (section 6)

# 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()

6.1 fig8 (fig.trade1.paper.png)

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

## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]

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

6.2 fig9 (Fig. 6.2 (tradeextra_paper.png))

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

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

 ggsave(filename = file.path(pfadfigs2, 
                             "fig9.png"), 
plot = fig.tradeextra_paper, 
width = 30, 
height = 30, 
units = "cm")

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

6.3 fig10 (Fig. 6.3. (tradeall.paper.png))

trade.lmer = trade.df %>%
  mutate(vowel_coda = interaction(vowel, coda)) %>%
  lmer(peak_prop ~ displacement * Dialect +
       (displacement  | vowel_coda) + 
         (displacement | speaker), . )

anova(trade.lmer)

## Type III Analysis of Variance Table with Satterthwaite's method
##                       Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## displacement         0.51780 0.51780     1 10.521 29.1014 0.0002549 ***
## Dialect              0.31896 0.31896     1 37.442 17.9261 0.0001432 ***
## displacement:Dialect 0.05293 0.05293     1 49.149  2.9748 0.0908545 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

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

## $emtrends
##  Dialect displacement.trend     SE   df lower.CL upper.CL
##  BRE                  0.325 0.0736 13.6    0.166    0.483
##  USE                  0.440 0.0850 24.2    0.264    0.615
## 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast  estimate     SE   df t.ratio p.value
##  BRE - USE   -0.115 0.0676 45.8  -1.702  0.0955
## 
## Degrees-of-freedom method: kenward-roger

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")

7 Inter-gestural analysis of tongue tip and velum (section 7)

7.1 fig11 (Fig. 7.1. fig.asynch.paper.png)

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")

7.2 fig12 (Fig. 7.2. (fig.velum_tongue.paper.png))

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()

## `summarise()` has grouped output by 'vowel', 'coda'. You can override using the
## `.groups` argument.

 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) 

## `summarise()` has grouped output by 'language', 'coda', 'vowel'. You can
## override using the `.groups` argument.

 velum_tongue.paper

ggsave(filename = file.path(pfadfigs2,
                            "fig12.png"), 
     plot = velum_tongue.paper, 
       width = 30, 
      height = 30, 
       units = "cm")

7.3 Statistics

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

Statistics: both coda and dialect sig. with no interaction.

#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()

## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## coda    0.54528 0.272642     2 752.71  127.30 < 2.2e-16 ***
## Dialect 0.02780 0.027801     1  40.72   12.98  0.000848 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 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

##  coda emmean      SE   df lower.CL upper.CL
##  n    0.1218 0.00658 50.7   0.1086   0.1351
##  nd   0.0724 0.00662 51.6   0.0591   0.0857
##  nz   0.0798 0.00697 61.1   0.0658   0.0937
## 
## Results are averaged over the levels of: Dialect 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

# 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))

## # A tibble: 2 × 6
##   contrast estimate      SE    df t.ratio p.value
##   <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 n - nd     0.0495 0.00329  808.    15.0       0
## 2 n - nz     0.0421 0.00390  698.    10.8       0

# 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

##  Dialect emmean      SE   df lower.CL upper.CL
##  BRE     0.1135 0.00769 43.7   0.0980   0.1290
##  USE     0.0691 0.00990 43.0   0.0492   0.0891
## 
## Results are averaged over the levels of: coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

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

Statistics: both vowel and coda, not dialect sig.

#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()

## Type III Analysis of Variance Table with Satterthwaite's method
##        Sum Sq  Mean Sq NumDF  DenDF F value    Pr(>F)    
## vowel 0.21985 0.109927     2 972.14  47.728 < 2.2e-16 ***
## coda  0.19078 0.095391     2  42.00  41.417 1.161e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# 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

##  coda  emmean      SE   df lower.CL upper.CL
##  n     0.0288 0.00715 42.0   0.0144   0.0433
##  nd   -0.0450 0.00399 43.1  -0.0530  -0.0369
##  nz   -0.0457 0.00670 42.8  -0.0592  -0.0322
## 
## Results are averaged over the levels of: vowel 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

# 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))

## # A tibble: 2 × 6
##   contrast estimate      SE    df t.ratio p.value
##   <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 n - nd     0.0738 0.00814  42.1    9.07       0
## 2 n - nz     0.0746 0.00997  42.2    7.48       0

# 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

##  vowel  emmean      SE   df lower.CL upper.CL
##  ʌ     -0.0359 0.00467 97.3 -0.04512 -0.02660
##  ɛ     -0.0250 0.00421 65.2 -0.03342 -0.01661
##  ɪ     -0.0010 0.00416 62.1 -0.00932  0.00731
## 
## Results are averaged over the levels of: coda 
## Degrees-of-freedom method: kenward-roger 
## Confidence level used: 0.95

# 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))

## # A tibble: 3 × 6
##   contrast estimate      SE    df t.ratio p.value
##   <chr>       <dbl>   <dbl> <dbl>   <dbl>   <dbl>
## 1 ʌ - ɛ     -0.0108 0.00388  973.   -2.80  0.0158
## 2 ʌ - ɪ     -0.0349 0.00387  973.   -9.01  0     
## 3 ɛ - ɪ     -0.0240 0.00331  972.   -7.25  0

8 Discussion (section 7.3)

Schematic outline in final conclusion for the paper

library(emuR)

## 
## Attaching package: 'emuR'

## The following object is masked from 'package:Matrix':
## 
##     norm

## The following object is masked from 'package:base':
## 
##     norm

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)

8.1 fig13 (Fig. 7.3 (scheme.png))

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")

9 Appendices

9.1 Appendix C. fig.C.1 (Fig. C1 (figUSE_dial.png))

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() 

## `summarise()` has grouped output by 'Region2', 'coda', 'vowel'. You can
## override using the `.groups` argument.

This is exactly the same as a1 except for /æ/ before coda /n, nd, nz/

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() 

## `summarise()` has grouped output by 'Region2', 'coda', 'vowel'. You can
## override using the `.groups` argument.

Bind a1 and a2 together and rename /æ/ as /æ̃/

a3 = rbind(a1, a2)
a3 = a3 %>% 
  mutate(vowel = 
           case_when(vowel == "æ" ~ "æ̃", 
                     TRUE ~ vowel))

Aggregated plots by speaker separately for the two categories above

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")

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

 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) 

## `summarise()` has grouped output by 'language', 'coda', 'vowel'. You can
## override using the `.groups` argument.

 figappb

ggsave(filename = file.path(pfadfigs, "figD.1.png"), 
       plot = figappb, 
       width = 45, 
       height = 30, 
       units = "cm")

9.3 Appendix E. fig.E.1 (Fig. E1 (figVN.png))

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")

9.4 Appendix F

9.4.1 fig.F.1 (Fig. F1 (velum_tongue2.png))

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()

## `summarise()` has grouped output by 'vowel', 'coda'. You can override using the
## `.groups` argument.

 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) 

## `summarise()` has grouped output by 'language', 'coda', 'vowel'. You can
## override using the `.groups` argument.

velum_tongue2

ggsave(filename = file.path(pfadfigs,
                            "figF.1.png"), 
     plot = velum_tongue2, 
       width = 30, 
      height = 30, 
       units = "cm")

9.4.2 fig.F.2 (Fig. F2. veltongueint.png)

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")